MUON/MUONmassPlot_ESD.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 "TH2.h"
   9 #include "TParticle.h"
  10 #include "TTree.h"
  11 #include <Riostream.h>
  12
  13 // STEER includes
  14 #include "AliRun.h"
  15 #include "AliRunLoader.h"
  16 #include "AliHeader.h"
  17 #include "AliLoader.h"
  18 #include "AliStack.h"
  19 #include "AliESD.h"
  20
  21 // MUON includes
  22 #include "AliESDMuonTrack.h"
  23 #endif
  24 //
  25 // Macro MUONmassPlot.C for ESD
  26 // Ch. Finck, Subatech, April. 2004
  27 //
  28
  29 // macro to make invariant mass plots
  30 // for combinations of 2 muons with opposite charges,
  31 // from root file "MUON.tracks.root" containing the result of track reconstruction.
  32 // Histograms are stored on the "MUONmassPlot.root" file.
  33 // introducing TLorentzVector for parameter calculations (Pt, P,rap,etc...)
  34 // using Invariant Mass for rapidity.
  35
  36 // Arguments:
  37 //   FirstEvent (default 0)
  38 //   LastEvent (default 0)
  39 //   ResType (default 553)
  40 //      553 for Upsilon, anything else for J/Psi
  41 //   Chi2Cut (default 100)
  42 //      to keep only tracks with chi2 per d.o.f. < Chi2Cut
  43 //   PtCutMin (default 1)
  44 //      to keep only tracks with transverse momentum > PtCutMin
  45 //   PtCutMax (default 10000)
  46 //      to keep only tracks with transverse momentum < PtCutMax
  47 //   massMin (default 9.17 for Upsilon)
  48 //      &  massMax (default 9.77 for Upsilon)
  49 //         to calculate the reconstruction efficiency for resonances with invariant mass
  50 //         massMin < mass < massMax.
  51
  52 // Add parameters and histograms for analysis
  53
  54 Bool_t MUONmassPlot(char* filename = "galice.root", Int_t FirstEvent = 0, Int_t LastEvent = 10000,
  55                   char* esdFileName = "AliESDs.root", Int_t ResType = 553,
  56                   Float_t Chi2Cut = 100., Float_t PtCutMin = 1., Float_t PtCutMax = 10000.,
  57                   Float_t massMin = 9.17,Float_t massMax = 9.77)
  58 {
  59   cout << "MUONmassPlot " << endl;
  60   cout << "FirstEvent " << FirstEvent << endl;
  61   cout << "LastEvent " << LastEvent << endl;
  62   cout << "ResType " << ResType << endl;
  63   cout << "Chi2Cut " << Chi2Cut << endl;
  64   cout << "PtCutMin " << PtCutMin << endl;
  65   cout << "PtCutMax " << PtCutMax << endl;
  66   cout << "massMin " << massMin << endl;
  67   cout << "massMax " << massMax << endl;
  68
  69
  70   //Reset ROOT and connect tree file
  71   gROOT->Reset();
  72
  73
  74   // File for histograms and histogram booking
  75   TFile *histoFile = new TFile("MUONmassPlot.root", "RECREATE");
  76   TH1F *hPtMuon = new TH1F("hPtMuon", "Muon Pt (GeV/c)", 100, 0., 20.);
  77   TH1F *hPtMuonPlus = new TH1F("hPtMuonPlus", "Muon+ Pt (GeV/c)", 100, 0., 20.);
  78   TH1F *hPtMuonMinus = new TH1F("hPtMuonMinus", "Muon- Pt (GeV/c)", 100, 0., 20.);
  79   TH1F *hPMuon = new TH1F("hPMuon", "Muon P (GeV/c)", 100, 0., 200.);
  80   TH1F *hChi2PerDof = new TH1F("hChi2PerDof", "Muon track chi2/d.o.f.", 100, 0., 20.);
  81   TH1F *hInvMassAll = new TH1F("hInvMassAll", "Mu+Mu- invariant mass (GeV/c2)", 480, 0., 12.);
  82   TH1F *hInvMassBg = new TH1F("hInvMassBg", "Mu+Mu- invariant mass BG(GeV/c2)", 480, 0., 12.);
  83 TH2F *hInvMassAll_vs_Pt = new TH2F("hInvMassAll_vs_Pt","hInvMassAll_vs_Pt",480,0.,12.,80,0.,20.);
  84 TH2F *hInvMassBgk_vs_Pt = new TH2F("hInvMassBgk_vs_Pt","hInvMassBgk_vs_Pt",480,0.,12.,80,0.,20.);
  85 TH1F *hInvMassRes;
  86
  87   if (ResType == 553) {
  88     hInvMassRes = new TH1F("hInvMassRes", "Mu+Mu- invariant mass (GeV/c2) around Upsilon", 60, 8., 11.);
  89   } else {
  90     hInvMassRes = new TH1F("hInvMassRes", "Mu+Mu- invariant mass (GeV/c2) around J/Psi", 80, 0., 5.);
  91   }
  92
  93   TH1F *hNumberOfTrack = new TH1F("hNumberOfTrack","nb of track /evt ",20,-0.5,19.5);
  94   TH1F *hRapMuon = new TH1F("hRapMuon"," Muon Rapidity",50,-4.5,-2);
  95   TH1F *hRapResonance = new TH1F("hRapResonance"," Resonance Rapidity",50,-4.5,-2);
  96   TH1F *hPtResonance = new TH1F("hPtResonance", "Resonance Pt (GeV/c)", 100, 0., 20.);
  97   TH2F *hThetaPhiPlus = new TH2F("hThetaPhiPlus", "Theta vs Phi +", 760, -190., 190., 400, 160., 180.);
  98   TH2F *hThetaPhiMinus = new TH2F("hThetaPhiMinus", "Theta vs Phi -", 760, -190., 190., 400, 160., 180.);
  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 thetaX, thetaY, pYZ;
 108   Double_t fPxRec1, fPyRec1, fPzRec1, fE1;
 109   Double_t fPxRec2, fPyRec2, fPzRec2, fE2;
 110   Int_t fCharge, fCharge2;
 111
 112   Int_t ntrackhits, nevents;
 113   Double_t fitfmin;
 114
 115
 116   TLorentzVector fV1, fV2, fVtot;
 117
 118   // open run loader and load gAlice, kinematics and header
 119   AliRunLoader* runLoader = AliRunLoader::Open(filename);
 120   if (!runLoader) {
 121     Error("MUONmass_ESD", "getting run loader from file %s failed",
 122             filename);
 123     return kFALSE;
 124   }
 125
 126   runLoader->LoadgAlice();
 127   gAlice = runLoader->GetAliRun();
 128   if (!gAlice) {
 129     Error("MUONmass_ESD", "no galice object found");
 130     return kFALSE;
 131   }
 132
 133
 134   // open the ESD file
 135   TFile* esdFile = TFile::Open(esdFileName);
 136   if (!esdFile || !esdFile->IsOpen()) {
 137     Error("MUONmass_ESD", "opening ESD file %s failed", esdFileName);
 138     return kFALSE;
 139   }
 140
 141   runLoader->LoadHeader();
 142   nevents = runLoader->GetNumberOfEvents();
 143
 144   // Loop over events
 145   for (Int_t iEvent = FirstEvent; iEvent <= TMath::Min(LastEvent, nevents - 1); iEvent++) {
 146
 147     // get current event
 148     runLoader->GetEvent(iEvent);
 149
 150     // get the event summary data
 151     char esdName[256];
 152     sprintf(esdName, "ESD%d", iEvent);
 153     AliESD* esd = (AliESD*) esdFile->Get(esdName);
 154     if (!esd) {
 155       Error("MUONmass_ESD", "no ESD object found for event %d", iEvent);
 156       return kFALSE;
 157     }
 158
 159     Int_t nTracks = (Int_t)esd->GetNumberOfMuonTracks() ; //
 160
 161     //    printf("\n Nb of events analysed: %d\r",iEvent);
 162     //   cout << " number of tracks: " << nrectracks  <<endl;
 163
 164     // loop over all reconstructed tracks (also first track of combination)
 165     for (Int_t iTrack = 0; iTrack <  nTracks;  iTrack++) {
 166
 167       AliESDMuonTrack* muonTrack = esd->GetMuonTrack(iTrack);
 168
 169       thetaX = muonTrack->GetThetaX();
 170       thetaY = muonTrack->GetThetaY();
 171
 172       pYZ     =  1./TMath::Abs(muonTrack->GetInverseBendingMomentum());
 173       fPzRec1  = - pYZ / TMath::Sqrt(1.0 + TMath::Tan(thetaY)*TMath::Tan(thetaX));
 174       fPxRec1  = fPzRec1 * TMath::Tan(thetaX);
 175       fPyRec1  = fPzRec1 * TMath::Tan(thetaY);
 176       fCharge = Int_t(TMath::Sign(1.,muonTrack->GetInverseBendingMomentum()));
 177
 178       fE1 = TMath::Sqrt(muonMass * muonMass + fPxRec1 * fPxRec1 + fPyRec1 * fPyRec1 + fPzRec1 * fPzRec1);
 179       fV1.SetPxPyPzE(fPxRec1, fPyRec1, fPzRec1, fE1);
 180
 181       ntrackhits = muonTrack->GetNHit();
 182       fitfmin    = muonTrack->GetChi2();
 183
 184       // transverse momentum
 185       Float_t pt1 = fV1.Pt();
 186
 187       // total momentum
 188       Float_t p1 = fV1.P();
 189
 190       // Rapidity
 191       Float_t rapMuon1 = fV1.Rapidity();
 192
 193       // chi2 per d.o.f.
 194       Float_t ch1 =  fitfmin / (2.0 * ntrackhits - 5);
 195 //      printf(" px %f py %f pz %f NHits %d  Norm.chi2 %f charge %d\n",
 196 //           fPxRec1, fPyRec1, fPzRec1, ntrackhits, ch1, fCharge);
 197
 198       // condition for good track (Chi2Cut and PtCut)
 199
 200       if ((ch1 < Chi2Cut) && (pt1 > PtCutMin) && (pt1 < PtCutMax)) {
 201
 202         // fill histos hPtMuon and hChi2PerDof
 203         hPtMuon->Fill(pt1);
 204         hPMuon->Fill(p1);
 205         hChi2PerDof->Fill(ch1);
 206         hRapMuon->Fill(rapMuon1);
 207         if (fCharge > 0) {
 208           hPtMuonPlus->Fill(pt1);
 209           hThetaPhiPlus->Fill(TMath::ATan2(fPyRec1,fPxRec1)*180./TMath::Pi(),TMath::ATan2(pt1,fPzRec1)*180./3.1415);
 210         } else {
 211           hPtMuonMinus->Fill(pt1);
 212           hThetaPhiMinus->Fill(TMath::ATan2(fPyRec1,fPxRec1)*180./TMath::Pi(),TMath::ATan2(pt1,fPzRec1)*180./3.1415);
 213         }
 214         // loop over second track of combination
 215         for (Int_t iTrack2 = iTrack + 1; iTrack2 < nTracks; iTrack2++) {
 216
 217           AliESDMuonTrack* muonTrack = esd->GetMuonTrack(iTrack2);
 218
 219           thetaX = muonTrack->GetThetaX();
 220           thetaY = muonTrack->GetThetaY();
 221
 222           pYZ     =  1./TMath::Abs(muonTrack->GetInverseBendingMomentum());
 223           fPzRec2  = - pYZ / TMath::Sqrt(1.0 + TMath::Tan(thetaY)*TMath::Tan(thetaX));
 224           fPxRec2  = fPzRec2 * TMath::Tan(thetaX);
 225           fPyRec2  = fPzRec2 * TMath::Tan(thetaY);
 226           fCharge2 = Int_t(TMath::Sign(1.,muonTrack->GetInverseBendingMomentum()));
 227
 228           fE2 = TMath::Sqrt(muonMass * muonMass + fPxRec2 * fPxRec2 + fPyRec2 * fPyRec2 + fPzRec2 * fPzRec2);
 229           fV2.SetPxPyPzE(fPxRec2, fPyRec2, fPzRec2, fE2);
 230
 231           ntrackhits = muonTrack->GetNHit();
 232           fitfmin    = muonTrack->GetChi2();
 233
 234           // transverse momentum
 235           Float_t pt2 = fV2.Pt();
 236
 237           // chi2 per d.o.f.
 238           Float_t ch2 = fitfmin  / (2.0 * ntrackhits - 5);
 239
 240           // condition for good track (Chi2Cut and PtCut)
 241           if ((ch2 < Chi2Cut) && (pt2 > PtCutMin)  && (pt2 < PtCutMax)) {
 242
 243             // condition for opposite charges
 244             if ((fCharge * fCharge2) == -1) {
 245
 246               // invariant mass
 247               fVtot = fV1 + fV2;
 248               Float_t invMass = fVtot.M();
 249
 250               // fill histos hInvMassAll and hInvMassRes
 251               hInvMassAll->Fill(invMass);
 252               hInvMassRes->Fill(invMass);
 253               hInvMassAll_vs_Pt->Fill(invMass,fVtot.Pt());
 254               if (invMass > massMin && invMass < massMax) {
 255                 EventInMass++;
 256                 hRapResonance->Fill(fVtot.Rapidity());
 257                 hPtResonance->Fill(fVtot.Pt());
 258               }
 259
 260             } // if (fCharge * fCharge2) == -1)
 261           } // if ((ch2 < Chi2Cut) && (pt2 > PtCutMin) && (pt2 < PtCutMax))
 262         } //  for (Int_t iTrack2 = iTrack + 1; iTrack2 < iTrack; iTrack2++)
 263       } // if (ch1 < Chi2Cut) && (pt1 > PtCutMin)&& (pt1 < PtCutMax) )
 264     } // for (Int_t iTrack = 0; iTrack < nrectracks; iTrack++)
 265
 266     hNumberOfTrack->Fill(nTracks);
 267     esdFile->Delete();
 268   } // for (Int_t iEvent = FirstEvent;
 269
 270 // Loop over events for bg event
 271
 272   Double_t thetaPlus,  phiPlus;
 273   Double_t thetaMinus, phiMinus;
 274   Float_t PtMinus, PtPlus;
 275
 276   for (Int_t iEvent = 0; iEvent < hInvMassAll->Integral(); iEvent++) {
 277
 278     hThetaPhiPlus->GetRandom2(phiPlus, thetaPlus);
 279     hThetaPhiMinus->GetRandom2(phiMinus,thetaMinus);
 280     PtPlus = hPtMuonPlus->GetRandom();
 281     PtMinus = hPtMuonMinus->GetRandom();
 282
 283     fPxRec1  = PtPlus * TMath::Cos(TMath::Pi()/180.*phiPlus);
 284     fPyRec1  = PtPlus * TMath::Sin(TMath::Pi()/180.*phiPlus);
 285     fPzRec1  = PtPlus / TMath::Tan(TMath::Pi()/180.*thetaPlus);
 286
 287     fE1 = TMath::Sqrt(muonMass * muonMass + fPxRec1 * fPxRec1 + fPyRec1 * fPyRec1 + fPzRec1 * fPzRec1);
 288     fV1.SetPxPyPzE(fPxRec1, fPyRec1, fPzRec1, fE1);
 289
 290     fPxRec2  = PtMinus * TMath::Cos(TMath::Pi()/180.*phiMinus);
 291     fPyRec2  = PtMinus * TMath::Sin(TMath::Pi()/180.*phiMinus);
 292     fPzRec2  = PtMinus / TMath::Tan(TMath::Pi()/180.*thetaMinus);
 293
 294     fE2 = TMath::Sqrt(muonMass * muonMass + fPxRec2 * fPxRec2 + fPyRec2 * fPyRec2 + fPzRec2 * fPzRec2);
 295     fV2.SetPxPyPzE(fPxRec2, fPyRec2, fPzRec2, fE2);
 296
 297     // invariant mass
 298     fVtot = fV1 + fV2;
 299
 300     // fill histos hInvMassAll and hInvMassRes
 301     hInvMassBg->Fill(fVtot.M());
 302     hInvMassBgk_vs_Pt->Fill( fVtot.M(), fVtot.Pt() );
 303   }
 304
 305   histoFile->Write();
 306   histoFile->Close();
 307
 308   cout << "MUONmassPlot " << endl;
 309   cout << "FirstEvent " << FirstEvent << endl;
 310   cout << "LastEvent " << LastEvent << endl;
 311   cout << "ResType " << ResType << endl;
 312   cout << "Chi2Cut " << Chi2Cut << endl;
 313   cout << "PtCutMin " << PtCutMin << endl;
 314   cout << "PtCutMax " << PtCutMax << endl;
 315   cout << "massMin " << massMin << endl;
 316   cout << "massMax " << massMax << endl;
 317   cout << "EventInMass " << EventInMass << endl;
 318
 319   return kTRUE;
 320 }
 321