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="gif/alilego.gif">
  26 */
  27 //End_Html
  28 //
  29 //////////////////////////////////////////////////////////////
  30
  31 #include "TMath.h"
  32 #include "TGeant3.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    const Int_t mpart = 48;
  74    Float_t orig[3], pmom[3];
  75    Float_t t, cost, sint, cosp, sinp;
  76
  77 // --- Set to 0 radiation length, absorption length and g/cm2 ---
  78    fTotRadl = 0;
  79    fTotAbso = 0;
  80    fTotGcm2 = 0;
  81
  82    fCurTheta = (fThetaMin+(fThetaBin-0.5)*(fThetaMax-fThetaMin)/fNtheta)*kDegrad;
  83    fCurPhi   = (fPhiMin+(fPhiBin-0.5)*(fPhiMax-fPhiMin)/fNphi)*kDegrad;
  84    cost      = TMath::Cos(fCurTheta);
  85    sint      = TMath::Sin(fCurTheta);
  86    cosp      = TMath::Cos(fCurPhi);
  87    sinp      = TMath::Sin(fCurPhi);
  88
  89    pmom[0] = cosp*sint;
  90    pmom[1] = sinp*sint;
  91    pmom[2] = cost;
  92
  93 // --- Where to start
  94    orig[0] = orig[1] = orig[2] = 0;
  95    Float_t dalicz = 3000;
  96    if (fRadMin > 0) {
  97       t = PropagateCylinder(orig,pmom,fRadMin,dalicz);
  98       orig[0] = pmom[0]*t;
  99       orig[1] = pmom[1]*t;
 100       orig[2] = pmom[2]*t;
 101       if (TMath::Abs(orig[2]) > fZMax) return;
 102    }
 103
 104 // --- We do start here
 105    Float_t polar[3]={0.,0.,0.};
 106    Int_t ntr;
 107    gAlice->SetTrack(1, 0, mpart, pmom, orig, polar, 0, "LEGO ray", ntr);
 108 }
 109
 110 //___________________________________________
 111 void AliLego::Init(Int_t ntheta,Float_t themin,Float_t themax,
 112           Int_t nphi,Float_t phimin,Float_t phimax,Float_t rmin,Float_t rmax,
 113                   Float_t zmax)
 114 {
 115 // specify the angular limits and the size of the rectangular box
 116    fNtheta   = ntheta;
 117    fThetaMin = themin;
 118    fThetaMax = themax;
 119    fNphi     = nphi;
 120    fPhiMin   = phimin;
 121    fPhiMax   = phimax;
 122    fRadMin   = rmin;
 123    fRadMax   = rmax;
 124    fZMax     = zmax;
 125    Float_t etamin = -TMath::Log(TMath::Tan(TMath::Min((Double_t)fThetaMax*kDegrad/2,TMath::Pi()/2-1.e-10)));
 126    Float_t etamax = -TMath::Log(TMath::Tan(TMath::Max((Double_t)fThetaMin*kDegrad/2,              1.e-10)));
 127
 128    fHistRadl = new TH2F("hradl","Radiation length map",    nphi,phimin,phimax,ntheta,themin,themax);
 129    fHistAbso = new TH2F("habso","Interaction length map",  nphi,phimin,phimax,ntheta,themin,themax);
 130    fHistGcm2 = new TH2F("hgcm2","g/cm2 length map",        nphi,phimin,phimax,ntheta,themin,themax);
 131    fHistReta = new TH2F("hetar","Radiation length vs. eta",nphi,phimin,phimax,ntheta,etamin,etamax);
 132
 133 }
 134
 135 //___________________________________________
 136 Float_t AliLego::PropagateCylinder(Float_t *x, Float_t *v, Float_t r, Float_t z)
 137 {
 138 // Propagate to cylinder from inside
 139
 140    Double_t hnorm, sz, t, t1, t2, t3, sr;
 141    Double_t d[3];
 142    const Float_t kSmall  = 1e-8;
 143    const Float_t kSmall2 = kSmall*kSmall;
 144
 145 // ---> Find intesection with Z planes
 146    d[0]  = v[0];
 147    d[1]  = v[1];
 148    d[2]  = v[2];
 149    hnorm = TMath::Sqrt(1/(d[0]*d[0]+d[1]*d[1]+d[2]*d[2]));
 150    d[0] *= hnorm;
 151    d[1] *= hnorm;
 152    d[2] *= hnorm;
 153    if (d[2] > kSmall)       sz = (z-x[2])/d[2];
 154    else if (d[2] < -kSmall) sz = -(z+x[2])/d[2];
 155    else                     sz = 1.e10;  // ---> Direction parallel to X-Y, no intersection
 156
 157 // ---> Intersection with cylinders
 158 //      Intersection point (x,y,z)
 159 //      (x,y,z) is on track :    x=X(1)+t*D(1)
 160 //                               y=X(2)+t*D(2)
 161 //                               z=X(3)+t*D(3)
 162 //      (x,y,z) is on cylinder : x**2 + y**2 = R**2
 163 //
 164 //      (D(1)**2+D(2)**2)*t**2
 165 //      +2.*(X(1)*D(1)+X(2)*D(2))*t
 166 //      +X(1)**2+X(2)**2-R**2=0
 167 // ---> Solve second degree equation
 168    t1 = d[0]*d[0] + d[1]*d[1];
 169    if (t1 <= kSmall2) {
 170       t = sz;  // ---> Track parallel to the z-axis, take distance to planes
 171    } else {
 172       t2 = x[0]*d[0] + x[1]*d[1];
 173       t3 = x[0]*x[0] + x[1]*x[1];
 174       // ---> It should be positive, but there may be numerical problems
 175       sr = (t2 +TMath::Sqrt(TMath::Max(t2*t2-(t3-r*r)*t1,0.)))/t1;
 176       // ---> Find minimum distance between planes and cylinder
 177       t  = TMath::Min(sz,sr);
 178    }
 179    return t;
 180 }
 181
 182 //___________________________________________
 183 void AliLego::Run()
 184 {
 185    // loop on phi,theta bins
 186    AliMC* pMC=AliMC::GetMC();
 187    pMC->InitLego();
 188    Float_t thed, phid, eta;
 189    for (fPhiBin=1; fPhiBin<=fNphi; fPhiBin++) {
 190       printf("AliLego Generating rays in phi bin:%d\n",fPhiBin);
 191       for (fThetaBin=1; fThetaBin<=fNtheta; fThetaBin++) {
 192          pMC->Gtrigi();
 193          pMC->Gtrigc();
 194          GenerateKinematics();
 195          pMC->Gtreve();
 196
 197          thed = fCurTheta*kRaddeg;
 198          phid = fCurPhi*kRaddeg;
 199          eta  = -TMath::Log(TMath::Tan(TMath::Max(
 200                      TMath::Min((Double_t)fCurTheta/2,TMath::Pi()/2-1.e-10),1.e-10)));
 201          fHistRadl->Fill(phid,thed,fTotRadl);
 202          fHistAbso->Fill(phid,thed,fTotAbso);
 203          fHistGcm2->Fill(phid,thed,fTotGcm2);
 204          fHistReta->Fill(phid,eta,fTotRadl);
 205          gAlice->FinishEvent();
 206       }
 207    }
 208    // store histograms in current Root file
 209    fHistRadl->Write();
 210    fHistAbso->Write();
 211    fHistGcm2->Write();
 212    fHistReta->Write();
 213 }
 214
 215 //___________________________________________
 216 void AliLego::StepManager()
 217 {
 218 // called from AliRun::Stepmanager from gustep.
 219 // Accumulate the 3 parameters step by step
 220   AliMC* pMC = AliMC::GetMC();
 221
 222    Float_t t, tt;
 223    Float_t a,z,dens,radl,absl;
 224
 225    Float_t step  = pMC->TrackStep();
 226
 227    Float_t vect[3], pmom[4];
 228    pMC->TrackPosition(vect);
 229    pMC->TrackMomentum(pmom);
 230    pMC->CurrentMaterial(a,z,dens,radl,absl);
 231
 232    if (z < 1) return;
 233
 234 // --- See if we have to stop now
 235    if (TMath::Abs(vect[2]) > fZMax  ||
 236        vect[0]*vect[0] +vect[1]*vect[1] > fRadMax*fRadMax) {
 237        pMC->StopEvent();
 238    } else {
 239
 240 // --- See how long we have to go
 241       t  = PropagateCylinder(vect,pmom,fRadMax,fZMax);
 242       tt = TMath::Min(step,t);
 243
 244       fTotAbso += tt/absl;
 245       fTotRadl += tt/radl;
 246       fTotGcm2 += tt*dens;
 247    }
 248 }
 249