]> git.uio.no Git - u/mrichter/AliRoot.git/blobdiff - MUON/MUONefficiency.C
Make the Scan method public
[u/mrichter/AliRoot.git] / MUON / MUONefficiency.C
index 42b9016dc7b48a6dfebf67b9aae70d8a0fe96560..6ed1bccf300d1bd477090f7cd2777afdfe975b6e 100644 (file)
 
 /* $Id$ */
 
-// Macro (upgraded version of MUONmassPlot_ESD.C, better handling of Jpsi) to make : 
-// 1) Ntuple (Ktuple) containing Upsilon kinematics variables (from kinematics.root files) 
-// 2) Ntuple (ESDtuple) containing Upsilon kinematics variables from reconstruction and 
-// combinations of 2 muons with opposite charges (ESDtupleBck will be used later)
-// 3) Some QA histograms
-// Ntuple are stored in the file MUONefficiency.root and  ESD tree and QA histograms in AliESDs.root
+/// \ingroup macros
+/// \file MUONefficiency.C
+/// \brief add brief description
+///
+/// Macro (upgraded version of MUONmassPlot_ESD.C, better handling of Jpsi) to make : 
+/// - Ntuple (Ktuple) containing Upsilon kinematics variables (from kinematics.root files) 
+/// - Ntuple (ESDtuple) containing Upsilon kinematics variables from reconstruction and 
+///   combinations of 2 muons with opposite charges (ESDtupleBck will be used later)
+/// - Some QA histograms
+/// Ntuple are stored in the file MUONefficiency.root and  ESD tree and QA histograms in AliESDs.root
+///
+/// \author Christophe Suire, IPN Orsay
 
-// Christophe Suire, IPN Orsay
 
 
+#if !defined(__CINT__) || defined(__MAKECINT__)
 
-// Arguments:
-//   FirstEvent (default 0)
-//   LastEvent (default 1.e6)
-//   ResType (default 553)
-//      553 for Upsilon, 443 for J/Psi
-//   Chi2Cut (default 100)
-//      to keep only tracks with chi2 per d.o.f. < Chi2Cut
-
+// MUON includes
+#include "AliMUONCDB.h"
+#include "AliMUONTrackParam.h"
+#include "AliMUONTrackExtrap.h"
+#include "AliMUONESDInterface.h"
 
+// STEER includes
+#include "AliRun.h"
+#include "AliRunLoader.h"
+#include "AliHeader.h"
+#include "AliLoader.h"
+#include "AliStack.h"
+#include "AliESDEvent.h"
+#include "AliESDVertex.h"
+#include "AliCDBManager.h"
+#include "AliESDMuonTrack.h"
 
-#if !defined(__CINT__) || defined(__MAKECINT__)
 // ROOT includes
 #include "TTree.h"
 #include "TNtuple.h"
-#include "TBranch.h"
-#include "TClonesArray.h"
 #include "TLorentzVector.h"
 #include "TFile.h"
 #include "TH1.h"
 #include "TH2.h"
 #include "TParticle.h"
-#include "TTree.h"
 #include "TString.h"
 #include <Riostream.h>
+#include <TGeoManager.h>
+#include <TROOT.h>
+#include <TF1.h>
+#include <TMath.h>
 
-// STEER includes
-#include "AliRun.h"
-#include "AliRunLoader.h"
-#include "AliHeader.h"
-#include "AliLoader.h"
-#include "AliStack.h"
-#include "AliMagF.h"
-#include "AliESD.h"
-
-// MUON includes
-#include "AliESDMuonTrack.h"
 #endif
 
