MUON/MUONmassPlot_NewIO.C

   1 #if !defined(__CINT__) || defined(__MAKECINT__)
   2 // ROOT includes
   3 #include "TBranch.h"
   4 #include "TClonesArray.h"
   5 #include "TLorentzVector.h"
   6 #include "TFile.h"
   7 #include "TH1.h"
   8 #include "TParticle.h"
   9 #include "TTree.h"
  10
  11 // STEER includes
  12 #include "AliRun.h"
  13 #include "AliRunLoader.h"
  14 #include "AliHeader.h"
  15 #include "AliLoader.h"
  16 #include "AliStack.h"
  17
  18 // MUON includes
  19 #include "AliMUON.h"
  20 #include "AliMUONData.h"
  21 #include "AliMUONHit.h"
  22 #include "AliMUONConstants.h"
  23 #include "AliMUONDigit.h"
  24 #include "AliMUONRawCluster.h"
  25 #include "AliMUONGlobalTrigger.h"
  26 #include "AliMUONLocalTrigger.h"
  27 #include "AliMUONTrack.h"
  28 #include "AliMUONTrackParam.h"
  29 #include "AliMUONTrackExtrap.h"
  30 #include "AliESDMuonTrack.h"
  31 #endif
  32 //
  33 // Macro MUONmassPlot.C for new I/O
  34 // Ch. Finck, Subatech, Jan. 2004
  35 //
  36
  37 // macro to make invariant mass plots
  38 // for combinations of 2 muons with opposite charges,
  39 // from root file "MUON.tracks.root" containing the result of track reconstruction.
  40 // Histograms are stored on the "MUONmassPlot.root" file.
  41 // introducing TLorentzVector for parameter calculations (Pt, P,rap,etc...)
  42 // using Invariant Mass for rapidity.
  43
  44 // Arguments:
  45 //   FirstEvent (default 0)
  46 //   LastEvent (default 0)
  47 //   ResType (default 553)
  48 //      553 for Upsilon, anything else for J/Psi
  49 //   Chi2Cut (default 100)
  50 //      to keep only tracks with chi2 per d.o.f. < Chi2Cut
  51 //   PtCutMin (default 1)
  52 //      to keep only tracks with transverse momentum > PtCutMin
  53 //   PtCutMax (default 10000)
  54 //      to keep only tracks with transverse momentum < PtCutMax
  55 //   massMin (default 9.17 for Upsilon)
  56 //      &  massMax (default 9.77 for Upsilon)
  57 //         to calculate the reconstruction efficiency for resonances with invariant mass
  58 //         massMin < mass < massMax.
  59
  60 // Add parameters and histograms for analysis
  61
  62 void MUONmassPlot(char* filename="galice.root", Int_t FirstEvent = 0, Int_t LastEvent = 0, Int_t ResType = 553,
  63                   Float_t Chi2Cut = 100., Float_t PtCutMin = 1., Float_t PtCutMax = 10000.,
  64                   Float_t massMin = 9.17,Float_t massMax = 9.77)
  65 {
  66   cout << "MUONmassPlot " << endl;
  67   cout << "FirstEvent " << FirstEvent << endl;
  68   cout << "LastEvent " << LastEvent << endl;
  69   cout << "ResType " << ResType << endl;
  70   cout << "Chi2Cut " << Chi2Cut << endl;
  71   cout << "PtCutMin " << PtCutMin << endl;
  72   cout << "PtCutMax " << PtCutMax << endl;
  73   cout << "massMin " << massMin << endl;
  74   cout << "massMax " << massMax << endl;
  75
  76
  77   //Reset ROOT and connect tree file
  78   gROOT->Reset();
  79
  80
  81   // File for histograms and histogram booking
  82   TFile *histoFile = new TFile("MUONmassPlot.root", "RECREATE");
  83   TH1F *hPtMuon = new TH1F("hPtMuon", "Muon Pt (GeV/c)", 100, 0., 20.);
  84   TH1F *hPMuon = new TH1F("hPMuon", "Muon P (GeV/c)", 100, 0., 200.);
  85   TH1F *hChi2PerDof = new TH1F("hChi2PerDof", "Muon track chi2/d.o.f.", 100, 0., 20.);
  86   TH1F *hInvMassAll = new TH1F("hInvMassAll", "Mu+Mu- invariant mass (GeV/c2)", 480, 0., 12.);
  87   TH1F *hInvMassRes;
  88
  89   if (ResType == 553) {
  90     hInvMassRes = new TH1F("hInvMassRes", "Mu+Mu- invariant mass (GeV/c2) around Upsilon", 60, 8., 11.);
  91   } else {
  92     hInvMassRes = new TH1F("hInvMassRes", "Mu+Mu- invariant mass (GeV/c2) around J/Psi", 80, 0., 5.);
  93   }
  94
  95   TH1F *hNumberOfTrack = new TH1F("hNumberOfTrack","nb of track /evt ",20,-0.5,19.5);
  96   TH1F *hRapMuon = new TH1F("hRapMuon"," Muon Rapidity",50,-4.5,-2);
  97   TH1F *hRapResonance = new TH1F("hRapResonance"," Resonance Rapidity",50,-4.5,-2);
  98   TH1F *hPtResonance = new TH1F("hPtResonance", "Resonance Pt (GeV/c)", 100, 0., 20.);
  99
 100
 101   // settings
 102   Int_t EventInMass = 0;
 103   Float_t muonMass = 0.105658389;
 104 //   Float_t UpsilonMass = 9.46037;
 105 //   Float_t JPsiMass = 3.097;
 106
 107   Double_t bendingSlope, nonBendingSlope, pYZ;
 108   Double_t fPxRec1, fPyRec1, fPzRec1, fZRec1, fE1;
 109   Double_t fPxRec2, fPyRec2, fPzRec2, fZRec2, fE2;
 110   Int_t fCharge, fCharge2;
 111
 112   Int_t ntrackhits, nevents;
 113   Double_t fitfmin;
 114
 115   TClonesArray * recTracksArray;
 116   TLorentzVector fV1, fV2, fVtot;
 117
 118   // Creating Run Loader and openning file containing Hits
 119   AliRunLoader * RunLoader = AliRunLoader::Open(filename,"MUONFolder","READ");
 120   if (RunLoader == 0x0) {
 121     printf(">>> Error : Error Opening %s file \n",filename);
 122     return;
 123   }
 124
 125   AliLoader * MUONLoader = RunLoader->GetLoader("MUONLoader");
 126   MUONLoader->LoadTracks("READ");
 127
 128   // Creating MUON data container
 129   AliMUONData muondata(MUONLoader,"MUON","MUON");
 130
 131   nevents = RunLoader->GetNumberOfEvents();
 132
 133   AliMUONTrack * rectrack;
 134   AliMUONTrackParam *trackParam;
 135
 136   // Loop over events
 137   for (Int_t ievent = FirstEvent; ievent <= TMath::Min(LastEvent, nevents - 1); ievent++) {
 138
 139     // get current event
 140     RunLoader->GetEvent(ievent);
 141
 142     muondata.SetTreeAddress("RT");
 143     muondata.GetRecTracks();
 144     recTracksArray = muondata.RecTracks();
 145
 146     Int_t nrectracks = (Int_t) recTracksArray->GetEntriesFast(); //
 147
 148     printf("\n Nb of events analysed: %d\r",ievent);
 149     //   cout << " number of tracks: " << nrectracks  <<endl;
 150
 151     // loop over all reconstructed tracks (also first track of combination)
 152      for (Int_t irectracks = 0; irectracks <  nrectracks;  irectracks++) {
 153
 154       rectrack = (AliMUONTrack*) recTracksArray->At(irectracks);
 155
 156       trackParam = (AliMUONTrackParam*)(rectrack->GetTrackParamAtHit()->First());
 157       AliMUONTrackExtrap::ExtrapToVertex(trackParam,0.,0.,0.);
 158       bendingSlope   = trackParam->GetBendingSlope();
 159       nonBendingSlope = trackParam->GetNonBendingSlope();
 160
 161       pYZ     = 1/TMath::Abs(trackParam->GetInverseBendingMomentum());
 162       fPzRec1  = - pYZ / TMath::Sqrt(1.0 + bendingSlope * bendingSlope); // spectro. (z<0)
 163       fPxRec1  = fPzRec1 * nonBendingSlope;
 164       fPyRec1  = fPzRec1 * bendingSlope;
 165       fZRec1   = trackParam->GetZ();
 166       fCharge = Int_t(TMath::Sign(1., trackParam->GetInverseBendingMomentum()));
 167
 168       fE1 = TMath::Sqrt(muonMass * muonMass + fPxRec1 * fPxRec1 + fPyRec1 * fPyRec1 + fPzRec1 * fPzRec1);
 169       fV1.SetPxPyPzE(fPxRec1, fPyRec1, fPzRec1, fE1);
 170
 171       ntrackhits = rectrack->GetNTrackHits();
 172       fitfmin = rectrack->GetFitFMin();
 173
 174       // transverse momentum
 175       Float_t pt1 = fV1.Pt();
 176
 177       // total momentum
 178       Float_t p1 = fV1.P();
 179
 180       // Rapidity
 181       Float_t rapMuon1 = fV1.Rapidity();
 182
 183       // chi2 per d.o.f.
 184       Float_t ch1 =  fitfmin / (2.0 * ntrackhits - 5);
 185 //      printf(" px %f py %f pz %f NHits %d  Norm.chi2 %f charge %d\n",
 186 //           fPxRec1, fPyRec1, fPzRec1, ntrackhits, ch1, fCharge);
 187
 188       // condition for good track (Chi2Cut and PtCut)
 189
 190       if ((ch1 < Chi2Cut) && (pt1 > PtCutMin) && (pt1 < PtCutMax)) {
 191
 192         // fill histos hPtMuon and hChi2PerDof
 193         hPtMuon->Fill(pt1);
 194         hPMuon->Fill(p1);
 195         hChi2PerDof->Fill(ch1);
 196         hRapMuon->Fill(rapMuon1);
 197
 198         // loop over second track of combination
 199         for (Int_t irectracks2 = irectracks + 1; irectracks2 < nrectracks; irectracks2++) {
 200
 201           rectrack = (AliMUONTrack*) recTracksArray->At(irectracks2);
 202
 203           trackParam = (AliMUONTrackParam*)(rectrack->GetTrackParamAtHit()->First());
 204           AliMUONTrackExtrap::ExtrapToVertex(trackParam,0.,0.,0.);
 205           bendingSlope    = trackParam->GetBendingSlope();
 206           nonBendingSlope = trackParam->GetNonBendingSlope();
 207
 208           pYZ      = 1/TMath::Abs(trackParam->GetInverseBendingMomentum());
 209           fPzRec2  = - pYZ / TMath::Sqrt(1.0 + bendingSlope * bendingSlope); // spectro. (z<0)
 210           fPxRec2  = fPzRec2 * nonBendingSlope;
 211           fPyRec2  = fPzRec2 * bendingSlope;
 212           fZRec2   = trackParam->GetZ();
 213           fCharge2 = Int_t(TMath::Sign(1., trackParam->GetInverseBendingMomentum()));
 214
 215           fE2 = TMath::Sqrt(muonMass * muonMass + fPxRec2 * fPxRec2 + fPyRec2 * fPyRec2 + fPzRec2 * fPzRec2);
 216           fV2.SetPxPyPzE(fPxRec2, fPyRec2, fPzRec2, fE2);
 217
 218           ntrackhits = rectrack->GetNTrackHits();
 219           fitfmin = rectrack->GetFitFMin();
 220
 221           // transverse momentum
 222           Float_t pt2 = fV2.Pt();
 223
 224           // chi2 per d.o.f.
 225           Float_t ch2 = fitfmin  / (2.0 * ntrackhits - 5);
 226
 227           // condition for good track (Chi2Cut and PtCut)
 228           if ((ch2 < Chi2Cut) && (pt2 > PtCutMin)  && (pt2 < PtCutMax)) {
 229
 230             // condition for opposite charges
 231             if ((fCharge * fCharge2) == -1) {
 232
 233               // invariant mass
 234               fVtot = fV1 + fV2;
 235               Float_t invMass = fVtot.M();
 236
 237               // fill histos hInvMassAll and hInvMassRes
 238               hInvMassAll->Fill(invMass);
 239               hInvMassRes->Fill(invMass);
 240
 241               if (invMass > massMin && invMass < massMax) {
 242                 EventInMass++;
 243                 hRapResonance->Fill(fVtot.Rapidity());
 244                 hPtResonance->Fill(fVtot.Pt());
 245               }
 246
 247             } // if (fCharge * fCharge2) == -1)
 248           } // if ((ch2 < Chi2Cut) && (pt2 > PtCutMin) && (pt2 < PtCutMax))
 249         } //  for (Int_t irectracks2 = irectracks + 1; irectracks2 < irectracks; irectracks2++)
 250       } // if (ch1 < Chi2Cut) && (pt1 > PtCutMin)&& (pt1 < PtCutMax) )
 251     } // for (Int_t irectracks = 0; irectracks < nrectracks; irectracks++)
 252
 253     hNumberOfTrack->Fill(nrectracks);
 254   } // for (Int_t ievent = FirstEvent;
 255
 256   histoFile->Write();
 257   histoFile->Close();
 258
 259   cout << "MUONmassPlot " << endl;
 260   cout << "FirstEvent " << FirstEvent << endl;
 261   cout << "LastEvent " << LastEvent << endl;
 262   cout << "ResType " << ResType << endl;
 263   cout << "Chi2Cut " << Chi2Cut << endl;
 264   cout << "PtCutMin " << PtCutMin << endl;
 265   cout << "PtCutMax " << PtCutMax << endl;
 266   cout << "massMin " << massMin << endl;
 267   cout << "massMax " << massMax << endl;
 268   cout << "EventInMass " << EventInMass << endl;
 269 }
 270