MUON/MUONTestAbso.C

   1 void MUONTestAbso (Int_t evNumber1=0,Int_t evNumber2=0)
   2 {
   3   // Useful to test the absorber correction in the reconstruction procedure
   4   // (mean energy loss + Branson correction)
   5   // The AliMUONTrackParam class from the reconstruction is directly checked
   6   // in this macro using the GEANT particle momentum downstream of the absorber.
   7   // Histograms are saved in file MUONTestAbso.root.
   8   // Use DrawTestAbso.C to display control histograms.
   9
  10
  11   // Dynamically link some shared libs
  12   if (gClassTable->GetID("AliRun") < 0) {
  13     gROOT->LoadMacro("loadlibs.C");
  14     loadlibs();
  15   }
  16
  17   // Connect the Root Galice file containing Geometry, Kine and Hits
  18
  19   TFile *file = (TFile*)gROOT->GetListOfFiles()->FindObject("galice.root");
  20
  21
  22   if (!file) {
  23     printf("\n Creating galice.root \n");
  24     file = new TFile("galice.root");
  25   } else {
  26     printf("\n galice.root found in file list");
  27   }
  28
  29 // Get AliRun object from file or create it if not on file
  30   if (!gAlice) {
  31     gAlice = (AliRun*)(file->Get("gAlice"));
  32     if (gAlice) printf("AliRun object found on file\n");
  33     if (!gAlice) {
  34       printf("\n Create new gAlice object");
  35       gAlice = new AliRun("gAlice","Alice test program");
  36     }
  37   }
  38
  39
  40   //  Create some histograms
  41
  42   TH1F *hInvMassRes = new TH1F("hInvMassRes", "Mu+Mu- invariant mass (GeV/c2) around Upsilon", 120, 8., 11.);
  43   TH1F *hDeltaP1 = new TH1F("hDeltaP1", " delta P for muons r < rLimit", 100, -10., 10.);
  44   TH1F *hDeltaAngle1 = new TH1F("hDeltaAngle1", " delta angle for muons r < rLimit", 100, 0., 0.5);
  45   //
  46   TH1F *hDeltaP2 = new TH1F("hDeltaP2", " delta P for muons r > rLimit", 100, -10., 10.);
  47   TH1F *hDeltaAngle2 = new TH1F("hDeltaAngle2", " delta angle for muons r > rLimit", 100, 0., 0.5);
  48
  49   // Define constant parameters
  50   Float_t massMuon = 0.105658389;  // muon mass
  51   Float_t massUpsilon = 9.46037;   // Upsilon mass
  52   Double_t PI = 3.141592654;
  53   Double_t zEndAbso = 503;         // z position of the end of the absorber
  54   Double_t rLimit = zEndAbso * TMath::Tan(3*PI/180); // 3 deg. limit (different absorber composition)
  55   Double_t printLevel = 0;
  56
  57   if (printLevel >= 1)
  58   cout <<" **** z end Absorber : "<< zEndAbso <<" rLimit : "<< rLimit << endl;
  59
  60   Double_t pX[2],pY[2],pZ[2],pTot[2],theta[2],radius[2]; // extrapolated parameters of particles from chamber 1 to vertex
  61
  62   for (Int_t nev=0; nev<= evNumber2; nev++) {  // loop over events
  63     if (printLevel >= 0) {
  64       cout << "  "<<endl;
  65       cout << " **** Event # " << nev <<endl;
  66     }
  67
  68     Int_t nparticles = gAlice->GetEvent(nev);
  69     if (printLevel >= 1) {
  70       cout << " Total number of nparticles  " << nparticles <<endl;
  71     }
  72     if (nev < evNumber1) continue;
  73     if (nparticles <= 0) return;
  74
  75     Double_t PxGen[2],PyGen[2],PzGen[2],PTotGen[2],ThetaGen[2]; // parameters of generated  particles at the vertex
  76
  77     for (int iPart = 0; iPart < nparticles  ; iPart++) {  // loop over particles
  78
  79       TParticle *particle;
  80       particle = gAlice->Particle(iPart);
  81
  82       if ( particle->GetPdgCode() == kMuonPlus ) {
  83         PxGen[0] = particle->Px();
  84         PyGen[0] = particle->Py();
  85         PzGen[0] = particle->Pz();
  86         PTotGen[0] = TMath::Sqrt(PxGen[0]*PxGen[0]+PyGen[0]*PyGen[0]+PzGen[0]*PzGen[0]);
  87         ThetaGen[0] = TMath::ATan2(TMath::Sqrt(PxGen[0]*PxGen[0]+PyGen[0]*PyGen[0]),PzGen[0])*180/PI;
  88         if (printLevel >= 1) {
  89           cout << " Particle id : "<<particle->GetPdgCode()<<" px,py,pz,theta : "<<PxGen[0]<<" "<<PyGen[0]<<" "<<PzGen[0]<<" "<<ThetaGen[0]<<endl;
  90         }
  91       }
  92       if ( particle->GetPdgCode() == kMuonMinus ) {
  93         PxGen[1] = particle->Px();
  94         PyGen[1] = particle->Py();
  95         PzGen[1] = particle->Pz();
  96         PTotGen[1] = TMath::Sqrt(PxGen[1]*PxGen[1]+PyGen[1]*PyGen[1]+PzGen[1]*PzGen[1]);
  97         ThetaGen[1] = TMath::ATan2(TMath::Sqrt(PxGen[1]*PxGen[1]+PyGen[1]*PyGen[1]),PzGen[1])*180/PI;
  98         if (printLevel >= 1)
  99           cout << " Particle #: "<<particle->GetPdgCode()<<" px,py,pz,theta : "<<PxGen[1]<<" "<<PyGen[1]<<" "<<PzGen[1]<<" "<<ThetaGen[1]<<endl;
 100       }
 101     }
 102
 103
 104     AliMUON *MUON  = (AliMUON*) gAlice->GetModule("MUON");
 105
 106
 107     TTree *TH = gAlice->TreeH();
 108     Int_t ntracks =(Int_t) TH->GetEntries();
 109
 110     for (Int_t i = 0; i < 2; i++){
 111       pX[i] = pY[i] = pZ[i] = pTot[i] = theta[i] = radius[i] = 0;
 112     }
 113
 114     for (Int_t track=0; track<ntracks;track++) { // loop over tracks
 115       gAlice->ResetHits();
 116       Int_t nbytes += TH->GetEvent(track);
 117       if (MUON)  {
 118         for(AliMUONHit* mHit=(AliMUONHit*)MUON->FirstHit(-1);
 119             mHit;
 120             mHit=(AliMUONHit*)MUON->NextHit())  // loop over GEANT muon hits
 121           {
 122             Int_t   nch   = mHit->Chamber();  // chamber number
 123
 124             if (nch != 1) continue;
 125
 126             Int_t ipart = mHit->Particle();
 127             Double_t x     = mHit->X();        // x-pos of hit in chamber #1
 128             Double_t y     = mHit->Y();        // y-pos of hit in chamber #1
 129             Double_t z     = mHit->Z();        // z-pos of hit in chamber #1
 130
 131
 132             if (ipart == kMuonPlus || ipart == kMuonMinus) {
 133               Double_t px0=mHit->Px();
 134               Double_t py0=mHit->Py();
 135               Double_t pz0=mHit->Pz();
 136               Double_t nonBendingSlope=px0/pz0;
 137               Double_t bendingSlope=py0/pz0;
 138
 139               // create an AliMUONTrackParam object in chamber #1
 140               AliMUONTrackParam *trackParam = new AliMUONTrackParam();
 141               Double_t bendingMometum = TMath::Sqrt(pz0*pz0+px0*px0);
 142               Double_t inverseBendingMomentum = 1/bendingMometum;
 143
 144               trackParam->SetBendingCoor(y);
 145               trackParam->SetNonBendingCoor(x);
 146               trackParam->SetInverseBendingMomentum(inverseBendingMomentum);
 147               trackParam->SetBendingSlope(bendingSlope);
 148               trackParam->SetNonBendingSlope(nonBendingSlope);
 149               trackParam->SetZ(z);
 150
 151               if (printLevel >= 1) {
 152                 cout <<" **** before extrap " <<endl;
 153                 cout << " px, py, pz = " <<px0<<" "<<py0<<" "<<pz0<<endl;
 154                 trackParam->Dump();
 155               }
 156               trackParam->ExtrapToVertex(); // extrapolate track parameters through the absorber
 157               if (printLevel >= 1) {
 158                 cout <<" **** after extrap " <<endl;
 159                 trackParam->Dump();
 160               }
 161
 162               Int_t iMuon;
 163               if (ipart == kMuonPlus) iMuon = 0;
 164               else iMuon = 1;
 165
 166               // calculate track parameters extrapolated to vertex.
 167               bendingSlope = trackParam->GetBendingSlope();
 168               nonBendingSlope = trackParam->GetNonBendingSlope();
 169               Double_t pYZ = 1/TMath::Abs(trackParam->GetInverseBendingMomentum());
 170               pZ[iMuon] = pYZ/TMath::Sqrt(1+bendingSlope*bendingSlope);
 171               pX[iMuon] = pZ[iMuon] * nonBendingSlope;
 172               pY[iMuon] = pZ[iMuon] * bendingSlope;
 173               pTot[iMuon] = TMath::Sqrt(pYZ*pYZ+pX[iMuon]*pX[iMuon]);
 174               theta[iMuon] = TMath::ATan2(TMath::Sqrt(pX[iMuon]*pX[iMuon]+pY[iMuon]*pY[iMuon]),pZ[iMuon])*180/PI;
 175               radius[iMuon] = TMath::Sqrt(x*x+y*y);
 176
 177               if (printLevel >= 1)
 178                 cout << " px, py, pz, theta, radius = " <<pX[iMuon]<<" "<<pY[iMuon]<<" "<<pZ[iMuon]<<" "<<theta[iMuon]<<" "<<radius[iMuon]<<endl;
 179               delete trackParam;
 180
 181             } // if MuonPlus or MuonMinus
 182           } // loop over MUON hits
 183       } // if MUON
 184     } // loop over tracks
 185
 186     if (pX[0] != 0 && pX[1] != 0) { // if track 1 & 2 exist
 187       Float_t massResonance = MuPlusMuMinusMass(pX[0],pY[0],pZ[0],pX[1],pY[1],pZ[1],massMuon);
 188       hInvMassRes->Fill(massResonance);
 189       for (Int_t i = 0; i<2; i++){
 190         Double_t cosAngle = pX[i]*PxGen[i]+pY[i]*PyGen[i]+pZ[i]*PzGen[i];
 191         cosAngle = cosAngle/(pTot[i]*PTotGen[i]);
 192         Double_t deltaAngle = TMath::ACos(cosAngle)*180/PI;
 193         if (radius[i] < rLimit) {
 194           hDeltaP1->Fill(pTot[i]-PTotGen[i]);
 195           hDeltaAngle1->Fill(deltaAngle);
 196         }else{
 197           hDeltaP2->Fill(pTot[i]-PTotGen[i]);
 198           hDeltaAngle2->Fill(deltaAngle);
 199         }
 200       }
 201     }
 202   } // loop over event
 203
 204   // save histograms in MUONTestAbso.root
 205   TFile *histoFile = new TFile("MUONTestAbso.root", "RECREATE");
 206   hInvMassRes->Write();
 207   hDeltaP1->Write();
 208   hDeltaAngle1->Write();
 209   hDeltaP2->Write();
 210   hDeltaAngle2->Write();
 211   histoFile->Close();
 212
 213 }
 214
 215
 216 Float_t MuPlusMuMinusMass(Float_t Px1, Float_t Py1, Float_t Pz1, Float_t Px2, Float_t Py2, Float_t Pz2, Float_t MuonMass)
 217 {
 218   // return invariant mass for particle 1 & 2 whose masses are equal to MuonMass
 219
 220   Float_t e1 = TMath::Sqrt(MuonMass * MuonMass + Px1 * Px1 + Py1 * Py1 + Pz1 * Pz1);
 221   Float_t e2 = TMath::Sqrt(MuonMass * MuonMass + Px2 * Px2 + Py2 * Py2 + Pz2 * Pz2);
 222   return (TMath::Sqrt(2.0 * (MuonMass * MuonMass + e1 * e2 - Px1 * Px2 - Py1 * Py2 - Pz1 * Pz2)));
 223 }