]> git.uio.no Git - u/mrichter/AliRoot.git/blobdiff - PWGHF/vertexingHF/AliAnalysisTaskSEDvsMultiplicity.cxx
Fix for D0 filtering (A.Dainese)
[u/mrichter/AliRoot.git] / PWGHF / vertexingHF / AliAnalysisTaskSEDvsMultiplicity.cxx
index e4ffbc7b400420e22bddc602520e600fa39f0f22..521ed239633b8020cebc35ccaadef4a02d6bbc16 100644 (file)
@@ -30,7 +30,7 @@
 #include <TH2F.h>     
 #include <TH3F.h>
 #include <THnSparse.h>
-#include <TProfile.h>  
+#include <TProfile.h>
 #include "AliAnalysisManager.h"
 #include "AliRDHFCuts.h"
 #include "AliRDHFCutsDplustoKpipi.h"
@@ -62,6 +62,15 @@ AliAnalysisTaskSE(),
   fHistNtrCorrEta1vsNtrRawEta1(0),
   fHistNtrVsZvtx(0),
   fHistNtrCorrVsZvtx(0),
+  fHistNtrVsNchMC(0),
+  fHistNtrCorrVsNchMC(0),
+  fHistNtrVsNchMCPrimary(0),
+  fHistNtrCorrVsNchMCPrimary(0),
+  fHistNtrVsNchMCPhysicalPrimary(0),
+  fHistNtrCorrVsNchMCPhysicalPrimary(0),
+  fHistGenPrimaryParticlesInelGt0(0),
+  fHistNchMCVsNchMCPrimaryVsNchMCPhysicalPrimary(0),
+  fHistNtrUnCorrEvSel(0),
   fHistNtrCorrEvSel(0),
   fHistNtrCorrEvWithCand(0),
   fHistNtrCorrEvWithD(0),
@@ -83,7 +92,8 @@ AliAnalysisTaskSE(),
   fReadMC(kFALSE),
   fMCOption(0),
   fUseBit(kTRUE),
