Make the Scan method public
[u/mrichter/AliRoot.git] / MUON / MUONmassPlot_ESD.C
index c507878..5b30efc 100644 (file)
 #include "TParticle.h"
 #include "TTree.h"
 #include <Riostream.h>
+#include <TGeoManager.h>
+#include <TROOT.h>
 
 // STEER includes
-#include "AliRun.h"
-#include "AliRunLoader.h"
-#include "AliHeader.h"
-#include "AliLoader.h"
-#include "AliStack.h"
-#include "AliMagF.h"
-#include "AliESD.h"
+#include "AliLog.h"
+#include "AliCDBManager.h"
+#include "AliESDEvent.h"
+#include "AliESDVertex.h"
+#include "AliESDMuonTrack.h"
 
 // MUON includes
-#include "AliESDMuonTrack.h"
+#include "AliMUONCDB.h"
+#include "AliMUONTrackParam.h"
+#include "AliMUONTrackExtrap.h"
+#include "AliMUONESDInterface.h"
 #endif
-//
-// Macro MUONmassPlot.C for ESD
-// Ch. Finck, Subatech, April. 2004
-//
-
-// macro to make invariant mass plots
-// for combinations of 2 muons with opposite charges,
-// from root file "MUON.tracks.root" containing the result of track reconstruction.
-// Histograms are stored on the "MUONmassPlot.root" file.
-// introducing TLorentzVector for parameter calculations (Pt, P,rap,etc...)
-// using Invariant Mass for rapidity.
-
-// Arguments:
-//   FirstEvent (default 0)
-//   LastEvent (default 0)
-//   ResType (default 553)
-//      553 for Upsilon, anything else for J/Psi
-//   Chi2Cut (default 100)
-//      to keep only tracks with chi2 per d.o.f. < Chi2Cut
-//   PtCutMin (default 1)
-//      to keep only tracks with transverse momentum > PtCutMin
-//   PtCutMax (default 10000)
-//      to keep only tracks with transverse momentum < PtCutMax
-//   massMin (default 9.17 for Upsilon) 
-//      &  massMax (default 9.77 for Upsilon) 
-//         to calculate the reconstruction efficiency for resonances with invariant mass
-//         massMin < mass < massMax.
-
-// Add parameters and histograms for analysis 
-
-Bool_t MUONmassPlot(char* filename = "galice.root", Int_t FirstEvent = 0, Int_t LastEvent = 10000,
-                 char* esdFileName = "AliESDs.root", 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)
+
+/// \ingroup macros
+/// \file MUONmassPlot_ESD.C
+/// \brief Macro MUONmassPlot_ESD.C for ESD
+///
+/// \author Ch. Finck, Subatech, April. 2004
+///
+///
+/// Macro to make invariant mass plots
+/// for combinations of 2 muons with opposite charges,
+/// from root file "MUON.tracks.root" containing the result of track reconstruction.
+/// Histograms are stored on the "MUONmassPlot.root" file.
+/// introducing TLorentzVector for parameter calculations (Pt, P,rap,etc...)
+/// using Invariant Mass for rapidity.
+///
+/// Add parameters and histograms for analysis 
+
+Bool_t MUONmassPlot(const char* esdFileName = "AliESDs.root", const char* geoFilename = "geometry.root",
+                   const char* ocdbPath = "local://$ALICE_ROOT/OCDB",
+                   Int_t FirstEvent = 0, Int_t LastEvent = 10000, Int_t ExtrapToVertex = -1,
+                   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)
 {
+  /// \param FirstEvent (default 0)
+  /// \param LastEvent (default 10000)
+  /// \param ExtrapToVertex (default -1)
+  ///   -      <0: no extrapolation;
+  ///   -      =0: extrapolation to (0,0,0);
+  ///   -      >0: extrapolation to ESDVertex if available, else to (0,0,0)
+  /// \param ResType    553 for Upsilon, anything else for J/Psi (default 553)
+  /// \param Chi2Cut    to keep only tracks with chi2 per d.o.f. < Chi2Cut (default 100)
+  /// \param PtCutMin   to keep only tracks with transverse momentum > PtCutMin (default 1)
+  /// \param PtCutMax   to keep only tracks with transverse momentum < PtCutMax (default 10000)
+  /// \param massMin   (default 9.17 for Upsilon) 
+  /// \param massMax   (default 9.77 for Upsilon);  
+  ///         to calculate the reconstruction efficiency for resonances with invariant mass
+  ///         massMin < mass < massMax.
+
   cout << "MUONmassPlot " << endl;
   cout << "FirstEvent " << FirstEvent << endl;
   cout << "LastEvent " << LastEvent << endl;
@@ -81,9 +87,10 @@ Bool_t MUONmassPlot(char* filename = "galice.root", Int_t FirstEvent = 0, Int_t
   TH1F *hChi2PerDof = new TH1F("hChi2PerDof", "Muon track chi2/d.o.f.", 100, 0., 20.);
   TH1F *hInvMassAll = new TH1F("hInvMassAll", "Mu+Mu- invariant mass (GeV/c2)", 480, 0., 12.);
   TH1F *hInvMassBg = new TH1F("hInvMassBg", "Mu+Mu- invariant mass BG(GeV/c2)", 480, 0., 12.);
-TH2F *hInvMassAll_vs_Pt = new TH2F("hInvMassAll_vs_Pt","hInvMassAll_vs_Pt",480,0.,12.,80,0.,20.);
-TH2F *hInvMassBgk_vs_Pt = new TH2F("hInvMassBgk_vs_Pt","hInvMassBgk_vs_Pt",480,0.,12.,80,0.,20.);
-TH1F *hInvMassRes;
+  TH2F *hInvMassAll_vs_Pt = new TH2F("hInvMassAll_vs_Pt","hInvMassAll_vs_Pt",480,0.,12.,80,0.,20.);
+  TH2F *hInvMassBgk_vs_Pt = new TH2F("hInvMassBgk_vs_Pt","hInvMassBgk_vs_Pt",480,0.,12.,80,0.,20.);
+  TH1F *hInvMassRes;
+  TH1F *hPrimaryVertex = new TH1F("hPrimaryVertex","SPD reconstructed Z vertex",150,-15,15);
 
   if (ResType == 553) {
     hInvMassRes = new TH1F("hInvMassRes", "Mu+Mu- invariant mass (GeV/c2) around Upsilon", 60, 8., 11.);
@@ -101,68 +108,87 @@ TH1F *hInvMassRes;
 
   // settings
   Int_t EventInMass = 0;
+  Int_t EventInMassMatch = 0;
+  Int_t NbTrigger = 0;
+
   Float_t muonMass = 0.105658389;
 //   Float_t UpsilonMass = 9.46037;
 //   Float_t JPsiMass = 3.097;
 
-  Double_t thetaX, thetaY, pYZ;
+  Int_t fCharge1, fCharge2;
   Double_t fPxRec1, fPyRec1, fPzRec1, fE1;
   Double_t fPxRec2, fPyRec2, fPzRec2, fE2;
-  Int_t fCharge, fCharge2;
 
-  Int_t ntrackhits, nevents;
+  Int_t ntrackhits;
   Double_t fitfmin;
-
+  Double_t fZVertex=0;
+  Double_t fYVertex=0;
+  Double_t fXVertex=0;
+  Double_t errXVtx=0;
+  Double_t errYVtx=0;
  
   TLorentzVector fV1, fV2, fVtot;
-
-  // set off mag field 
-  AliMagF::SetReadField(kFALSE);
-
-  // open run loader and load gAlice, kinematics and header
-  AliRunLoader* runLoader = AliRunLoader::Open(filename);
-  if (!runLoader) {
-    Error("MUONmass_ESD", "getting run loader from file %s failed", 
-           filename);
-    return kFALSE;
-  }
-
-  if (!gAlice) {
-    Error("MUONmass_ESD", "no galice object found");
-    return kFALSE;
+  
+  // Import TGeo geometry (needed by AliMUONTrackExtrap::ExtrapToVertex)
+  if (!gGeoManager) {
+    TGeoManager::Import(geoFilename);
+    if (!gGeoManager) {
+      Error("MUONmass_ESD", "getting geometry from file %s failed", geoFilename);
+      return kFALSE;
+    }
   }
   
-
   // open the ESD file
   TFile* esdFile = TFile::Open(esdFileName);
   if (!esdFile || !esdFile->IsOpen()) {
     Error("MUONmass_ESD", "opening ESD file %s failed", esdFileName);
     return kFALSE;
   }
-  
-  AliESD* esd = new AliESD();
+  AliESDEvent* esd = new AliESDEvent();
   TTree* tree = (TTree*) esdFile->Get("esdTree");
   if (!tree) {
-    Error("CheckESD", "no ESD tree found");
+    Error("MUONmass_ESD", "no ESD tree found");
+    return kFALSE;
+  }
+  esd->ReadFromTree(tree);
+  
+  // get run number
+  if (tree->GetEvent(0) <= 0) {
+    Error("MUONmass_ESD", "no ESD object found for event 0");
     return kFALSE;
   }
-  tree->SetBranchAddress("ESD", &esd);
+  Int_t runNumber = esd->GetRunNumber();
+  
+  // load necessary data from OCDB
+  AliCDBManager::Instance()->SetDefaultStorage(ocdbPath);
+  AliCDBManager::Instance()->SetRun(runNumber);
+  if (!AliMUONCDB::LoadField()) return kFALSE;
+  
+  // set the magnetic field for track extrapolations
+  AliMUONTrackExtrap::SetField();
+  
+  AliMUONTrackParam trackParam;
 
-  runLoader->LoadHeader();
-  nevents = runLoader->GetNumberOfEvents();
-        
   // Loop over events
+  Int_t nevents = (Int_t)tree->GetEntries();
   for (Int_t iEvent = FirstEvent; iEvent <= TMath::Min(LastEvent, nevents - 1); iEvent++) {
 
-    // get current event
-    runLoader->GetEvent(iEvent);
-   
     // get the event summary data
-    tree->GetEvent(iEvent);
-    if (!esd) {
-      Error("CheckESD", "no ESD object found for event %d", iEvent);
+    if (tree->GetEvent(iEvent) <= 0) {
+      Error("MUONmass_ESD", "no ESD object found for event %d", iEvent);
       return kFALSE;
     }
+    
+    // get the SPD reconstructed vertex (vertexer) and fill the histogram
+    AliESDVertex* Vertex = (AliESDVertex*) esd->GetVertex();
+    if (Vertex->GetNContributors()) {
+      fZVertex = Vertex->GetZv();
+      fYVertex = Vertex->GetYv();
+      fXVertex = Vertex->GetXv();
+      errXVtx = Vertex->GetXRes();
+      errYVtx = Vertex->GetYRes();
+    }
+    hPrimaryVertex->Fill(fZVertex);
 
     Int_t nTracks = (Int_t)esd->GetNumberOfMuonTracks() ; 
 
@@ -172,20 +198,26 @@ TH1F *hInvMassRes;
     // loop over all reconstructed tracks (also first track of combination)
     for (Int_t iTrack = 0; iTrack <  nTracks;  iTrack++) {
 
-      AliESDMuonTrack* muonTrack = esd->GetMuonTrack(iTrack);
-
-      thetaX = muonTrack->GetThetaX();
-      thetaY = muonTrack->GetThetaY();
-
-      pYZ     =  1./TMath::Abs(muonTrack->GetInverseBendingMomentum());
-      fPzRec1  = - pYZ / TMath::Sqrt(1.0 + TMath::Tan(thetaY)*TMath::Tan(thetaY));
-      fPxRec1  = fPzRec1 * TMath::Tan(thetaX);
-      fPyRec1  = fPzRec1 * TMath::Tan(thetaY);
-      fCharge = Int_t(TMath::Sign(1.,muonTrack->GetInverseBendingMomentum()));
-
-      fE1 = TMath::Sqrt(muonMass * muonMass + fPxRec1 * fPxRec1 + fPyRec1 * fPyRec1 + fPzRec1 * fPzRec1);
-      fV1.SetPxPyPzE(fPxRec1, fPyRec1, fPzRec1, fE1);
-
+      // skip ghosts
+      if (!esd->GetMuonTrack(iTrack)->ContainTrackerData()) continue;
+      
+      AliESDMuonTrack* muonTrack = new AliESDMuonTrack(*(esd->GetMuonTrack(iTrack)));
+
+      // extrapolate to vertex if required and available
+      if (ExtrapToVertex > 0 && Vertex->GetNContributors()) {
+       AliMUONESDInterface::GetParamAtFirstCluster(*muonTrack, trackParam);
+       AliMUONTrackExtrap::ExtrapToVertex(&trackParam, fXVertex, fYVertex, fZVertex, errXVtx, errYVtx);
+       AliMUONESDInterface::SetParamAtVertex(trackParam, *muonTrack); // put the new parameters in this copy of AliESDMuonTrack
+      } else if ((ExtrapToVertex > 0 && !Vertex->GetNContributors()) || ExtrapToVertex == 0){
+       AliMUONESDInterface::GetParamAtFirstCluster(*muonTrack, trackParam);
+       AliMUONTrackExtrap::ExtrapToVertex(&trackParam, 0., 0., 0., 0., 0.);
+       AliMUONESDInterface::SetParamAtVertex(trackParam, *muonTrack); // put the new parameters in this copy of AliESDMuonTrack
+      }
+
+      fCharge1 = Int_t(TMath::Sign(1.,muonTrack->GetInverseBendingMomentum()));
+      
+      muonTrack->LorentzP(fV1);
+      
       ntrackhits = muonTrack->GetNHit();
       fitfmin    = muonTrack->GetChi2();
 
@@ -201,43 +233,49 @@ TH1F *hInvMassRes;
       // chi2 per d.o.f.
       Float_t ch1 =  fitfmin / (2.0 * ntrackhits - 5);
 //      printf(" px %f py %f pz %f NHits %d  Norm.chi2 %f charge %d\n", 
-//          fPxRec1, fPyRec1, fPzRec1, ntrackhits, ch1, fCharge);
+//          fPxRec1, fPyRec1, fPzRec1, ntrackhits, ch1, fCharge1);
 
       // condition for good track (Chi2Cut and PtCut)
 
-      if ((ch1 < Chi2Cut) && (pt1 > PtCutMin) && (pt1 < PtCutMax)) {
+//      if ((ch1 < Chi2Cut) && (pt1 > PtCutMin) && (pt1 < PtCutMax)) {
 
        // fill histos hPtMuon and hChi2PerDof
        hPtMuon->Fill(pt1);
        hPMuon->Fill(p1);
        hChi2PerDof->Fill(ch1);
        hRapMuon->Fill(rapMuon1);
-       if (fCharge > 0) {
+       if (fCharge1 > 0) {
          hPtMuonPlus->Fill(pt1);
-         hThetaPhiPlus->Fill(TMath::ATan2(fPyRec1,fPxRec1)*180./TMath::Pi(),TMath::ATan2(pt1,fPzRec1)*180./3.1415);
+         hThetaPhiPlus->Fill(fV1.Phi()*180./TMath::Pi(),fV1.Theta()*180./TMath::Pi());
        } else {
          hPtMuonMinus->Fill(pt1);
-         hThetaPhiMinus->Fill(TMath::ATan2(fPyRec1,fPxRec1)*180./TMath::Pi(),TMath::ATan2(pt1,fPzRec1)*180./3.1415);
+         hThetaPhiMinus->Fill(fV1.Phi()*180./TMath::Pi(),fV1.Theta()*180./TMath::Pi());
        }
        // loop over second track of combination
        for (Int_t iTrack2 = iTrack + 1; iTrack2 < nTracks; iTrack2++) {
          
-         AliESDMuonTrack* muonTrack = esd->GetMuonTrack(iTrack2);
-
-         thetaX = muonTrack->GetThetaX();
-         thetaY = muonTrack->GetThetaY();
-
-         pYZ     =  1./TMath::Abs(muonTrack->GetInverseBendingMomentum());
-         fPzRec2  = - pYZ / TMath::Sqrt(1.0 + TMath::Tan(thetaY)*TMath::Tan(thetaY));
-         fPxRec2  = fPzRec2 * TMath::Tan(thetaX);
-         fPyRec2  = fPzRec2 * TMath::Tan(thetaY);
-         fCharge2 = Int_t(TMath::Sign(1.,muonTrack->GetInverseBendingMomentum()));
+         // skip ghosts
+         if (!esd->GetMuonTrack(iTrack2)->ContainTrackerData()) continue;
+         
+         AliESDMuonTrack* muonTrack2 = new AliESDMuonTrack(*(esd->GetMuonTrack(iTrack2)));
+         
+         // extrapolate to vertex if required and available
+         if (ExtrapToVertex > 0 && Vertex->GetNContributors()) {
+           AliMUONESDInterface::GetParamAtFirstCluster(*muonTrack2, trackParam);
+           AliMUONTrackExtrap::ExtrapToVertex(&trackParam, fXVertex, fYVertex, fZVertex, errXVtx, errYVtx);
+           AliMUONESDInterface::SetParamAtVertex(trackParam, *muonTrack2); // put the new parameters in this copy of AliESDMuonTrack
+         } else if ((ExtrapToVertex > 0 && !Vertex->GetNContributors()) || ExtrapToVertex == 0){
+           AliMUONESDInterface::GetParamAtFirstCluster(*muonTrack2, trackParam);
+           AliMUONTrackExtrap::ExtrapToVertex(&trackParam, 0., 0., 0., 0., 0.);
+           AliMUONESDInterface::SetParamAtVertex(trackParam, *muonTrack2); // put the new parameters in this copy of AliESDMuonTrack
+         }
+         
+         fCharge2 = Int_t(TMath::Sign(1.,muonTrack2->GetInverseBendingMomentum()));
 
-         fE2 = TMath::Sqrt(muonMass * muonMass + fPxRec2 * fPxRec2 + fPyRec2 * fPyRec2 + fPzRec2 * fPzRec2);
-         fV2.SetPxPyPzE(fPxRec2, fPyRec2, fPzRec2, fE2);
+         muonTrack2->LorentzP(fV2);
 
-         ntrackhits = muonTrack->GetNHit();
-         fitfmin    = muonTrack->GetChi2();
+         ntrackhits = muonTrack2->GetNHit();
+         fitfmin    = muonTrack2->GetChi2();
 
          // transverse momentum
          Float_t pt2 = fV2.Pt();
@@ -246,10 +284,10 @@ TH1F *hInvMassRes;
          Float_t ch2 = fitfmin  / (2.0 * ntrackhits - 5);
 
          // condition for good track (Chi2Cut and PtCut)
-         if ((ch2 < Chi2Cut) && (pt2 > PtCutMin)  && (pt2 < PtCutMax)) {
+//       if ((ch2 < Chi2Cut) && (pt2 > PtCutMin)  && (pt2 < PtCutMax)) {
 
            // condition for opposite charges
-           if ((fCharge * fCharge2) == -1) {
+           if ((fCharge1 * fCharge2) == -1) {
 
              // invariant mass
              fVtot = fV1 + fV2;
@@ -259,16 +297,28 @@ TH1F *hInvMassRes;
              hInvMassAll->Fill(invMass);
              hInvMassRes->Fill(invMass);
              hInvMassAll_vs_Pt->Fill(invMass,fVtot.Pt());
+             Int_t ptTrig;
+             if (ResType == 553) 
+               ptTrig =  0x20;// mask for Hpt unlike sign pair
+             else 
+               ptTrig =  0x10;// mask for Lpt unlike sign pair
+
+             if (esd->GetTriggerMask() &  ptTrig) NbTrigger++; 
              if (invMass > massMin && invMass < massMax) {
                EventInMass++;
+               if (muonTrack2->GetMatchTrigger() && (esd->GetTriggerMask() & ptTrig))// match with trigger
+                 EventInMassMatch++;
+
                hRapResonance->Fill(fVtot.Rapidity());
                hPtResonance->Fill(fVtot.Pt());
              }
 
-           } // if (fCharge * fCharge2) == -1)
-         } // if ((ch2 < Chi2Cut) && (pt2 > PtCutMin) && (pt2 < PtCutMax))
+           } // if (fCharge1 * fCharge2) == -1)
+//       } // if ((ch2 < Chi2Cut) && (pt2 > PtCutMin) && (pt2 < PtCutMax))
+         delete muonTrack2;
        } //  for (Int_t iTrack2 = iTrack + 1; iTrack2 < iTrack; iTrack2++)
-      } // if (ch1 < Chi2Cut) && (pt1 > PtCutMin)&& (pt1 < PtCutMax) )
+//      } // if (ch1 < Chi2Cut) && (pt1 > PtCutMin)&& (pt1 < PtCutMax) )
+      delete muonTrack;
     } // for (Int_t iTrack = 0; iTrack < nrectracks; iTrack++)
 
     hNumberOfTrack->Fill(nTracks);
@@ -313,16 +363,10 @@ TH1F *hInvMassRes;
   histoFile->Write();
   histoFile->Close();
 
-  cout << "MUONmassPlot " << endl;
-  cout << "FirstEvent " << FirstEvent << endl;
-  cout << "LastEvent " << LastEvent << endl;
-  cout << "ResType " << ResType << endl;
-  cout << "Chi2Cut " << Chi2Cut << endl;
-  cout << "PtCutMin " << PtCutMin << endl;
-  cout << "PtCutMax " << PtCutMax << endl;
-  cout << "massMin " << massMin << endl;
-  cout << "massMax " << massMax << endl;
+  cout << endl;
   cout << "EventInMass " << EventInMass << endl;
+  cout << "NbTrigger " << NbTrigger << endl;
+  cout << "EventInMass match with trigger " << EventInMassMatch << endl;
 
   return kTRUE;
 }