correcting baryon mass subtraction for visible energy in MC
[u/mrichter/AliRoot.git] / JETAN / AliJetAODReader.cxx
index b62ee7f..21256b4 100644 (file)
 #include "AliJetAODReaderHeader.h"
 #include "AliAODEvent.h"
 #include "AliAODTrack.h"
+#include "AliAODMCParticle.h"
 #include "AliJetDummyGeo.h"
 #include "AliJetAODFillUnitArrayTracks.h"
 #include "AliJetAODFillUnitArrayEMCalDigits.h"
-#include "AliJetHadronCorrection.h"
+#include "AliJetHadronCorrectionv1.h"
 #include "AliJetUnitArray.h"
 
 ClassImp(AliJetAODReader)
 
 AliJetAODReader::AliJetAODReader():
     AliJetReader(),
-    fChain(0x0),
-    fTree(0x0),
     fAOD(0x0),
     fRef(new TRefArray),
     fDebug(0),
@@ -61,10 +60,7 @@ AliJetAODReader::AliJetAODReader():
     fGrid2(0),
     fGrid3(0),
     fGrid4(0),
-    fPtCut(0),
-    fHCorrection(0),
     fApplyElectronCorrection(kFALSE),
-    fEFlag(kFALSE),
     fApplyMIPCorrection(kTRUE),
     fApplyFractionHadronicCorrection(kFALSE),
     fFractionHadronicCorrection(0.3),
@@ -75,10 +71,8 @@ AliJetAODReader::AliJetAODReader():
     fDZ(0),
     fNeta(0),
     fNphi(0),
-    fArrayInitialised(0),
     fRefArray(0x0),
     fProcId(kFALSE)
-
 {
   // Constructor    
 }
@@ -88,7 +82,6 @@ AliJetAODReader::AliJetAODReader():
 AliJetAODReader::~AliJetAODReader()
 {
   // Destructor
-    delete fChain;
     delete fAOD;
     delete fRef;
     delete fTpcGrid;
@@ -126,9 +119,7 @@ void AliJetAODReader::OpenInputFiles()
        if (a>=naod) continue;
        
        if (strstr(name,pattern)){
-          char path[256];
-          sprintf(path,"%s/%s/aod.root",dirName,name);
-          fChain->AddFile(path);
+          fChain->AddFile(Form("%s/%s/aod.root",dirName,name));
           a++;
        }
    }
@@ -169,6 +160,98 @@ void AliJetAODReader::ConnectTree(TTree* tree, TObject* /*data*/) {
      }
 }
 
