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 //   NSigma (default 3)
  18 //      the number of combinations is counted around the resonance mass
  19 //      within +/- NSigma times the nominal sigma's
  20 //      (0.099 GeV for Upsilon, 0.0615 GeV for J/Psi)
  21 //   Chi2Cut (default 100)
  22 //      to keep only tracks with chi2 per d.o.f. < Chi2Cut
  23 //   PtCut (default 1)
  24 //      to keep only tracks with transverse momentum > PtCut
  25
  26 // IMPORTANT NOTICE FOR USERS:
  27 // under "root" or "root.exe", execute the following commands:
  28 // 1. "gSystem->SetIncludePath("-I$ALICE_ROOT/MUON -I$ALICE_ROOT/STEER")" to get the right path at compilation time
  29 // 2. ".x loadlibs.C" to load the shared libraries
  30 // 3. ".L MUONrecoNtuple.C+"
  31 // 4. ".x MUONmassPlot.C()" with the right arguments according to the list above
  32
  33 void MUONmassPlot(Int_t FirstEvent = 0, Int_t LastEvent = 0, Int_t ResType = 553, Float_t Nsig = 3., Float_t Chi2Cut = 100., Float_t PtCut = 1.)
  34 {
  35   cout << "MUONmassPlot" << endl;
  36   cout << "FirstEvent" << FirstEvent << endl;
  37   cout << "LastEvent" << LastEvent << endl;
  38   cout << "ResType" << ResType << endl;
  39   cout << "Nsig" << Nsig << endl;
  40   cout << "Chi2Cut" << Chi2Cut << endl;
  41   cout << "PtCut" << PtCut << endl;
  42
  43   //////////////////////////////////////////////////////////
  44   //   This file has been automatically generated
  45   //     (Thu Sep 21 14:53:11 2000 by ROOT version2.25/02)
  46   //   from TTree MUONtrackReco/MUONtrackReco
  47   //   found on file: MUONtrackReco.root
  48   //////////////////////////////////////////////////////////
  49
  50
  51   //Reset ROOT and connect tree file
  52   gROOT->Reset();
  53   TFile *f = (TFile*)gROOT->GetListOfFiles()->FindObject("MUONtrackReco.root");
  54   if (!f) {
  55     f = new TFile("MUONtrackReco.root");
  56   }
  57   TTree *MUONtrackReco = (TTree*)gDirectory->Get("MUONtrackReco");
  58   MUONtrackReco->SetMakeClass(1);
  59
  60 //Declaration of leaves types
  61   Int_t           fEvent;
  62   UInt_t          fUniqueID;
  63   UInt_t          fBits;
  64   Int_t           Tracks_;
  65   Int_t           Tracks_fCharge[5];
  66   Float_t         Tracks_fPxRec[5];
  67   Float_t         Tracks_fPyRec[5];
  68   Float_t         Tracks_fPzRec[5];
  69   Float_t         Tracks_fZRec[5];
  70   Float_t         Tracks_fZRec1[5];
  71   Int_t           Tracks_fNHits[5];
  72   Float_t         Tracks_fChi2[5];
  73   Float_t         Tracks_fPxGen[5];
  74   Float_t         Tracks_fPyGen[5];
  75   Float_t         Tracks_fPzGen[5];
  76   UInt_t          Tracks_fUniqueID[5];
  77   UInt_t          Tracks_fBits[5];
  78
  79   //Set branch addresses
  80   //MUONtrackReco->SetBranchAddress("Header",&Header);
  81   MUONtrackReco->SetBranchAddress("fEvent",&fEvent);
  82   MUONtrackReco->SetBranchAddress("fUniqueID",&fUniqueID);
  83   MUONtrackReco->SetBranchAddress("fBits",&fBits);
  84 //   MUONtrackReco->SetBranchAddress("Tracks_",&Tracks_);
  85   MUONtrackReco->SetBranchAddress("Tracks",&Tracks_);
  86   MUONtrackReco->SetBranchAddress("Tracks.fCharge",Tracks_fCharge);
  87   MUONtrackReco->SetBranchAddress("Tracks.fPxRec",Tracks_fPxRec);
  88   MUONtrackReco->SetBranchAddress("Tracks.fPyRec",Tracks_fPyRec);
  89   MUONtrackReco->SetBranchAddress("Tracks.fPzRec",Tracks_fPzRec);
  90   MUONtrackReco->SetBranchAddress("Tracks.fZRec",Tracks_fZRec);
  91   MUONtrackReco->SetBranchAddress("Tracks.fZRec1",Tracks_fZRec1);
  92   MUONtrackReco->SetBranchAddress("Tracks.fNHits",Tracks_fNHits);
  93   MUONtrackReco->SetBranchAddress("Tracks.fChi2",Tracks_fChi2);
  94   MUONtrackReco->SetBranchAddress("Tracks.fPxGen",Tracks_fPxGen);
  95   MUONtrackReco->SetBranchAddress("Tracks.fPyGen",Tracks_fPyGen);
  96   MUONtrackReco->SetBranchAddress("Tracks.fPzGen",Tracks_fPzGen);
  97   MUONtrackReco->SetBranchAddress("Tracks.fUniqueID",Tracks_fUniqueID);
  98   MUONtrackReco->SetBranchAddress("Tracks.fBits",Tracks_fBits);
  99
 100 //     This is the loop skeleton
 101 //       To read only selected branches, Insert statements like:
 102 // MUONtrackReco->SetBranchStatus("*",0);  // disable all branches
 103 // TTreePlayer->SetBranchStatus("branchname",1);  // activate branchname
 104
 105   Int_t nentries = MUONtrackReco->GetEntries();
 106
 107   Int_t nbytes = 0;
 108   //   for (Int_t i=0; i<nentries;i++) {
 109   //      nbytes += MUONtrackReco->GetEntry(i);
 110   //   }
 111
 112   /////////////////////////////////////////////////////////////////
 113   // Here comes the specialized part for MUONmassPlot
 114   /////////////////////////////////////////////////////////////////
 115
 116   // File for histograms and histogram booking
 117   TFile *histoFile = new TFile("MUONmassPlot.root", "RECREATE");
 118   TH1F *hPtMuon = new TH1F("hPtMuon", "Muon Pt (GeV/c)", 100, 0., 20.);
 119   TH1F *hChi2PerDof = new TH1F("hChi2PerDof", "Muon track chi2/d.o.f.", 100, 0., 20.);
 120   TH1F *hInvMassAll = new TH1F("hInvMassAll", "Mu+Mu- invariant mass (GeV/c2)", 240, 0., 12.);
 121   if (ResType = 553) TH1F *hInvMassRes = new TH1F("hInvMassRes", "Mu+Mu- invariant mass (GeV/c2) around Upsilon", 60, 8., 11.);
 122   else TH1F *hInvMassRes = new TH1F("hInvMassRes", "Mu+Mu- invariant mass (GeV/c2) around J/Psi", 80, 1., 5.);
 123
 124   // Loop over events
 125   for (Int_t event = FirstEvent; event <= TMath::Min(LastEvent, nentries - 1); event++) {
 126     // get current event
 127     nbytes += MUONtrackReco->GetEntry(event);
 128     // loop over all reconstructed tracks (also first track of combination)
 129     for (Int_t t1 = 0; t1 < Tracks_; t1++) {
 130       // transverse momentum
 131       Float_t pt1 = TMath::Sqrt(Tracks_fPxRec[t1] * Tracks_fPxRec[t1] + Tracks_fPyRec[t1] * Tracks_fPyRec[t1]);
 132       // chi2 per d.o.f.
 133       Float_t ch1 = Tracks_fChi2[t1] / (2.0 * Tracks_fNHits[t1] - 5);
 134       // condition for good track (Chi2Cut and PtCut)
 135       if ((ch1 < Chi2Cut) && (pt1 > PtCut)) {
 136         // fill histos hPtMuon and hChi2PerDof
 137         hPtMuon->Fill(pt1);
 138         hChi2PerDof->Fill(ch1);
 139         // loop over second track of combination
 140         for (Int_t t2 = t1 + 1; t2 < Tracks_; t2++) {
 141           // transverse momentum
 142           Float_t pt2 = TMath::Sqrt(Tracks_fPxRec[t2] * Tracks_fPxRec[t2] + Tracks_fPyRec[t2] * Tracks_fPyRec[t2]);
 143           // chi2 per d.o.f.
 144           Float_t ch2 = Tracks_fChi2[t2] / (2.0 * Tracks_fNHits[t2] - 5);
 145           // condition for good track (Chi2Cut and PtCut)
 146           if ((ch2 < Chi2Cut) && (pt2 > PtCut)) {
 147             // condition for opposite charges
 148             if ((Tracks_fCharge[t1] * Tracks_fCharge[t2]) == -1) {
 149               // invariant mass
 150               Float_t invMass = MuPlusMuMinusMass(Tracks_fPxRec[t1], Tracks_fPyRec[t1], Tracks_fPzRec[t1], Tracks_fPxRec[t2], Tracks_fPyRec[t2], Tracks_fPzRec[t2]);
 151               // fill histos hInvMassAll and hInvMassRes
 152               hInvMassAll->Fill(invMass);
 153               hInvMassRes->Fill(invMass);
 154             } //if ((Tracks_fCharge[t1] * Tracks_fCharge[t2]) == -1)
 155           } // if ((Tracks_fChi2[t2] < Chi2Cut) && (pt2 > PtCut))
 156         } // for (Int_t t2 = t1 + 1; t2 < Tracks_; t2++)
 157       } // if ((Tracks_fChi2[t1] < Chi2Cut) && (pt1 > PtCut))
 158     } // for (Int_t t1 = 0; t1 < Tracks_; t1++)
 159   } // for (Int_t event = FirstEvent;
 160
 161   histoFile->Write();
 162   histoFile->Close();
 163 }
 164
 165 Float_t MuPlusMuMinusMass(Float_t Px1, Float_t Py1, Float_t Pz1, Float_t Px2, Float_t Py2, Float_t Pz2)
 166 {
 167   Float_t muonMass = 0.10566;
 168   Float_t e1 = TMath::Sqrt(muonMass * muonMass + Px1 * Px1 + Py1 * Py1 + Pz1 * Pz1);
 169   Float_t e2 = TMath::Sqrt(muonMass * muonMass + Px2 * Px2 + Py2 * Py2 + Pz2 * Pz2);
 170   return (TMath::Sqrt(2.0 * (muonMass * muonMass + e1 * e2 - Px1 * Px2 - Py1 * Py2 - Pz1 * Pz2)));
 171 }