void aligntestmacro() 
{
  //Constants-------------------------------------
  const Double_t meter=1000.;  // 1 m = 1000 mm
  const Double_t cm=10.;       // 1 * cm = 10 * mm
  const Double_t mm=1.;        // 1 mm is default unit
  const Double_t micron=1e-3;  // 1 micron = 10^-3 mm

  const Double_t rad = 1.;  // 1 rad is default angular unit
  const Double_t mrad = 1e-3;  // 1 * millirad = 10^-3 rad
  const Double_t degree=TMath::Pi()/180.;

  const Int_t NumberOfStrips = 770; // Number of strips on one readout side plus 2 (Strip 0 and 769 are "guard ring" so to speak)
  const Double_t StripPitch = 80 * micron;
  const Double_t ModuleHalfLenghtX = NumberOfStrips * StripPitch * 0.5;  // half the length of one SCT wafer in x direction (perpendicular to strips)
  const Double_t ModuleHalfLenghtY = 6.0 * cm;  // half the length of the SCT wafer in y direction (along to strips)
  const Double_t ModuleHalfLenghtZ = 142.5 * micron;  // half the length of the SCT wafer in z direction (normal to wafer plane)
  
  const Double_t InitialSCT_x = 8.8*cm;
  const Double_t InitialSCT_y = 20*cm;
  const Double_t InitialSCT_z =-50*cm;

  TGeoVolume *sct_wafer;
  TGeoVolume *sct;
  TGeoVolume *top;

  TGeoPhysicalNode* node; //physical node for shifting module ("alignment")

  TGeoCombiTrans *ph; //matrix for sct module position
  
  TGeoCombiTrans *combi1;
  TGeoCombiTrans *combi2;
  
  delete gRandom;
  gRandom=new TRandom3(196851);

  // Geometry
  new TGeoManager("SCT", "Residual");

  TGeoMaterial *SiMat = new TGeoMaterial("Si",28.086,14,2.42);
  TGeoMaterial *VacMat = new TGeoMaterial("Vac",0,0,0);
  TGeoMedium *Si = new TGeoMedium ("Si",1,SiMat);
  TGeoMedium *Vac = new TGeoMedium("Vac",1,VacMat);

  //Volumen aus Shape und Medium basteln
  top = gGeoManager->MakeBox("TOP",Vac,2 * meter, 2 * meter, 2 * meter);
  sct =gGeoManager->MakeBox("SCT",Vac,10 * cm,10 * cm,1 * cm);
  //sct->SetVisibility(kFALSE);
  sct_wafer = gGeoManager->MakeTrap("SCT_Wafer", Si, ModuleHalfLenghtZ, 0, 0, ModuleHalfLenghtY ,ModuleHalfLenghtX ,ModuleHalfLenghtX,
                                                                                       0, ModuleHalfLenghtY ,ModuleHalfLenghtX ,ModuleHalfLenghtX,0);
  sct_wafer->SetLineColor(4);

  TGeoRotation *r_phi  = new TGeoRotation("rot1", 90 , 0, 90, 90 , 0, 0);
  TGeoRotation *stereo = new TGeoRotation("rot2", 90 ,  (40 * mrad) / degree, 90, 90 + (40 * mrad) / degree, 0, 0);

  combi1 = new TGeoCombiTrans(0, 0,  482.5 * micron, r_phi);
  combi2 = new TGeoCombiTrans(0, 0, -482.5 * micron, stereo);

  sct->AddNode(sct_wafer,1,combi1);
  sct->AddNode(sct_wafer,2,combi2);
  top->AddNode(sct,1); // This is new

  // set global graphics environment variables
  gGeoManager->SetTopVolume(top);
  gGeoManager->CloseGeometry();
	
  node = gGeoManager->MakePhysicalNode("/TOP_1/SCT_1");

  //--------------------------------------
  top->Draw(); //comment out this line
  //--------------------------------------
  
  TVector3 d1(InitialSCT_x,InitialSCT_y,InitialSCT_z);
  TVector3 dir_xy(InitialSCT_x,InitialSCT_y,0);
  TVector3 dir_xz(InitialSCT_x,0,InitialSCT_z);  
  
  //TVector3 ParticleDirection = d1.Unit();

  TVector3 ParticlePosition(0.,0.,0.);
  TVector3 dir;
  
  // Determine solid angle direction cone
  const Double_t Phi_0 = d1.Phi();
  const Double_t Phi_Plus = atan(ModuleHalfLenghtY / dir_xy.Mag());
  const Double_t PhiFrom = Phi_0 - Phi_Plus;
  const Double_t PhiTo = Phi_0 + Phi_Plus;  
  
  const Double_t Theta_0 =  acos(d1.CosTheta());
  const Double_t Theta_Plus = atan(ModuleHalfLenghtY / dir_xz.Mag());
  const Double_t ThetaFrom = Theta_0 - Theta_Plus;
  const Double_t ThetaTo = Theta_0 + Theta_Plus;  
 
  Double_t r,Phi,Theta; //coordinates for random direction
  
  for ( int k = 0; k < 10; k++)
  {
	  cout << "Track No. " << k << endl;

    //determine particle direction
    //equally distributed directions over solid angle of SCT module
	  r=1.;
	  Phi= gRandom->Rndm() * (PhiTo-PhiFrom) + PhiFrom;
	  Theta = acos( gRandom->Rndm() * (cos(ThetaFrom) - cos(ThetaTo)) + cos(ThetaTo) );		
	  dir.SetXYZ(r*sin(Theta)*cos(Phi),r*sin(Theta)*sin(Phi),r*cos(Theta));

	  GetHitPositionsLocal(InitialSCT_x, InitialSCT_y, InitialSCT_z, //position of module
			       0.,0.,0.,  //angles
			       ParticlePosition.X(), ParticlePosition.Y(), ParticlePosition.Z(),
			       dir.X(), dir.Y(), dir.Z());

  }
}