+
+Bool_t AliJetAODReader::AcceptAODMCParticle(AliAODMCParticle *mcP,Short_t flag){
+
+  //
+  // filter for charge and for charged and neutral, no detector
+  // response yet
+  // physical priamries already selected
+
+  Int_t   pdg     = TMath::Abs(mcP->GetPdgCode());
+
+  // exclude neutrinos anyway   
+  if((pdg == 12 || pdg == 14 || pdg == 16)) return kFALSE;
+
+  if(flag==AliJetAODReaderHeader::kAllMC)return kTRUE;
+  if(flag==AliJetAODReaderHeader::kChargedMC){
+    if(mcP->Charge()==0)return kFALSE;
+    return kTRUE;
+  }
+
+  return kTRUE;
+
+}
+
+Bool_t AliJetAODReader::FillMomentumArrayMC(){
+
+  // 
+  // This routine fetches the MC particles from the AOD
+  // Depending on the flag all particles except neurinos are use
+  // or only charged particles
+  //
+
+  TClonesArray *mcArray = (TClonesArray*)fAOD->FindListObject(AliAODMCParticle::StdBranchName());
+  if(!mcArray){
+    Printf("%s:%d No MC particle branch found",(char*)__FILE__,__LINE__);
+    return kFALSE;
+  }
+  
+  Int_t nMC = mcArray->GetEntriesFast();
+  
+  // get number of tracks in event (for the loop)
+  if(fDebug>0)printf("AOD MC tracks: %5d \t", nMC);
+  
+  // temporary storage of signal and pt cut flag
+  Int_t* sflag  = new Int_t[nMC];
+  Int_t* cflag  = new Int_t[nMC];
+  
+  // get cuts set by user
+  Float_t ptMin =  fReaderHeader->GetPtCut();
+  Float_t etaMin = fReaderHeader->GetFiducialEtaMin();
+  Float_t etaMax = fReaderHeader->GetFiducialEtaMax();  
+  Int_t mcTrack = 0;
+  Float_t pt, eta;
+  TVector3 p3;
+
+  Short_t readerFlag = ((AliJetAODReaderHeader*)fReaderHeader)->GetReadAODMC();
+
+
+  for (Int_t it = 0; it < nMC; it++) {
+    AliAODMCParticle *track = (AliAODMCParticle*)mcArray->At(it);
+    if(!track->IsPhysicalPrimary())continue;
+
+    p3.SetXYZ(track->Px(),track->Py(),track->Pz());
+    eta = p3.Eta();
+    if ( (eta > etaMax) || (eta < etaMin)) continue;      // checking eta cut
+    if(!AcceptAODMCParticle(track,readerFlag))continue;
+    pt = p3.Pt();
+
+
+    
+
+    if (mcTrack == 0){
+      fRef->Delete(); // make sure to delete before placement new...
+      new(fRef) TRefArray(TProcessID::GetProcessWithUID(track));
+    }
+    new ((*fMomentumArray)[mcTrack]) TLorentzVector(p3,p3.Mag());
+    sflag[mcTrack] = 1;
+    cflag[mcTrack] = ( pt > ptMin ) ? 1: 0;
+    fRef->Add(track);
+    mcTrack++;
+  }
+  if(fDebug>0)printf("Used MC tracks: %5d \n", mcTrack);
+  // set the signal flags
+  fSignalFlag.Set(mcTrack,sflag);
+  fCutFlag.Set(mcTrack,cflag);
+
+  delete [] sflag;
+  delete [] cflag;
+
+  return kTRUE;
+}
+
 //____________________________________________________________________________
 
 Bool_t AliJetAODReader::FillMomentumArray()
@@ -177,14 +260,18 @@ Bool_t AliJetAODReader::FillMomentumArray()
   ClearArray();
   fRef->Clear();
   fDebug = fReaderHeader->GetDebug();
-  
+
   if (!fAOD) {
       return kFALSE;
   }
-  
+
+  if(((AliJetAODReaderHeader*)fReaderHeader)->GetReadAODMC()){
+    return FillMomentumArrayMC();
+  }
+   
   // get number of tracks in event (for the loop)
   Int_t nt = fAOD->GetNTracks();
-  printf("AOD tracks: %5d \t", nt);
+  if(fDebug>0)printf("AOD tracks: %5d \t", nt);
   
   // temporary storage of signal and pt cut flag
   Int_t* sflag  = new Int_t[nt];
@@ -213,13 +300,17 @@ Bool_t AliJetAODReader::FillMomentumArray()
     if((filterMask>0)&&!(track->TestFilterBit(filterMask)))continue;
     if ( (eta > etaMax) || (eta < etaMin)) continue;      // checking eta cut
 
+    if (aodTrack == 0){
+      fRef->Delete(); // make sure to delete before placement new...
+      new(fRef) TRefArray(TProcessID::GetProcessWithUID(track));
+    }
     new ((*fMomentumArray)[aodTrack]) TLorentzVector(p3,p3.Mag());
     sflag[aodTrack] = (TMath::Abs(track->GetLabel()) < 10000) ? 1 : 0;
     cflag[aodTrack] = ( pt > ptMin ) ? 1: 0;
     aodTrack++;
     fRef->Add(track);
   }
-  printf("Used AOD tracks: %5d \n", aodTrack);
+  if(fDebug>0)printf("Used AOD tracks: %5d \n", aodTrack);
   // set the signal flags
   fSignalFlag.Set(aodTrack,sflag);
   fCutFlag.Set(aodTrack,cflag);