+Bool_t MUONefficiency(const char* filename = "generated/galice.root", const char* esdFileName = "AliESDs.root",
+                     const char* geoFilename = "generated/geometry.root", const char* ocdbPath = "local://$ALICE_ROOT/OCDB",
+                     Int_t ExtrapToVertex = -1, Int_t ResType = 553, Int_t FirstEvent = 0, Int_t LastEvent = 1000000 )
+{ 
+/// \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, 443 for J/Psi (default 553)  
+/// \param FirstEvent (default 0)
+/// \param LastEvent  (default 1.e6)
+/// \param Chi2Cut    to keep only tracks with chi2 per d.o.f. < Chi2Cut (default 100)  
 
-Bool_t MUONefficiency( Int_t ResType = 553, Int_t FirstEvent = 0, Int_t LastEvent = 1000000,
-                      char* esdFileName = "AliESDs.root", char* filename = "galice.root")
-{ // MUONefficiency starts
+
+  // MUONefficiency starts
 
   Double_t MUON_MASS = 0.105658369;
   Double_t UPSILON_MASS = 9.4603 ;
@@ -116,7 +129,6 @@ Bool_t MUONefficiency( Int_t ResType = 553, Int_t FirstEvent = 0, Int_t LastEven
   
   // Printing Level 
   Int_t PRINTLEVEL = 0 ;
-  Int_t SELECT =  0 ; // not used
   
   //for kinematic, i.e. reference tracks
   TNtuple *Ktuple = new TNtuple("Ktuple","Kinematics NTuple","ev:npart:id:idmo:idgdmo:p:pt:y:theta:pseudorap:vx:vy:vz");
@@ -141,7 +153,7 @@ Bool_t MUONefficiency( Int_t ResType = 553, Int_t FirstEvent = 0, Int_t LastEven
   
   hInvMassRes = new TH1F("hInvMassRes", "Mu+Mu- invariant mass (GeV/c2) around Resonance",(Int_t) (nBinsPerGev*3*countingRange*2),massResonance-3*countingRange,massResonance+3*countingRange);
 
-  
+  TH1F *hPrimaryVertex = new TH1F("hPrimaryVertex","SPD reconstructed Z vertex",150,-15,15);
   TH1F *hChi2PerDof = new TH1F("hChi2PerDof", "Muon track chi2/d.o.f.", 100, 0., 20.);
   TH1F *hNumberOfTrack = new TH1F("hNumberOfTrack","nb of track /evt ",20,-0.5,19.5);
   TH1F *hRapMuon = new TH1F("hRapMuon"," Muon Rapidity",50,-4.5,-2);
@@ -160,7 +172,12 @@ Bool_t MUONefficiency( Int_t ResType = 553, Int_t FirstEvent = 0, Int_t LastEven
   Int_t NbTrigger = 0;
   Int_t ptTrig = 0;
 
-  Double_t thetaX, thetaY, pYZ;
+  Double_t fXVertex=0;
+  Double_t fYVertex=0;
+  Double_t fZVertex=0;
+  Double_t errXVtx=0;
+  Double_t errYVtx=0;
+
   Double_t fPxRec1, fPyRec1, fPzRec1, fE1;
   Double_t fPxRec2, fPyRec2, fPzRec2, fE2;
   Int_t fCharge1, fCharge2;
@@ -171,8 +188,46 @@ Bool_t MUONefficiency( Int_t ResType = 553, Int_t FirstEvent = 0, Int_t LastEven
 
   TLorentzVector fV1, fV2, fVtot;
 
-  // set off mag field 
-  AliMagF::SetReadField(kFALSE);
+  // Import TGeo geometry (needed by AliMUONTrackExtrap::ExtrapToVertex)
+  if (!gGeoManager) {
+    TGeoManager::Import(geoFilename);
+    if (!gGeoManager) {
+      Error("MUONefficiency", "getting geometry from file %s failed", geoFilename);
+      return kFALSE;
+    }
+  }
+  
+  // open the ESD file
+  TFile* esdFile = TFile::Open(esdFileName);
+  if (!esdFile || !esdFile->IsOpen()) {
+    Error("MUONefficiency", "opening ESD file %s failed", esdFileName);
+    return kFALSE;
+  }
+
+  AliESDEvent* esd = new AliESDEvent();
+  TTree* tree = (TTree*) esdFile->Get("esdTree");
+  if (!tree) {
+    Error("MUONefficiency", "no ESD tree found");
+    return kFALSE;
+  } 
+  esd->ReadFromTree(tree);
+
+  // get run number
+  if (tree->GetEvent(0) <= 0) {
+    Error("MUONefficiency", "no ESD object found for event 0");
+    return kFALSE;
+  }
+  Int_t runNumber = esd->GetRunNumber();
+
+  // load necessary data from OCDB
+  AliCDBManager::Instance()->SetDefaultStorage(ocdbPath);
+  AliCDBManager::Instance()->SetSpecificStorage("GRP/GRP/Data",
+                         Form("local://%s",gSystem->pwd()));
+  AliCDBManager::Instance()->SetRun(runNumber);
+  if (!AliMUONCDB::LoadField()) return kFALSE;
+
+  // set the magnetic field for track extrapolations
+  AliMUONTrackExtrap::SetField();
 
   // open run loader and load gAlice, kinematics and header
   AliRunLoader* runLoader = AliRunLoader::Open(filename);
@@ -187,25 +242,16 @@ Bool_t MUONefficiency( Int_t ResType = 553, Int_t FirstEvent = 0, Int_t LastEven
     Error("MUONefficiency", "no galice object found");
     return kFALSE;
   }
-  
-  // open the ESD file
-  TFile* esdFile = TFile::Open(esdFileName);
-  if (!esdFile || !esdFile->IsOpen()) {
-    Error("MUONefficiency", "opening ESD file %s failed", esdFileName);
-    return kFALSE;
-  }
-  
-  AliESD* esd = new AliESD();
-  TTree* tree = (TTree*) esdFile->Get("esdTree");
-  if (!tree) {
-    Error("CheckESD", "no ESD tree found");
+
+  runLoader->LoadHeader();
+  if (runNumber != runLoader->GetHeader()->GetRun()) {
+    Error("MUONefficiency", "mismatch between run number from ESD and from runLoader");
     return kFALSE;
   }
-  tree->SetBranchAddress("ESD", &esd);
 
-  runLoader->LoadHeader();
   nevents = runLoader->GetNumberOfEvents();
-  
+  AliMUONTrackParam trackParam;
+
   // to access the particle  Stack
   runLoader->LoadKinematics("READ");
 
@@ -213,16 +259,8 @@ Bool_t MUONefficiency( Int_t ResType = 553, Int_t FirstEvent = 0, Int_t LastEven
 
   TParticle *particle; 
   
-  Int_t track1Id = 0 ;
-  Int_t track1PDGId = 0 ;
-  Int_t track1MotherId = 0 ;
-  Int_t track1MotherPDGId = 0 ;
   Int_t track1Trigger = 0 ;
   Float_t track1TriggerChi2 = 0 ;
-  Int_t track2Id = 0 ;
-  Int_t track2PDGId = 0 ;
-  Int_t track2MotherId = 0 ;
-  Int_t track2MotherPDGId = 0 ;
   Int_t track2Trigger = 0 ;
   Float_t track2TriggerChi2 = 0 ;
 
@@ -246,8 +284,7 @@ Bool_t MUONefficiency( Int_t ResType = 553, Int_t FirstEvent = 0, Int_t LastEven
     Int_t ntracks = theStack->GetNtrack();
   
     if (PRINTLEVEL || (iEvent%100==0)) printf("\n  >>> Event %d \n",iEvent);
-    if (PRINTLEVEL) cout << nprimarypart << " Particles generated (total is " << ntracks << ")"<< endl ;
-    
+    if (PRINTLEVEL) cout << nprimarypart << " Particles generated (total is " << ntracks << ")"<< endl ;    
     
     for(Int_t iparticle=0; iparticle<nparticles; iparticle++) { // Start loop over particles
       particle = theStack->Particle(iparticle);
@@ -291,13 +328,23 @@ Bool_t MUONefficiency( Int_t ResType = 553, Int_t FirstEvent = 0, Int_t LastEven
 
     
     // get the event summary data
-    tree->GetEvent(iEvent);
-    if (!esd) {
+    if (tree->GetEvent(iEvent) <= 0) {
       Error("CheckESD", "no ESD object found for event %d", iEvent);
       return kFALSE;
     }
-
-    Int_t triggerWord = esd->GetTrigger();
+    
+    // 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 triggerWord = esd->GetTriggerMask();
     Int_t nTracks = (Int_t)esd->GetNumberOfMuonTracks() ; 
 
     if (PRINTLEVEL > 0){
@@ -308,17 +355,20 @@ Bool_t MUONefficiency( Int_t ResType = 553, Int_t FirstEvent = 0, Int_t LastEven
     // loop over all reconstructed tracks (also first track of combination)
     for (Int_t iTrack = 0; iTrack <  nTracks;  iTrack++) {
 
-      AliESDMuonTrack* muonTrack = esd->GetMuonTrack(iTrack);
+      // skip ghosts
+      if (!esd->GetMuonTrack(iTrack)->ContainTrackerData()) continue;
 
-      //if (PRINTLEVEL > 5) cout << "1st muonTrack->GetTrackID() : " << track1Id << endl; 
+      AliESDMuonTrack* muonTrack = new AliESDMuonTrack(*(esd->GetMuonTrack(iTrack)));
 
-      if(SELECT && track1Id) {
-       particle = theStack->Particle(track1Id);
-       track1PDGId = particle->GetPdgCode() ;
-       track1MotherId = particle->GetFirstMother();
-       if (track1MotherId >=0 )
-         track1MotherPDGId = ((TParticle*) theStack->Particle(track1MotherId))->GetPdgCode();
-       if (PRINTLEVEL > 0) cout << "track1MotherPDGId = " << track1MotherPDGId << endl ;
+      // 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
       }
 
       // Trigger
@@ -329,17 +379,9 @@ Bool_t MUONefficiency( Int_t ResType = 553, Int_t FirstEvent = 0, Int_t LastEven
       else 
        track1TriggerChi2 = 0. ;
 
-      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);
       fCharge1 = Int_t(TMath::Sign(1.,muonTrack->GetInverseBendingMomentum()));
       
-      fE1 = TMath::Sqrt(MUON_MASS * MUON_MASS + fPxRec1 * fPxRec1 + fPyRec1 * fPyRec1 + fPzRec1 * fPzRec1);
-      fV1.SetPxPyPzE(fPxRec1, fPyRec1, fPzRec1, fE1);
+      muonTrack->LorentzP(fV1);
 
       ntrackhits = muonTrack->GetNHit();
       fitfmin    = muonTrack->GetChi2();
@@ -356,7 +398,7 @@ Bool_t MUONefficiency( Int_t ResType = 553, Int_t FirstEvent = 0, Int_t LastEven
       // chi2 per d.o.f.
       
       Float_t ch1 =  fitfmin / (2.0 * ntrackhits - 5);
-      if (PRINTLEVEL > 5 ) printf(" px %f py %f pz %f pt %f NHits %d  Norm.chi2 %f charge %d\n",fPxRec1, fPyRec1, fPzRec1, pt1, ntrackhits, ch1, fCharge1);
+      if (PRINTLEVEL > 5 ) printf(" px %f py %f pz %f pt %f NHits %d  Norm.chi2 %f charge %d\n",fV1.Px(), fV1.Py(), fV1.Pz(), pt1, ntrackhits, ch1, fCharge1);
       
       
       if ((ch1 < Chi2Cut) && (pt1 > PtCutMin) && (pt1 < PtCutMax)) { // condition for good track (Chi2Cut and PtCut)
@@ -370,47 +412,43 @@ Bool_t MUONefficiency( Int_t ResType = 553, Int_t FirstEvent = 0, Int_t LastEven
 
        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);
-         //Int_t track2Id = muonTrack->GetTrackID();
-         //if (PRINTLEVEL > 5) cout << "2nd muonTrack->GetTrackID() : " << track2Id << endl;
-
-         if(SELECT && track2Id) {
-           particle = theStack->Particle(track2Id);
-           track2PDGId = particle->GetPdgCode();
-           track2MotherId = particle->GetFirstMother();
-           if (track2MotherId >=0 ) 
-             track2MotherPDGId = ((TParticle*) theStack->Particle(track2MotherId))->GetPdgCode();
+         // 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
          }
 
-         track2Trigger = muonTrack->GetMatchTrigger();
+         track2Trigger = muonTrack2->GetMatchTrigger();
          if (track2Trigger) 
-           track2TriggerChi2 = muonTrack->GetChi2MatchTrigger();
+           track2TriggerChi2 = muonTrack2->GetChi2MatchTrigger();
          else 
            track2TriggerChi2 = 0. ;
 
-         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()));
+         fCharge2 = Int_t(TMath::Sign(1.,muonTrack2->GetInverseBendingMomentum()));
 
-         fE2 = TMath::Sqrt(MUON_MASS * MUON_MASS + 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();
@@ -418,64 +456,58 @@ Bool_t MUONefficiency( Int_t ResType = 553, Int_t FirstEvent = 0, Int_t LastEven
          // chi2 per d.o.f.
          Float_t ch2 = fitfmin  / (2.0 * ntrackhits - 5);
 
-         if (PRINTLEVEL > 5)  cout << "track1MotherId : "<< track1MotherId << "    track2MotherId : " << track2MotherId << endl ;
-         if (PRINTLEVEL > 5) cout << "track1MotherPDGId : " << track1MotherPDGId << "    track2MotherPDGId : "  << track2MotherPDGId << endl ;
 
-         // Select Condition
-         if (!SELECT || (track2MotherId == track1MotherId && track2MotherPDGId == ResType && TMath::Abs(track1PDGId)==13 && TMath::Abs(track2PDGId)==13 )) {
-
-           // condition for good track (Chi2Cut and PtCut)
-           if ((ch2 < Chi2Cut) && (pt2 > PtCutMin)  && (pt2 < PtCutMax)) {
+         // condition for good track (Chi2Cut and PtCut)
+         if ((ch2 < Chi2Cut) && (pt2 > PtCutMin)  && (pt2 < PtCutMax)) {
            
-             // condition for opposite charges
-             if ((fCharge1 * fCharge2) == -1) {
-
-               if (PRINTLEVEL > 8) cout << "---------> Now filling the Ntuple "  <<  endl ;
-               
-               // invariant mass
-               fVtot = fV1 + fV2;
-               Float_t invMass = fVtot.M();
-               
-               if (fCharge1 < 0){ //mu_minus is index 1 in the ntuple
-                 Float_t ESDFill[16] = {iEvent,triggerWord,fVtot.Pt(),fVtot.Rapidity(),fVtot.Theta()/TMath::Pi()*180,invMass,fV1.Pt(),fV1.Rapidity(),fV1.Theta()/TMath::Pi()*180,fCharge1,track1TriggerChi2,fV2.Pt(),fV2.Rapidity(),fV2.Theta()/TMath::Pi()*180,fCharge2,track2TriggerChi2};
-                 ESDtuple->Fill(ESDFill);
-               }
-               else{
-                 Float_t ESDFill[16] = {iEvent,triggerWord,fVtot.Pt(),fVtot.Rapidity(),fVtot.Theta()/TMath::Pi()*180,invMass,fV2.Pt(),fV2.Rapidity(),fV2.Theta()/TMath::Pi()*180,fCharge2,track2TriggerChi2,fV1.Pt(),fV1.Rapidity(),fV1.Theta()/TMath::Pi()*180,fCharge1,track1TriggerChi2};
-                 ESDtuple->Fill(ESDFill);
-               }
-               
-               // fill histos hInvMassAll and hInvMassRes
-               hInvMassAll->Fill(invMass);
-               hInvMassRes->Fill(invMass);
-               hInvMassAll_vs_Pt->Fill(invMass,fVtot.Pt());
-
-               //trigger info 
-               if (ResType == 553)
-                 ptTrig = 0x400;// mask for Hpt unlike sign pair
-               else if (ResType == 443)
-                 ptTrig = 0x800;// mask for Apt unlike sign pair
-               else 
-                 ptTrig = 0x200;// mask for Lpt unlike sign pair
+           // condition for opposite charges
+           if ((fCharge1 * fCharge2) == -1) {
+             
+             if (PRINTLEVEL > 8) cout << "---------> Now filling the Ntuple "  <<  endl ;
+             
+             // invariant mass
+             fVtot = fV1 + fV2;
+             Float_t invMass = fVtot.M();
+             
+             if (fCharge1 < 0){ //mu_minus is index 1 in the ntuple
+               Float_t ESDFill[16] = {iEvent,triggerWord,fVtot.Pt(),fVtot.Rapidity(),fVtot.Theta()/TMath::Pi()*180,invMass,fV1.Pt(),fV1.Rapidity(),fV1.Theta()/TMath::Pi()*180,fCharge1,track1TriggerChi2,fV2.Pt(),fV2.Rapidity(),fV2.Theta()/TMath::Pi()*180,fCharge2,track2TriggerChi2};
+               ESDtuple->Fill(ESDFill);
+             }
+             else{
+               Float_t ESDFill[16] = {iEvent,triggerWord,fVtot.Pt(),fVtot.Rapidity(),fVtot.Theta()/TMath::Pi()*180,invMass,fV2.Pt(),fV2.Rapidity(),fV2.Theta()/TMath::Pi()*180,fCharge2,track2TriggerChi2,fV1.Pt(),fV1.Rapidity(),fV1.Theta()/TMath::Pi()*180,fCharge1,track1TriggerChi2};
+               ESDtuple->Fill(ESDFill);
+             }
+             
+             // fill histos hInvMassAll and hInvMassRes
+             hInvMassAll->Fill(invMass);
+             hInvMassRes->Fill(invMass);
+             hInvMassAll_vs_Pt->Fill(invMass,fVtot.Pt());
+             
+             //trigger info 
+             if (ResType == 553)
+               ptTrig = 0x08;// mask for Hpt unlike sign pair
+             else if (ResType == 443)
+               ptTrig = 0x04;// mask for Lpt unlike sign pair
+             
+             
+             if (esd->GetTriggerMask() &  ptTrig) NbTrigger++;
+             
+             if (invMass > invMassMinInPeak && invMass < invMassMaxInPeak) {
+               EventInMass++;
+               hRapResonance->Fill(fVtot.Rapidity());
+               hPtResonance->Fill(fVtot.Pt());
                
-
-               if (esd->GetTrigger() &  ptTrig) NbTrigger++;
+               // match with trigger
+               if (muonTrack2->GetMatchTrigger()>=0 && (esd->GetTriggerMask() & ptTrig))  EventInMassMatch++;
                
-               if (invMass > invMassMinInPeak && invMass < invMassMaxInPeak) {
-                 EventInMass++;
-                 hRapResonance->Fill(fVtot.Rapidity());
-                 hPtResonance->Fill(fVtot.Pt());
-
-                 // match with trigger
-                 if (muonTrack->GetMatchTrigger() && (esd->GetTrigger() & ptTrig))  EventInMassMatch++;
-
-               }
-               
-             } // if (fCharge1 * fCharge2) == -1)
-           } // if ((ch2 < Chi2Cut) && (pt2 > PtCutMin) && (pt2 < PtCutMax))
-         } // if (track2MotherId == track1MotherId && track2MotherPDGId == ResType)
+             }
+             
+           } // 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) )
+      delete muonTrack;
     } // for (Int_t iTrack = 0; iTrack < nrectracks; iTrack++)
     
     hNumberOfTrack->Fill(nTracks);
@@ -537,6 +569,7 @@ Bool_t MUONefficiency( Int_t ResType = 553, Int_t FirstEvent = 0, Int_t LastEven
   ntupleFile->Close();
   
   TFile *histoFile = new TFile("MUONhistos.root", "RECREATE");  
+  hPrimaryVertex->Write();
   hPtMuon->Write();
   hPtMuonPlus->Write();
   hPtMuonMinus->Write();
@@ -566,7 +599,17 @@ Bool_t MUONefficiency( Int_t ResType = 553, Int_t FirstEvent = 0, Int_t LastEven
   cout << "Chi2Cut for muon tracks = " << Chi2Cut << endl;
   cout << "PtCutMin for muon tracks = " << PtCutMin << endl;
   cout << "PtCutMax for muon tracks = " << PtCutMax << endl;
-  cout << "Entries (unlike sign dimuons) in the mass range  ["<<invMassMinInPeak<<";"<<invMassMaxInPeak<<"] : " << EventInMass <<endl;
+  
+  cout << "Entries (unlike sign dimuons) : " << hInvMassAll->GetEntries();
+  
+  if (hInvMassAll->GetEntries() > 0) {
+    hInvMassAll->Fit("gaus","q0");
+    TF1* f1 = hInvMassAll->GetFunction("gaus");
+    cout << Form(". Rough sigma = %7.2f MeV/c2",f1->GetParameter(2)*1000.0);
+  }
+  
+  cout << endl << "Entries (unlike sign dimuons) in the mass range  ["<<invMassMinInPeak<<";"<<invMassMaxInPeak<<"] : " << EventInMass <<endl;
+  
   if (ptTrig==0x800) cout << "Unlike Pair - All Pt" ;   
   if (ptTrig==0x400) cout << "Unlike Pair - High Pt" ;   
   if (ptTrig==0x200) cout << "Unlike Pair - Low Pt" ; 
@@ -577,4 +620,3 @@ Bool_t MUONefficiency( Int_t ResType = 553, Int_t FirstEvent = 0, Int_t LastEven
 
   return kTRUE;
 }
-