-  fRefMult(9.5),
+  fSubtractTrackletsFromDau(kFALSE),
+  fRefMult(9.26),
   fPdgMeson(411)
 {
    // Default constructor
@@ -103,6 +113,15 @@ AliAnalysisTaskSEDvsMultiplicity::AliAnalysisTaskSEDvsMultiplicity(const char *n
   fHistNtrCorrEta1vsNtrRawEta1(0),
   fHistNtrVsZvtx(0),
   fHistNtrCorrVsZvtx(0),
+  fHistNtrVsNchMC(0),
+  fHistNtrCorrVsNchMC(0),
+  fHistNtrVsNchMCPrimary(0),
+  fHistNtrCorrVsNchMCPrimary(0),
+  fHistNtrVsNchMCPhysicalPrimary(0),
+  fHistNtrCorrVsNchMCPhysicalPrimary(0),
+  fHistGenPrimaryParticlesInelGt0(0),
+  fHistNchMCVsNchMCPrimaryVsNchMCPhysicalPrimary(0),
+  fHistNtrUnCorrEvSel(0),
   fHistNtrCorrEvSel(0),
   fHistNtrCorrEvWithCand(0),
   fHistNtrCorrEvWithD(0),
@@ -124,7 +143,8 @@ AliAnalysisTaskSEDvsMultiplicity::AliAnalysisTaskSEDvsMultiplicity(const char *n
   fReadMC(kFALSE),
   fMCOption(0),
   fUseBit(kTRUE),
-  fRefMult(9.5),
+  fSubtractTrackletsFromDau(kFALSE),
+  fRefMult(9.26),
   fPdgMeson(pdgMeson)
 {
   // 
@@ -235,19 +255,34 @@ void AliAnalysisTaskSEDvsMultiplicity::UserCreateOutputObjects()
   fOutput->SetOwner();
   fOutput->SetName("OutputHistos");
 
-  fHistNtrCorrEvSel = new TH1F("hNtrCorrEvSel","Tracklets multiplicity for selected events; Tracklets ; Entries",200,0.,200.);
-  fHistNtrCorrEvWithCand = new TH1F("hNtrCorrEvWithCand", "Tracklets multiplicity for events with D candidates; Tracklets ; Entries",200,0.,200.);// Total multiplicity
-  fHistNtrCorrEvWithD = new TH1F("hNtrCorrEvWithD", "Tracklets multiplicity for events with D in mass region ; Tracklets ; Entries",200,0.,200.); // 
+  fHistNtrUnCorrEvSel = new TH1F("hNtrUnCorrEvSel","Uncorrected tracklets multiplicity for selected events; Tracklets ; Entries",200,-0.5,199.5);
+  fHistNtrCorrEvSel = new TH1F("hNtrCorrEvSel","Corrected tracklets multiplicity for selected events; Tracklets ; Entries",200,-0.5,199.5);
+  fHistNtrCorrEvWithCand = new TH1F("hNtrCorrEvWithCand", "Tracklets multiplicity for events with D candidates; Tracklets ; Entries",200,-0.5,199.5);// Total multiplicity
+  fHistNtrCorrEvWithD = new TH1F("hNtrCorrEvWithD", "Tracklets multiplicity for events with D in mass region ; Tracklets ; Entries",200,-0.5,199.5); // 
   fHistNtrEta16vsNtrEta1 = new TH2F("hNtrEta16vsNtrEta1","Uncorrected Eta1.6 vs Eta1.0; Ntracklets #eta<1.0; Ntracklets #eta<1.6",200,-0.5,199.5,200,-0.5,199.5); //eta 1.6 vs eta 1.0 histogram 
-  fHistNtrCorrEta1vsNtrRawEta1 = new TH2F("hNtrCorrEta1vsNtrRawEta1","Corrected Eta1 vs Eta1.0; Ntracklets #eta<1.0 corrected; Ntracklets #eta<1",200,-0.,200.,200,-0.5,199.5); //eta 1.6 vs eta 1.0 histogram 
-  fHistNtrVsZvtx = new TH2F("hNtrVsZvtx","Ntracklet vs VtxZ; VtxZ;N_{tracklet};",300,-15,15,200,0,200.); // 
-  fHistNtrCorrVsZvtx = new TH2F("hNtrCorrVsZvtx","Ntracklet vs VtxZ; VtxZ;N_{tracklet};",300,-15,15,200,0,200.); // 
+  fHistNtrCorrEta1vsNtrRawEta1 = new TH2F("hNtrCorrEta1vsNtrRawEta1","Corrected Eta1 vs Eta1.0; Ntracklets #eta<1.0 corrected; Ntracklets #eta<1",200,-0.5,199.5,200,-0.5,199.5); //eta 1.6 vs eta 1.0 histogram 
+  fHistNtrVsZvtx = new TH2F("hNtrVsZvtx","Ntracklet vs VtxZ; VtxZ;N_{tracklet};",300,-15,15,200,-0.5,199.5); // 
+  fHistNtrCorrVsZvtx = new TH2F("hNtrCorrVsZvtx","Ntracklet vs VtxZ; VtxZ;N_{tracklet};",300,-15,15,200,-0.5,199.5); // 
+
+  fHistNtrVsNchMC = new TH2F("hNtrVsNchMC","Ntracklet vs NchMC; Nch;N_{tracklet};",200,-0.5,199.5,200,-0.5,199.5); // 
+  fHistNtrCorrVsNchMC = new TH2F("hNtrCorrVsNchMC","Ntracklet vs Nch; Nch;N_{tracklet};",200,-0.5,199.5,200,-0.5,199.5); // 
+  
+  fHistNtrVsNchMCPrimary = new TH2F("hNtrVsNchMCPrimary","Ntracklet vs Nch (Primary); Nch (Primary);N_{tracklet};",200,-0.5,199.5,200,-0.5,199.5); // 
+  fHistNtrCorrVsNchMCPrimary = new TH2F("hNtrCorrVsNchMCPrimary","Ntracklet vs Nch (Primary); Nch(Primary) ;N_{tracklet};",200,-0.5,199.5,200,-0.5,199.5); // 
+
+  fHistNtrVsNchMCPhysicalPrimary = new TH2F("hNtrVsNchMCPhysicalPrimary","Ntracklet vs Nch (Physical Primary); Nch (Physical Primary);N_{tracklet};",200,-0.5,199.5,200,-0.5,199.5); // 
+  fHistNtrCorrVsNchMCPhysicalPrimary = new TH2F("hNtrCorrVsMCPhysicalPrimary","Ntracklet vs Nch (Physical Primary); Nch (Physical Primary);N_{tracklet};",200,-0.5,199.5,200,-0.5,199.5); // 
   
+  fHistGenPrimaryParticlesInelGt0 = new TH1F("hGenPrimaryParticlesInelGt0","Multiplcity of generated charged particles ; Nparticles ; Entries",200,-0.5,199.5);
 
+  fHistNchMCVsNchMCPrimaryVsNchMCPhysicalPrimary = new TH3F("fHistNchMCVsNchMCPrimaryVsNchMCPhysicalPrimary", "MC: Nch (Physical Primary) vs Nch (Primary) vs Nch (Generated); Nch (Generated); Nch (Primary); Nch (Physical Primary)",200,-0.5,199.5,200,-0.5,199.5,200,-0.5,199.5);
+
+  fHistNtrUnCorrEvSel->Sumw2();
   fHistNtrCorrEvSel->Sumw2();
   fHistNtrCorrEvWithCand->Sumw2();
   fHistNtrCorrEvWithD->Sumw2();
-
+  fHistGenPrimaryParticlesInelGt0->Sumw2();
+  fOutput->Add(fHistNtrUnCorrEvSel);
   fOutput->Add(fHistNtrCorrEvSel);
   fOutput->Add(fHistNtrCorrEvWithCand);
   fOutput->Add(fHistNtrCorrEvWithD);
@@ -256,6 +291,15 @@ void AliAnalysisTaskSEDvsMultiplicity::UserCreateOutputObjects()
   fOutput->Add(fHistNtrVsZvtx);
   fOutput->Add(fHistNtrCorrVsZvtx);
 
+  fOutput->Add(fHistNtrVsNchMC);
+  fOutput->Add(fHistNtrCorrVsNchMC);
+  fOutput->Add(fHistNtrVsNchMCPrimary);
+  fOutput->Add(fHistNtrCorrVsNchMCPrimary);
+  fOutput->Add(fHistNtrVsNchMCPhysicalPrimary);
+  fOutput->Add(fHistNtrCorrVsNchMCPhysicalPrimary);
+  fOutput->Add(fHistGenPrimaryParticlesInelGt0);
+  fOutput->Add(fHistNchMCVsNchMCPrimaryVsNchMCPhysicalPrimary);
+
   
   fHistNEvents = new TH1F("fHistNEvents", "number of events ",11,-0.5,10.5);
   fHistNEvents->GetXaxis()->SetBinLabel(1,"nEvents total");
@@ -274,15 +318,15 @@ void AliAnalysisTaskSEDvsMultiplicity::UserCreateOutputObjects()
   fHistNEvents->SetMinimum(0);
   fOutput->Add(fHistNEvents);
 
-  fPtVsMassVsMult=new TH3F("hPtVsMassvsMult", "D candidates: p_{t} vs mass vs tracklets multiplicity; Tracklets; Mass M [GeV/c^{2}]; p_{t} [GeV/c]",200,0.,200.,fNMassBins,fLowmasslimit,fUpmasslimit,48,0.,24.);
+  fPtVsMassVsMult=new TH3F("hPtVsMassvsMult", "D candidates: p_{t} vs mass vs tracklets multiplicity; Tracklets; Mass M [GeV/c^{2}]; p_{t} [GeV/c]",200,-0.5,199.5,fNMassBins,fLowmasslimit,fUpmasslimit,48,0.,24.);
  
-  fPtVsMassVsMultNoPid=new TH3F("hPtVsMassvsMultNoPid", "D candidates: p_{t} vs mass vs tracklets multiplicity; Tracklets; Mass M [GeV/c^{2}]; p_{t} [GeV/c]",200,0.,200.,fNMassBins,fLowmasslimit,fUpmasslimit,48,0.,24.); 
+  fPtVsMassVsMultNoPid=new TH3F("hPtVsMassvsMultNoPid", "D candidates: p_{t} vs mass vs tracklets multiplicity; Tracklets; Mass M [GeV/c^{2}]; p_{t} [GeV/c]",200,-0.5,199.5,fNMassBins,fLowmasslimit,fUpmasslimit,48,0.,24.); 
 
   fPtVsMassVsMultUncorr=new TH3F("hPtVsMassvsMultUncorr", "D candidates: p_{t} vs mass vs tracklets multiplicity; Tracklets; Mass M [GeV/c^{2}]; p_{t} [GeV/c]",200,-0.5,199.5,fNMassBins,fLowmasslimit,fUpmasslimit,48,0.,24.);
 
-  fPtVsMassVsMultPart=new TH3F("hPtVsMassvsMultPart", "D candidates: p_{t} vs mass vs tracklets multiplicity; Tracklets; Mass M [GeV/c^{2}]; p_{t} [GeV/c]",200,0.,200.,fNMassBins,fLowmasslimit,fUpmasslimit,48,0.,24.);
+  fPtVsMassVsMultPart=new TH3F("hPtVsMassvsMultPart", "D candidates: p_{t} vs mass vs tracklets multiplicity; Tracklets; Mass M [GeV/c^{2}]; p_{t} [GeV/c]",200,-0.5,199.5,fNMassBins,fLowmasslimit,fUpmasslimit,48,0.,24.);
 
-  fPtVsMassVsMultAntiPart=new TH3F("hPtVsMassvsMultAntiPart", "D candidates: p_{t} vs mass vs tracklets multiplicity; Tracklets; Mass M [GeV/c^{2}]; p_{t} [GeV/c]",200,0.,200.,fNMassBins,fLowmasslimit,fUpmasslimit,48,0.,24.);
+  fPtVsMassVsMultAntiPart=new TH3F("hPtVsMassvsMultAntiPart", "D candidates: p_{t} vs mass vs tracklets multiplicity; Tracklets; Mass M [GeV/c^{2}]; p_{t} [GeV/c]",200,-0.5,199.5,fNMassBins,fLowmasslimit,fUpmasslimit,48,0.,24.);
 
   fOutput->Add(fPtVsMassVsMult);
   fOutput->Add(fPtVsMassVsMultUncorr);
@@ -433,6 +477,24 @@ void AliAnalysisTaskSEDvsMultiplicity::UserExec(Option_t */*option*/)
       printf("AliAnalysisTaskSEDvsMultiplicity::UserExec: MC header branch not found!\n");
       return;
      }
+  
+
+    Int_t nChargedMC=AliVertexingHFUtils::GetGeneratedMultiplicityInEtaRange(arrayMC,-1.0,1.0);
+    Int_t nChargedMCPrimary=AliVertexingHFUtils::GetGeneratedPrimariesInEtaRange(arrayMC,-1.0,1.0);
+    Int_t nChargedMCPhysicalPrimary=AliVertexingHFUtils::GetGeneratedPhysicalPrimariesInEtaRange(arrayMC,-1.0,1.0);
+    if(nChargedMCPhysicalPrimary>0){ // INEL>0 for |eta|<1
+      fHistGenPrimaryParticlesInelGt0->Fill(nChargedMCPhysicalPrimary);
+    }
+    fHistNtrVsNchMC->Fill(nChargedMC,countTreta1);
+    fHistNtrCorrVsNchMC->Fill(nChargedMC,countTreta1corr);
+
+    fHistNtrVsNchMCPrimary->Fill(nChargedMCPrimary,countTreta1);
+    fHistNtrCorrVsNchMCPrimary->Fill(nChargedMCPrimary,countTreta1corr);
+
+    fHistNtrVsNchMCPhysicalPrimary->Fill(nChargedMCPhysicalPrimary,countTreta1);
+    fHistNtrCorrVsNchMCPhysicalPrimary->Fill(nChargedMCPhysicalPrimary,countTreta1corr);
+
+    fHistNchMCVsNchMCPrimaryVsNchMCPhysicalPrimary->Fill(nChargedMC,nChargedMCPrimary,nChargedMCPhysicalPrimary);
   }
   
   Int_t nCand = arrayCand->GetEntriesFast(); 
@@ -453,7 +515,7 @@ void AliAnalysisTaskSEDvsMultiplicity::UserExec(Option_t */*option*/)
     Double_t rapid=d->Y(fPdgMeson);
     Bool_t isFidAcc=fRDCutsAnalysis->IsInFiducialAcceptance(ptCand,rapid);
     if(!isFidAcc) continue;
-    Int_t passAllCuts=fRDCutsAnalysis->IsSelected(d,AliRDHFCuts::kCandidate,aod);
+    Int_t passAllCuts=fRDCutsAnalysis->IsSelected(d,AliRDHFCuts::kAll,aod);
     Int_t passTopolCuts=fRDCutsAnalysis->GetIsSelectedCuts();
     if(passTopolCuts==0) continue;
     nSelectedNoPID++;
@@ -462,6 +524,16 @@ void AliAnalysisTaskSEDvsMultiplicity::UserExec(Option_t */*option*/)
       nSelectedPID++;
       fHistNEvents->Fill(10);
     }
+    Double_t multForCand=countTreta1corr;
+    if(fSubtractTrackletsFromDau){
+      for(Int_t iDau=0; iDau<nDau; iDau++){
+       AliAODTrack *t = (AliAODTrack*)d->GetDaughter(iDau);
+       if(!t) continue;
+       if(t->HasPointOnITSLayer(0) && t->HasPointOnITSLayer(1)){
+         if(multForCand>0) multForCand-=1;
+       }
+      }
+    }
     Bool_t isPrimary=kTRUE;
     Int_t labD=-1;
     Double_t trueImpParXY=9999.;
@@ -488,7 +560,7 @@ void AliAnalysisTaskSEDvsMultiplicity::UserExec(Option_t */*option*/)
     for(Int_t iHyp=0; iHyp<2; iHyp++){
       if(mass[iHyp]<0.) continue; // for D+ and D* we have 1 mass hypothesis
       Double_t invMass=mass[iHyp];
-      Double_t arrayForSparse[5]={invMass,ptCand,impparXY,dlen,countTreta1corr};
+      Double_t arrayForSparse[5]={invMass,ptCand,impparXY,dlen,multForCand};
 
       if(fReadMC){
        
@@ -508,7 +580,7 @@ void AliAnalysisTaskSEDvsMultiplicity::UserExec(Option_t */*option*/)
            }else if(fPdgMeson==413){
              trueImpParXY=0.; /// FIXME
            }
-           Double_t arrayForSparseTrue[5]={invMass,ptCand,trueImpParXY,dlen,countTreta1corr};
+           Double_t arrayForSparseTrue[5]={invMass,ptCand,trueImpParXY,dlen,multForCand};
            if(fillHisto && passAllCuts){
              fHistMassPtImpPar[2]->Fill(arrayForSparse);
              fHistMassPtImpPar[3]->Fill(arrayForSparseTrue);
@@ -530,17 +602,22 @@ void AliAnalysisTaskSEDvsMultiplicity::UserExec(Option_t */*option*/)
        if(iHyp==1 && !(passTopolCuts&2)) continue; // candidate not passing as D0bar
       }
 
-      fPtVsMassVsMultNoPid->Fill(countTreta1corr,invMass,ptCand);
+      fPtVsMassVsMultNoPid->Fill(multForCand,invMass,ptCand);
+
+      if(fPdgMeson==421){
+       if(iHyp==0 && !(passAllCuts&1)) continue; // candidate not passing as D0
+       if(iHyp==1 && !(passAllCuts&2)) continue; // candidate not passing as D0bar
+      }
       if(passAllCuts){
-       fPtVsMassVsMult->Fill(countTreta1corr,invMass,ptCand);                
+       fPtVsMassVsMult->Fill(multForCand,invMass,ptCand);                    
        fPtVsMassVsMultUncorr->Fill(countTreta1,invMass,ptCand);
        // Add separation between part antipart
        if(fPdgMeson==411){
-         if(d->GetCharge()>0) fPtVsMassVsMultPart->Fill(countTreta1corr,invMass,ptCand);
-         else fPtVsMassVsMultAntiPart->Fill(countTreta1corr,invMass,ptCand);
+         if(d->GetCharge()>0) fPtVsMassVsMultPart->Fill(multForCand,invMass,ptCand);
+         else fPtVsMassVsMultAntiPart->Fill(multForCand,invMass,ptCand);
        }else if(fPdgMeson==421){
-         if(passTopolCuts&1) fPtVsMassVsMultPart->Fill(countTreta1corr,invMass,ptCand);
-         if(passTopolCuts&2) fPtVsMassVsMultAntiPart->Fill(countTreta1corr,invMass,ptCand);
+         if(passAllCuts&1) fPtVsMassVsMultPart->Fill(multForCand,invMass,ptCand);
+         if(passAllCuts&2) fPtVsMassVsMultAntiPart->Fill(multForCand,invMass,ptCand);
        }else if(fPdgMeson==413){
          // FIXME ADD Dstar!!!!!!!!
        }
@@ -555,6 +632,7 @@ void AliAnalysisTaskSEDvsMultiplicity::UserExec(Option_t */*option*/)
   }
   fCounter->StoreCandidates(aod,nSelectedNoPID,kTRUE);
   fCounter->StoreCandidates(aod,nSelectedPID,kFALSE);
+  fHistNtrUnCorrEvSel->Fill(countTreta1);
   fHistNtrCorrEvSel->Fill(countTreta1corr);
   if(nSelectedPID>0) fHistNtrCorrEvWithCand->Fill(countTreta1corr);
   if(nSelectedInMassPeak>0) fHistNtrCorrEvWithD->Fill(countTreta1corr);
@@ -669,4 +747,3 @@ TProfile* AliAnalysisTaskSEDvsMultiplicity::GetEstimatorHistogram(const AliVEven
 
   return fMultEstimatorAvg[period];
 }
-