@@ -294,6 +385,10 @@ void AliJetAODReader::SetFractionHadronicCorrection(Double_t val)
 //____________________________________________________________________________
 void AliJetAODReader::CreateTasks(TChain* tree)
 {
+  //
+  // For reader task initialization
+  //
+
   fDebug = fReaderHeader->GetDebug();
   fDZ = fReaderHeader->GetDZ();
   fTree = tree;
@@ -329,7 +424,7 @@ void AliJetAODReader::CreateTasks(TChain* tree)
   fFillUAFromEMCalDigits->SetApplyFractionHadronicCorrection(fApplyFractionHadronicCorrection);
   fFillUAFromEMCalDigits->SetFractionHadronicCorrection(fFractionHadronicCorrection);
   fFillUAFromEMCalDigits->SetApplyElectronCorrection(fApplyElectronCorrection);
-
+  // Add the task to global task
   fFillUnitArray->Add(fFillUAFromTracks);
   fFillUnitArray->Add(fFillUAFromEMCalDigits);
   fFillUAFromTracks->SetActive(kFALSE);
@@ -355,6 +450,7 @@ Bool_t AliJetAODReader::ExecTasks(Bool_t procid, TRefArray* refArray)
 
   // clear momentum array
   ClearArray();
+  fRef->Clear();
 
   fDebug = fReaderHeader->GetDebug();
   fOpt = fReaderHeader->GetDetector();
@@ -369,6 +465,9 @@ Bool_t AliJetAODReader::ExecTasks(Bool_t procid, TRefArray* refArray)
     fFillUAFromTracks->SetActive(kTRUE);
     fFillUAFromTracks->SetUnitArray(fUnitArray);
     fFillUAFromTracks->SetRefArray(fRefArray);
+    fFillUAFromTracks->SetReferences(fRef);
+    fFillUAFromTracks->SetSignalFlag(fSignalFlag);
+    fFillUAFromTracks->SetCutFlag(fCutFlag);
     fFillUAFromTracks->SetProcId(fProcId);
     //    fFillUAFromTracks->ExecuteTask("tpc"); // => Temporarily changed
     fFillUAFromTracks->Exec("tpc");
@@ -376,6 +475,8 @@ Bool_t AliJetAODReader::ExecTasks(Bool_t procid, TRefArray* refArray)
       fNumCandidate = fFillUAFromTracks->GetMult();
       fNumCandidateCut = fFillUAFromTracks->GetMultCut();
     }
+    fSignalFlag = fFillUAFromTracks->GetSignalFlag();
+    fCutFlag = fFillUAFromTracks->GetCutFlag();
   }
 
   // Digits only or Digits+TPC
@@ -433,13 +534,6 @@ void AliJetAODReader::InitParameters()
 {
   // Initialise parameters
   fOpt = fReaderHeader->GetDetector();
-  //  fHCorrection    = 0;          // For hadron correction
-  fHadCorr        = 0;          // For hadron correction
-  if(fEFlag==kFALSE){
-    if(fOpt==0 || fOpt==1)
-      fECorrection    = 0;        // For electron correction
-    else fECorrection = 1;        // For electron correction
-  }
   fNumUnits       = fGeom->GetNCells();      // Number of cells in EMCAL
   if(fDebug>1) printf("\n EMCal parameters initiated ! \n");
 }
