STEER/AliLego.cxx

   1 //////////////////////////////////////////////////////////////
   2 //////////////////////////////////////////////////////////////
   3 //
   4 //  Utility class to evaluate the material budget from
   5 //  a given radius to the surface of an arbitrary cylinder
   6 //  along radial directions from the centre:
   7 //
   8 //   - radiation length
   9 //   - Interaction length
  10 //   - g/cm2
  11 //
  12 //  Geantinos are shot in the bins in the fNtheta bins in theta
  13 //  and fNphi bins in phi with specified rectangular limits.
  14 //  The statistics are accumulated per
  15 //    fRadMin < r < fRadMax    and  <0 < z < fZMax
  16 //
  17 //  To activate this option, you can do:
  18 //      Root > gAlice.RunLego();
  19 //  or  Root > .x menu.C  then select menu item "RunLego"
  20 //  Note that when calling gAlice->RunLego, an optional list
  21 //  of arguments may be specified.
  22 //
  23 //Begin_Html
  24 /*
  25 <img src="picts/alilego.gif">
  26 */
  27 //End_Html
  28 //
  29 //////////////////////////////////////////////////////////////
  30
  31 #include "TMath.h"
  32 #include "AliLego.h"
  33 #include "AliRun.h"
  34 #include "AliConst.h"
  35
  36 ClassImp(AliLego)
  37
  38
  39 //___________________________________________
  40 AliLego::AliLego()
  41 {
  42    fHistRadl = 0;
  43    fHistAbso = 0;
  44    fHistGcm2 = 0;
  45    fHistReta = 0;
  46 }
  47
  48 //___________________________________________
  49 AliLego::AliLego(const char *name, const char *title)
  50         : TNamed(name,title)
  51 {
  52    fHistRadl = 0;
  53    fHistAbso = 0;
  54    fHistGcm2 = 0;
  55    fHistReta = 0;
  56 }
  57
  58 //___________________________________________
  59 AliLego::~AliLego()
  60 {
  61    delete fHistRadl;
  62    delete fHistAbso;
  63    delete fHistGcm2;
  64    delete fHistReta;
  65 }
  66
  67 //___________________________________________
  68 void AliLego::GenerateKinematics()
  69 {
  70 // Create a geantino with kinematics corresponding to the current
  71 // bins in theta and phi.
  72
  73   //
  74   // Rootinos are 0
  75    const Int_t mpart = 0;
  76    Float_t orig[3], pmom[3];
  77    Float_t t, cost, sint, cosp, sinp;
  78
  79 // --- Set to 0 radiation length, absorption length and g/cm2 ---
  80    fTotRadl = 0;
  81    fTotAbso = 0;
  82    fTotGcm2 = 0;
  83
  84    fCurTheta = (fThetaMin+(fThetaBin-0.5)*(fThetaMax-fThetaMin)/fNtheta)*kDegrad;
  85    fCurPhi   = (fPhiMin+(fPhiBin-0.5)*(fPhiMax-fPhiMin)/fNphi)*kDegrad;
  86    cost      = TMath::Cos(fCurTheta);
  87    sint      = TMath::Sin(fCurTheta);
  88    cosp      = TMath::Cos(fCurPhi);
  89    sinp      = TMath::Sin(fCurPhi);
  90
  91    pmom[0] = cosp*sint;
  92    pmom[1] = sinp*sint;
  93    pmom[2] = cost;
  94
  95 // --- Where to start
  96    orig[0] = orig[1] = orig[2] = 0;
  97    Float_t dalicz = 3000;
  98    if (fRadMin > 0) {
  99       t = PropagateCylinder(orig,pmom,fRadMin,dalicz);
 100       orig[0] = pmom[0]*t;
 101       orig[1] = pmom[1]*t;
 102       orig[2] = pmom[2]*t;
 103       if (TMath::Abs(orig[2]) > fZMax) return;
 104    }
 105
 106 // --- We do start here
 107    Float_t polar[3]={0.,0.,0.};
 108    Int_t ntr;
 109    gAlice->SetTrack(1, 0, mpart, pmom, orig, polar, 0, "LEGO ray", ntr);
 110 }
 111
 112 //___________________________________________
 113 void AliLego::Init(Int_t ntheta,Float_t themin,Float_t themax,
 114           Int_t nphi,Float_t phimin,Float_t phimax,Float_t rmin,Float_t rmax,
 115                   Float_t zmax)
 116 {
 117 // specify the angular limits and the size of the rectangular box
 118    fNtheta   = ntheta;
 119    fThetaMin = themin;
 120    fThetaMax = themax;
 121    fNphi     = nphi;
 122    fPhiMin   = phimin;
 123    fPhiMax   = phimax;
 124    fRadMin   = rmin;
 125    fRadMax   = rmax;
 126    fZMax     = zmax;
 127    Float_t etamin = -TMath::Log(TMath::Tan(TMath::Min((Double_t)fThetaMax*kDegrad/2,TMath::Pi()/2-1.e-10)));
 128    Float_t etamax = -TMath::Log(TMath::Tan(TMath::Max((Double_t)fThetaMin*kDegrad/2,              1.e-10)));
 129
 130    fHistRadl = new TH2F("hradl","Radiation length map",    nphi,phimin,phimax,ntheta,themin,themax);
 131    fHistAbso = new TH2F("habso","Interaction length map",  nphi,phimin,phimax,ntheta,themin,themax);
 132    fHistGcm2 = new TH2F("hgcm2","g/cm2 length map",        nphi,phimin,phimax,ntheta,themin,themax);
 133    fHistReta = new TH2F("hetar","Radiation length vs. eta",nphi,phimin,phimax,ntheta,etamin,etamax);
 134
 135 }
 136
 137 //___________________________________________
 138 Float_t AliLego::PropagateCylinder(Float_t *x, Float_t *v, Float_t r, Float_t z)
 139 {
 140 // Propagate to cylinder from inside
 141
 142    Double_t hnorm, sz, t, t1, t2, t3, sr;
 143    Double_t d[3];
 144    const Float_t kSmall  = 1e-8;
 145    const Float_t kSmall2 = kSmall*kSmall;
 146
 147 // ---> Find intesection with Z planes
 148    d[0]  = v[0];
 149    d[1]  = v[1];
 150    d[2]  = v[2];
 151    hnorm = TMath::Sqrt(1/(d[0]*d[0]+d[1]*d[1]+d[2]*d[2]));
 152    d[0] *= hnorm;
 153    d[1] *= hnorm;
 154    d[2] *= hnorm;
 155    if (d[2] > kSmall)       sz = (z-x[2])/d[2];
 156    else if (d[2] < -kSmall) sz = -(z+x[2])/d[2];
 157    else                     sz = 1.e10;  // ---> Direction parallel to X-Y, no intersection
 158
 159 // ---> Intersection with cylinders
 160 //      Intersection point (x,y,z)
 161 //      (x,y,z) is on track :    x=X(1)+t*D(1)
 162 //                               y=X(2)+t*D(2)
 163 //                               z=X(3)+t*D(3)
 164 //      (x,y,z) is on cylinder : x**2 + y**2 = R**2
 165 //
 166 //      (D(1)**2+D(2)**2)*t**2
 167 //      +2.*(X(1)*D(1)+X(2)*D(2))*t
 168 //      +X(1)**2+X(2)**2-R**2=0
 169 // ---> Solve second degree equation
 170    t1 = d[0]*d[0] + d[1]*d[1];
 171    if (t1 <= kSmall2) {
 172       t = sz;  // ---> Track parallel to the z-axis, take distance to planes
 173    } else {
 174       t2 = x[0]*d[0] + x[1]*d[1];
 175       t3 = x[0]*x[0] + x[1]*x[1];
 176       // ---> It should be positive, but there may be numerical problems
 177       sr = (t2 +TMath::Sqrt(TMath::Max(t2*t2-(t3-r*r)*t1,0.)))/t1;
 178       // ---> Find minimum distance between planes and cylinder
 179       t  = TMath::Min(sz,sr);
 180    }
 181    return t;
 182 }
 183
 184 //___________________________________________
 185 void AliLego::Run()
 186 {
 187    // loop on phi,theta bins
 188    gMC->InitLego();
 189    Float_t thed, phid, eta;
 190    for (fPhiBin=1; fPhiBin<=fNphi; fPhiBin++) {
 191       printf("AliLego Generating rays in phi bin:%d\n",fPhiBin);
 192       for (fThetaBin=1; fThetaBin<=fNtheta; fThetaBin++) {
 193          gMC->Gtrigi();
 194          gMC->Gtrigc();
 195          GenerateKinematics();
 196          gMC->Gtreve_root();
 197
 198          thed = fCurTheta*kRaddeg;
 199          phid = fCurPhi*kRaddeg;
 200          eta  = -TMath::Log(TMath::Tan(TMath::Max(
 201                      TMath::Min((Double_t)fCurTheta/2,TMath::Pi()/2-1.e-10),1.e-10)));
 202          fHistRadl->Fill(phid,thed,fTotRadl);
 203          fHistAbso->Fill(phid,thed,fTotAbso);
 204          fHistGcm2->Fill(phid,thed,fTotGcm2);
 205          fHistReta->Fill(phid,eta,fTotRadl);
 206          gAlice->FinishEvent();
 207       }
 208    }
 209    // store histograms in current Root file
 210    fHistRadl->Write();
 211    fHistAbso->Write();
 212    fHistGcm2->Write();
 213    fHistReta->Write();
 214 }
 215
 216 //___________________________________________
 217 void AliLego::StepManager()
 218 {
 219 // called from AliRun::Stepmanager from gustep.
 220 // Accumulate the 3 parameters step by step
 221
 222    Float_t t, tt;
 223    Float_t a,z,dens,radl,absl;
 224    Int_t i;
 225
 226    Float_t step  = gMC->TrackStep();
 227
 228    Float_t vect[3], dir[3];
 229    TLorentzVector pos, mom;
 230
 231    gMC->TrackPosition(pos);
 232    gMC->TrackMomentum(mom);
 233    gMC->CurrentMaterial(a,z,dens,radl,absl);
 234
 235    if (z < 1) return;
 236
 237 // --- See if we have to stop now
 238    if (TMath::Abs(pos[2]) > fZMax  ||
 239        pos[0]*pos[0] +pos[1]*pos[1] > fRadMax*fRadMax) {
 240        gMC->StopEvent();
 241    } else {
 242
 243 // --- See how long we have to go
 244       for(i=0;i<3;++i) {
 245         vect[i]=pos[i];
 246         dir[i]=mom[i];
 247       }
 248       t  = PropagateCylinder(vect,dir,fRadMax,fZMax);
 249       tt = TMath::Min(step,t);
 250
 251       fTotAbso += tt/absl;
 252       fTotRadl += tt/radl;
 253       fTotGcm2 += tt*dens;
 254    }
 255 }
 256