//__________________________________________________________________________
void GetHitPositionsLocal(Double_t SCT_Pos_x, Double_t SCT_Pos_y, Double_t SCT_Pos_z,
                          Double_t alpha, Double_t beta, Double_t gamma,
                          Double_t ParticlePos_x, Double_t ParticlePos_y, Double_t ParticlePos_z,
                          Double_t ParticleDir_x, Double_t ParticleDir_y, Double_t ParticleDir_z)
{
  Double_t *HitPoint;
  Double_t HitPointLocal[3];
  
  TGeoPhysicalNode *node=gGeoManager->GetPhysicalNode(0);
  //const char* nodename=node->GetName();
  //cout<<nodename<<endl;

  if (node->IsAligned())
  {
		ph = (TGeoCombiTrans*)node->GetNode()->GetMatrix();
  }	else {
		ph = new TGeoCombiTrans();
  }

  ph->Clear();
  ph->RotateZ(gamma);
  ph->RotateY(beta);
  ph->RotateX(alpha);
  ph->SetTranslation(SCT_Pos_x,SCT_Pos_y,SCT_Pos_z);
  	
  node->Align(ph);
	
  gGeoManager->InitTrack(ParticlePos_x, ParticlePos_y, ParticlePos_z, ParticleDir_x, ParticleDir_y, ParticleDir_z);

  TGeoNode *NextBoundary = gGeoManager->FindNextBoundary();
  if (NextBoundary)
  {
    cout<<"GetHitPosition 1st boundary:"<<endl;
    NextBoundary->Print();
    gGeoManager->Step(kTRUE,kTRUE);
    HitPoint = gGeoManager->GetCurrentPoint();
    cout<<"Hitpoint global: "<<HitPoint[0]<<", "<<HitPoint[1]<<", "<<HitPoint[2]<<endl; //should be around 88, 200, -500
  }
}