@@ -449,21 +543,21 @@ void AliJetAODReader::InitUnitArray()
 {
   //Initialises unit arrays
   Int_t nElements = fTpcGrid->GetNEntries();
-  Float_t eta = 0., phi = 0., Deta = 0., Dphi = 0.;
+  Float_t eta = 0., phi = 0., deltaEta = 0., deltaPhi = 0.;
   if(fArrayInitialised) fUnitArray->Delete();
 
   if(fTpcGrid->GetGridType()==0)
     { // Fill the following quantities :
-      // Good track ID, (Eta,Phi) position ID, eta, phi, energy, px, py, pz, Deta, Dphi,
+      // Good track ID, (Eta,Phi) position ID, eta, phi, energy, px, py, pz, deltaEta, deltaPhi,
       // detector flag, in/out jet, pt cut, mass, cluster ID)
       for(Int_t nBin = 1; nBin < nElements+1; nBin++)
         {
           //      fTpcGrid->GetEtaPhiFromIndex2(nBin,eta,phi);
           fTpcGrid->GetEtaPhiFromIndex2(nBin,phi,eta);
           phi = ((phi < 0) ? phi + 2. * TMath::Pi() : phi);
-          Deta = fTpcGrid->GetDeta();
-          Dphi = fTpcGrid->GetDphi();
-          new ((*fUnitArray)[nBin-1]) AliJetUnitArray(nBin-1,0,eta,phi,0.,Deta,Dphi,kTpc,kOutJet,kPtSmaller,kPtSmaller,kBad,0.,-1);
+          deltaEta = fTpcGrid->GetDeta();
+          deltaPhi = fTpcGrid->GetDphi();
+          new ((*fUnitArray)[nBin-1]) AliJetUnitArray(nBin-1,0,eta,phi,0.,deltaEta,deltaPhi,kTpc,kOutJet,kPtSmaller,kPtSmaller,kBad,0.,-1);
         }
     }
 
@@ -519,60 +613,60 @@ void AliJetAODReader::InitUnitArray()
              fGeom->EtaPhiFromIndex(nBin, eta, phi); // From EMCal geometry 
              // fEmcalGrid->GetEtaPhiFromIndex2(nBin,phi,eta); // My function from Grid
              phi = ((phi < 0) ? phi + 2. * TMath::Pi() : phi);
-             Deta = fEmcalGrid->GetDeta(); // Modify with the exact detector values
-             Dphi = fEmcalGrid->GetDphi(); // Modify with the exact detector values
-             new ((*fUnitArray)[nBin]) AliJetUnitArray(nBin,0,eta,phi,0.,Deta,Dphi,kTpc,kOutJet,kPtSmaller,kPtSmaller,kBad,0.,-1);
+             deltaEta = fEmcalGrid->GetDeta(); // Modify with the exact detector values
+             deltaPhi = fEmcalGrid->GetDphi(); // Modify with the exact detector values
+             new ((*fUnitArray)[nBin]) AliJetUnitArray(nBin,0,eta,phi,0.,deltaEta,deltaPhi,kTpc,kOutJet,kPtSmaller,kPtSmaller,kBad,0.,-1);
            } 
          else {
            if(nBin>=fNumUnits && nBin<fNumUnits+nElements){
              fTpcGrid->GetEtaPhiFromIndex2(nBin+1-fNumUnits,phi,eta);
              phi = ((phi < 0) ? phi + 2. * TMath::Pi() : phi);
-             Deta = fTpcGrid->GetDeta();
-             Dphi = fTpcGrid->GetDphi();
-             new ((*fUnitArray)[nBin]) AliJetUnitArray(nBin,0,eta,phi,0.,Deta,Dphi,kTpc,kOutJet,kPtSmaller,kPtSmaller,kBad,0.,-1);
+             deltaEta = fTpcGrid->GetDeta();
+             deltaPhi = fTpcGrid->GetDphi();
+             new ((*fUnitArray)[nBin]) AliJetUnitArray(nBin,0,eta,phi,0.,deltaEta,deltaPhi,kTpc,kOutJet,kPtSmaller,kPtSmaller,kBad,0.,-1);
            }
            else {
              if(fDZ) {
                if(nBin>=fNumUnits+nElements && nBin<fNumUnits+nElements+nGaps){
                  if(nBin<fNumUnits+nElements+n0)
                    {
-                     Float_t phi = eta = 0.;
+                     phi = eta = 0.;
                      fGrid0->GetEtaPhiFromIndex2(nBin+1-(fNumUnits+nElements),phi,eta);
-                     Deta = fGrid0->GetDeta(); 
-                     Dphi = fGrid0->GetDphi(); 
-                     new ((*fUnitArray)[nBin]) AliJetUnitArray(nBin,0,eta,phi,0.,Deta,Dphi,kTpc,kOutJet,kPtSmaller,kPtSmaller,kBad,0.,-1);
+                     deltaEta = fGrid0->GetDeta(); 
+                     deltaPhi = fGrid0->GetDphi(); 
+                     new ((*fUnitArray)[nBin]) AliJetUnitArray(nBin,0,eta,phi,0.,deltaEta,deltaPhi,kTpc,kOutJet,kPtSmaller,kPtSmaller,kBad,0.,-1);
                    }
                  else if(nBin>=fNumUnits+nElements+n0 && nBin<fNumUnits+nElements+n0+n1)
                    {
-                     Float_t phi = eta = 0.;
+                     phi = eta = 0.;
                      fGrid1->GetEtaPhiFromIndex2(nBin+1-(fNumUnits+nElements+n0),phi,eta);
-                     Deta = fGrid1->GetDeta(); 
-                     Dphi = fGrid1->GetDphi(); 
-                     new ((*fUnitArray)[nBin]) AliJetUnitArray(nBin,0,eta,phi,0.,Deta,Dphi,kTpc,kOutJet,kPtSmaller,kPtSmaller,kBad,0.,-1);
+                     deltaEta = fGrid1->GetDeta(); 
+                     deltaPhi = fGrid1->GetDphi(); 
+                     new ((*fUnitArray)[nBin]) AliJetUnitArray(nBin,0,eta,phi,0.,deltaEta,deltaPhi,kTpc,kOutJet,kPtSmaller,kPtSmaller,kBad,0.,-1);
                    }
                  else if(nBin>=fNumUnits+nElements+n0+n1 && nBin<fNumUnits+nElements+n0+n1+n2)
                    {
-                     Float_t phi = eta = 0.;
+                     phi = eta = 0.;
                      fGrid2->GetEtaPhiFromIndex2(nBin+1-(fNumUnits+nElements+n0+n1),phi,eta);
-                     Deta = fGrid2->GetDeta(); 
-                     Dphi = fGrid2->GetDphi(); 
-                     new ((*fUnitArray)[nBin]) AliJetUnitArray(nBin,0,eta,phi,0.,Deta,Dphi,kTpc,kOutJet,kPtSmaller,kPtSmaller,kBad,0.,-1);
+                     deltaEta = fGrid2->GetDeta(); 
+                     deltaPhi = fGrid2->GetDphi(); 
+                     new ((*fUnitArray)[nBin]) AliJetUnitArray(nBin,0,eta,phi,0.,deltaEta,deltaPhi,kTpc,kOutJet,kPtSmaller,kPtSmaller,kBad,0.,-1);
                    }
                  else if(nBin>=fNumUnits+nElements+n0+n1+n2 && nBin<fNumUnits+nElements+n0+n1+n2+n3)
                    {
-                     Float_t phi = eta = 0.;
+                     phi = eta = 0.;
                      fGrid3->GetEtaPhiFromIndex2(nBin+1-(fNumUnits+nElements+n0+n1+n2),phi,eta);
-                     Deta = fGrid3->GetDeta(); 
-                     Dphi = fGrid3->GetDphi(); 
-                     new ((*fUnitArray)[nBin]) AliJetUnitArray(nBin,0,eta,phi,0.,Deta,Dphi,kTpc,kOutJet,kPtSmaller,kPtSmaller,kBad,0.,-1);
+                     deltaEta = fGrid3->GetDeta(); 
+                     deltaPhi = fGrid3->GetDphi(); 
+                     new ((*fUnitArray)[nBin]) AliJetUnitArray(nBin,0,eta,phi,0.,deltaEta,deltaPhi,kTpc,kOutJet,kPtSmaller,kPtSmaller,kBad,0.,-1);
                    }
                  else if(nBin>=fNumUnits+nElements+n0+n1+n2+n3 && nBin<fNumUnits+nElements+nGaps)
                    {
-                     Float_t phi = eta = 0.;
+                     phi = eta = 0.;
                      fGrid4->GetEtaPhiFromIndex2(nBin+1-(fNumUnits+nElements+n0+n1+n2+n3),phi,eta);
-                     Deta = fGrid4->GetDeta(); 
-                     Dphi = fGrid4->GetDphi(); 
-                     new ((*fUnitArray)[nBin]) AliJetUnitArray(nBin,0,eta,phi,0.,Deta,Dphi,kTpc,kOutJet,kPtSmaller,kPtSmaller,kBad,0.,-1);
+                     deltaEta = fGrid4->GetDeta(); 
+                     deltaPhi = fGrid4->GetDphi(); 
+                     new ((*fUnitArray)[nBin]) AliJetUnitArray(nBin,0,eta,phi,0.,deltaEta,deltaPhi,kTpc,kOutJet,kPtSmaller,kPtSmaller,kBad,0.,-1);
                    }
                }
              } // end if(fDZ)
@@ -582,3 +676,4 @@ void AliJetAODReader::InitUnitArray()
     } // end grid type == 1
   fArrayInitialised = 1;
 }
+