MUON/MUONmassPlot.C

   1 // macro to make invariant mass plots
   2 // for combinations of 2 muons with opposite charges,
   3 // from root file "MUONtrackReco.root" containing the result of track reconstruction,
   4 // generated by the macro "MUONrecoNtuple.C".
   5 // Histograms are stored on the "MUONmassPlot.root" file.
   6 // A model for macros using the Ntuple in the file "MUONtrackReco.root"
   7 // may be found in "MUONtrackRecoModel.C":
   8 // it has been obtained by reading with Root a file "MUONtrackReco.root",
   9 // and executing the command:
  10 // MUONtrackReco->MakeCode("MUONtrackRecoModel.C")
  11
  12 // Arguments:
  13 //   FirstEvent (default 0)
  14 //   LastEvent (default 0)
  15 //   ResType (default 553)
  16 //      553 for Upsilon, anything else for J/Psi
  17 //   Chi2Cut (default 100)
  18 //      to keep only tracks with chi2 per d.o.f. < Chi2Cut
  19 //   PtCutMin (default 1)
  20 //      to keep only tracks with transverse momentum > PtCutMin
  21 //   PtCutMax (default 10000)
  22 //      to keep only tracks with transverse momentum < PtCutMax
  23 //   massMin (default 9.17 for Upsilon)
  24 //      &  massMax (default 9.77 for Upsilon)
  25 //         to calculate the reconstruction efficiency for resonances with invariant mass
  26 //         massMin < mass < massMax.
  27
  28 // compile MUONrecoNtuple.C and load shared library in this macro
  29 // Add parameters and histograms for analysis
  30
  31 void MUONmassPlot(Int_t FirstEvent = 0, Int_t LastEvent = 0, Int_t ResType = 553, Float_t Chi2Cut = 100., Float_t PtCutMin = 1., Float_t PtCutMax = 10000.,Float_t massMin = 9.17,Float_t massMax = 9.77)
  32 {
  33   cout << "MUONmassPlot " << endl;
  34   cout << "FirstEvent " << FirstEvent << endl;
  35   cout << "LastEvent " << LastEvent << endl;
  36   cout << "ResType " << ResType << endl;
  37 //   cout << "Nsig " << Nsig << endl;
  38   cout << "Chi2Cut " << Chi2Cut << endl;
  39   cout << "PtCutMin " << PtCutMin << endl;
  40   cout << "PtCutMax " << PtCutMax << endl;
  41   cout << "massMin " << massMin << endl;
  42   cout << "massMax " << massMax << endl;
  43
  44
  45   //////////////////////////////////////////////////////////
  46   //   This file has been automatically generated
  47   //     (Thu Sep 21 14:53:11 2000 by ROOT version2.25/02)
  48   //   from TTree MUONtrackReco/MUONtrackReco
  49   //   found on file: MUONtrackReco.root
  50   //////////////////////////////////////////////////////////
  51
  52
  53   //Reset ROOT and connect tree file
  54   gROOT->Reset();
  55 // Dynamically link some shared libs
  56
  57   if (gClassTable->GetID("AliRun") < 0) {
  58     gSystem->SetIncludePath("-I$ALICE_ROOT/MUON -I$ALICE_ROOT/STEER");
  59     gROOT->LoadMacro("loadlibs.C");
  60     loadlibs();
  61 // compile MUONrecoNtuple and load shared library
  62     gSystem->CompileMacro("$(ALICE_ROOT)/macros/MUONrecoNtuple.C","kf");
  63   }
  64
  65   TFile *f = (TFile*)gROOT->GetListOfFiles()->FindObject("MUONtrackReco.root");
  66   if (!f) {
  67     f = new TFile("MUONtrackReco.root");
  68   }
  69   TTree *MUONtrackReco = (TTree*)gDirectory->Get("MUONtrackReco");
  70   MUONtrackReco->SetMakeClass(1);
  71
  72 //Declaration of leaves types
  73   Int_t           fEvent;
  74   UInt_t          fUniqueID;
  75   UInt_t          fBits;
  76   Int_t           Tracks_;
  77   Int_t           Tracks_fCharge[500];
  78   Float_t         Tracks_fPxRec[500];
  79   Float_t         Tracks_fPyRec[500];
  80   Float_t         Tracks_fPzRec[500];
  81   Float_t         Tracks_fZRec[500];
  82   Float_t         Tracks_fZRec1[500];
  83   Int_t           Tracks_fNHits[500];
  84   Float_t         Tracks_fChi2[500];
  85   Float_t         Tracks_fPxGen[500];
  86   Float_t         Tracks_fPyGen[500];
  87   Float_t         Tracks_fPzGen[500];
  88   UInt_t          Tracks_fUniqueID[500];
  89   UInt_t          Tracks_fBits[500];
  90
  91   //Set branch addresses
  92   //MUONtrackReco->SetBranchAddress("Header",&Header);
  93   MUONtrackReco->SetBranchAddress("fEvent",&fEvent);
  94   MUONtrackReco->SetBranchAddress("fUniqueID",&fUniqueID);
  95   MUONtrackReco->SetBranchAddress("fBits",&fBits);
  96 //   MUONtrackReco->SetBranchAddress("Tracks_",&Tracks_);
  97   MUONtrackReco->SetBranchAddress("Tracks",&Tracks_);
  98   MUONtrackReco->SetBranchAddress("Tracks.fCharge",Tracks_fCharge);
  99   MUONtrackReco->SetBranchAddress("Tracks.fPxRec",Tracks_fPxRec);
 100   MUONtrackReco->SetBranchAddress("Tracks.fPyRec",Tracks_fPyRec);
 101   MUONtrackReco->SetBranchAddress("Tracks.fPzRec",Tracks_fPzRec);
 102   MUONtrackReco->SetBranchAddress("Tracks.fZRec",Tracks_fZRec);
 103   MUONtrackReco->SetBranchAddress("Tracks.fZRec1",Tracks_fZRec1);
 104   MUONtrackReco->SetBranchAddress("Tracks.fNHits",Tracks_fNHits);
 105   MUONtrackReco->SetBranchAddress("Tracks.fChi2",Tracks_fChi2);
 106   MUONtrackReco->SetBranchAddress("Tracks.fPxGen",Tracks_fPxGen);
 107   MUONtrackReco->SetBranchAddress("Tracks.fPyGen",Tracks_fPyGen);
 108   MUONtrackReco->SetBranchAddress("Tracks.fPzGen",Tracks_fPzGen);
 109   MUONtrackReco->SetBranchAddress("Tracks.fUniqueID",Tracks_fUniqueID);
 110   MUONtrackReco->SetBranchAddress("Tracks.fBits",Tracks_fBits);
 111
 112 //     This is the loop skeleton
 113 //       To read only selected branches, Insert statements like:
 114 // MUONtrackReco->SetBranchStatus("*",0);  // disable all branches
 115 // TTreePlayer->SetBranchStatus("branchname",1);  // activate branchname
 116
 117   Int_t nentries = MUONtrackReco->GetEntries();
 118
 119   Int_t nbytes = 0;
 120   //   for (Int_t i=0; i<nentries;i++) {
 121   //      nbytes += MUONtrackReco->GetEntry(i);
 122   //   }
 123
 124   /////////////////////////////////////////////////////////////////
 125   // Here comes the specialized part for MUONmassPlot
 126   /////////////////////////////////////////////////////////////////
 127
 128   // File for histograms and histogram booking
 129   TFile *histoFile = new TFile("MUONmassPlot.root", "RECREATE");
 130   TH1F *hPtMuon = new TH1F("hPtMuon", "Muon Pt (GeV/c)", 100, 0., 20.);
 131   TH1F *hPMuon = new TH1F("hPMuon", "Muon P (GeV/c)", 100, 0., 200.);
 132   TH1F *hChi2PerDof = new TH1F("hChi2PerDof", "Muon track chi2/d.o.f.", 100, 0., 20.);
 133   TH1F *hInvMassAll = new TH1F("hInvMassAll", "Mu+Mu- invariant mass (GeV/c2)", 480, 0., 12.);
 134   if (ResType == 553) {
 135     TH1F *hInvMassRes = new TH1F("hInvMassRes", "Mu+Mu- invariant mass (GeV/c2) around Upsilon", 60, 8., 11.);
 136   } else {
 137     TH1F *hInvMassRes = new TH1F("hInvMassRes", "Mu+Mu- invariant mass (GeV/c2) around J/Psi", 80, 0., 5.);
 138   }
 139
 140   TH1F *hNumberOfTrack = new TH1F("hNumberOfTrack","nb of track /evt ",20,-0.5,19.5);
 141   TH1F *hRapMuon = new TH1F("hRapMuon"," Muon Rapidity",50,2.,4.5);
 142   TH1F *hRapResonance = new TH1F("hRapResonance"," Resonance Rapidity",50,2.,4.5);
 143   TH1F *hPtResonance = new TH1F("hPtResonance", "Resonance Pt (GeV/c)", 100, 0., 20.);
 144   Int_t EventInMass = 0;
 145   Float_t muonMass = 0.105658389;
 146   Float_t UpsilonMass = 9.46037;
 147   Float_t JPsiMass = 3.097;
 148   // Loop over events
 149   for (Int_t event = FirstEvent; event <= TMath::Min(LastEvent, nentries - 1); event++) {
 150     // get current event
 151     Int_t NumberOfTrack = 0;
 152     cout <<" event: " << event <<endl;
 153     cout << " number of tracks: " << Tracks_ <<endl;
 154     nbytes += MUONtrackReco->GetEntry(event);
 155     // loop over all reconstructed tracks (also first track of combination)
 156     for (Int_t t1 = 0; t1 < Tracks_; t1++) {
 157       // transverse momentum
 158       Float_t pt1 = TMath::Sqrt(Tracks_fPxRec[t1] * Tracks_fPxRec[t1] + Tracks_fPyRec[t1] * Tracks_fPyRec[t1]);
 159
 160       // total momentum
 161       Float_t p1 = TMath::Sqrt(Tracks_fPxRec[t1] * Tracks_fPxRec[t1] + Tracks_fPyRec[t1] * Tracks_fPyRec[t1] + Tracks_fPzRec[t1] * Tracks_fPzRec[t1] );
 162       // Rapidity
 163       Float_t rapMuon1 = rapParticle(Tracks_fPxRec[t1],Tracks_fPyRec[t1],Tracks_fPzRec[t1],muonMass);
 164
 165       // chi2 per d.o.f.
 166       Float_t ch1 = Tracks_fChi2[t1] / (2.0 * Tracks_fNHits[t1] - 5);
 167       printf(" px %f py %f pz %f NHits %d  Norm.chi2 %f \n",Tracks_fPxRec[t1],Tracks_fPyRec[t1],Tracks_fPzRec[t1],Tracks_fNHits[t1],ch1);
 168       // condition for good track (Chi2Cut and PtCut)
 169       if ((ch1 < Chi2Cut) && (pt1 > PtCutMin) && (pt1 < PtCutMax)) {
 170         // fill histos hPtMuon and hChi2PerDof
 171         NumberOfTrack++;
 172         hPtMuon->Fill(pt1);
 173         hPMuon->Fill(p1);
 174         hChi2PerDof->Fill(ch1);
 175         hRapMuon->Fill(rapMuon1);
 176         // loop over second track of combination
 177         for (Int_t t2 = t1 + 1; t2 < Tracks_; t2++) {
 178           // transverse momentum
 179           Float_t pt2 = TMath::Sqrt(Tracks_fPxRec[t2] * Tracks_fPxRec[t2] + Tracks_fPyRec[t2] * Tracks_fPyRec[t2]);
 180           // chi2 per d.o.f.
 181           Float_t ch2 = Tracks_fChi2[t2] / (2.0 * Tracks_fNHits[t2] - 5);
 182           // condition for good track (Chi2Cut and PtCut)
 183           if ((ch2 < Chi2Cut) && (pt2 > PtCutMin)  && (pt2 < PtCutMax)) {
 184             // condition for opposite charges
 185             if ((Tracks_fCharge[t1] * Tracks_fCharge[t2]) == -1) {
 186               // invariant mass
 187               Float_t invMass = MuPlusMuMinusMass(Tracks_fPxRec[t1], Tracks_fPyRec[t1], Tracks_fPzRec[t1], Tracks_fPxRec[t2], Tracks_fPyRec[t2], Tracks_fPzRec[t2],muonMass);
 188               // fill histos hInvMassAll and hInvMassRes
 189               hInvMassAll->Fill(invMass);
 190               hInvMassRes->Fill(invMass);
 191               Float_t pxRes = Tracks_fPxRec[t1]+Tracks_fPxRec[t2];
 192               Float_t pyRes = Tracks_fPyRec[t1]+Tracks_fPyRec[t2];
 193               Float_t pzRes = Tracks_fPzRec[t1]+Tracks_fPzRec[t2];
 194               if (invMass > massMin && invMass < massMax) {
 195                 EventInMass++;
 196                 Float_t rapRes  = 0;
 197                 if (ResType == 553) {
 198                   rapRes  =  rapParticle(pxRes,pyRes,pzRes,UpsilonMass);
 199                 } else {
 200                   rapRes  =  rapParticle(pxRes,pyRes,pzRes,JPsiMass);
 201                 }
 202                 hRapResonance->Fill(rapRes);
 203                 hPtResonance->Fill(TMath::Sqrt(pxRes*pxRes+pyRes*pyRes));
 204               }
 205
 206             } //if ((Tracks_fCharge[t1] * Tracks_fCharge[t2]) == -1)
 207           } // if ((Tracks_fChi2[t2] < Chi2Cut) && (pt2 > PtCutMin) && (pt2 < PtCutMax))
 208         } // for (Int_t t2 = t1 + 1; t2 < Tracks_; t2++)
 209       } // if ((Tracks_fChi2[t1] < Chi2Cut) && (pt1 > PtCutMin)&& (pt1 < PtCutMax) )
 210     } // for (Int_t t1 = 0; t1 < Tracks_; t1++)
 211     hNumberOfTrack->Fill(NumberOfTrack);
 212   } // for (Int_t event = FirstEvent;
 213
 214   histoFile->Write();
 215   histoFile->Close();
 216
 217   cout << "MUONmassPlot " << endl;
 218   cout << "FirstEvent " << FirstEvent << endl;
 219   cout << "LastEvent " << LastEvent << endl;
 220   cout << "ResType " << ResType << endl;
 221 //   cout << "Nsig " << Nsig << endl;
 222   cout << "Chi2Cut " << Chi2Cut << endl;
 223   cout << "PtCutMin " << PtCutMin << endl;
 224   cout << "PtCutMax " << PtCutMax << endl;
 225   cout << "massMin " << massMin << endl;
 226   cout << "massMax " << massMax << endl;
 227   cout << "EventInMass " << EventInMass << endl;
 228 }
 229
 230
 231 Float_t MuPlusMuMinusMass(Float_t Px1, Float_t Py1, Float_t Pz1, Float_t Px2, Float_t Py2, Float_t Pz2, Float_t muonMass)
 232 {
 233   // return invariant mass for particle 1 & 2 whose masses are equal to muonMass
 234
 235   Float_t e1 = TMath::Sqrt(muonMass * muonMass + Px1 * Px1 + Py1 * Py1 + Pz1 * Pz1);
 236   Float_t e2 = TMath::Sqrt(muonMass * muonMass + Px2 * Px2 + Py2 * Py2 + Pz2 * Pz2);
 237   return (TMath::Sqrt(2.0 * (muonMass * muonMass + e1 * e2 - Px1 * Px2 - Py1 * Py2 - Pz1 * Pz2)));
 238 }
 239
 240 Float_t rapParticle(Float_t Px, Float_t Py, Float_t Pz, Float_t Mass)
 241 {
 242   // return rapidity for particle
 243
 244   // Particle energy
 245   Float_t Energy = TMath::Sqrt(Px*Px+Py*Py+Pz*Pz+Mass*Mass);
 246   // Rapidity
 247   Float_t Rapidity = 0.5*TMath::Log((Energy+Pz)/(Energy-Pz));
 248   return Rapidity;
 249 }