]> git.uio.no Git - u/mrichter/AliRoot.git/blobdiff - PWGJE/UserTasks/AliAnalysisTaskPID.cxx
- Added commented code which can be used to scale feed-down from sigma+ decays
[u/mrichter/AliRoot.git] / PWGJE / UserTasks / AliAnalysisTaskPID.cxx
index 9a7941f7ce0b8bc55ffd3f19ae5c74d4297e77f9..d1ab0e851939696701836366f4adcdd19cf341c7 100644 (file)
@@ -53,8 +53,9 @@ AliAnalysisTaskPID::AliAnalysisTaskPID()
   , fInputFromOtherTask(kFALSE)
   , fDoPID(kTRUE)
   , fDoEfficiency(kTRUE)
-  , fDoPtResolution(kTRUE)
-  , fDoDeDxCheck(kTRUE)
+  , fDoPtResolution(kFALSE)
+  , fDoDeDxCheck(kFALSE)
+  , fDoBinZeroStudy(kFALSE)
   , fStoreCentralityPercentile(kFALSE)
   , fStoreAdditionalJetInformation(kFALSE)
   , fTakeIntoAccountMuons(kFALSE)
@@ -139,6 +140,10 @@ AliAnalysisTaskPID::AliAnalysisTaskPID()
   , fhEventsTriggerSel(0x0)
   , fhEventsTriggerSelVtxCut(0x0) 
   , fhEventsProcessedNoPileUpRejection(0x0)
+  , fChargedGenPrimariesTriggerSel(0x0)
+  , fChargedGenPrimariesTriggerSelVtxCut(0x0)
+  , fChargedGenPrimariesTriggerSelVtxCutZ(0x0)
+  , fChargedGenPrimariesTriggerSelVtxCutZPileUpRej(0x0)
   , fhMCgeneratedYieldsPrimaries(0x0)
   , fh2FFJetPtRec(0x0)
   , fh2FFJetPtGen(0x0)
@@ -189,8 +194,9 @@ AliAnalysisTaskPID::AliAnalysisTaskPID(const char *name)
   , fInputFromOtherTask(kFALSE)
   , fDoPID(kTRUE)
   , fDoEfficiency(kTRUE)
-  , fDoPtResolution(kTRUE)
-  , fDoDeDxCheck(kTRUE)
+  , fDoPtResolution(kFALSE)
+  , fDoDeDxCheck(kFALSE)
+  , fDoBinZeroStudy(kFALSE)
   , fStoreCentralityPercentile(kFALSE)
   , fStoreAdditionalJetInformation(kFALSE)
   , fTakeIntoAccountMuons(kFALSE)
@@ -275,6 +281,10 @@ AliAnalysisTaskPID::AliAnalysisTaskPID(const char *name)
   , fhEventsTriggerSel(0x0)
   , fhEventsTriggerSelVtxCut(0x0) 
   , fhEventsProcessedNoPileUpRejection(0x0)
+  , fChargedGenPrimariesTriggerSel(0x0)
+  , fChargedGenPrimariesTriggerSelVtxCut(0x0)
+  , fChargedGenPrimariesTriggerSelVtxCutZ(0x0)
+  , fChargedGenPrimariesTriggerSelVtxCutZPileUpRej(0x0)
   , fhMCgeneratedYieldsPrimaries(0x0)
   , fh2FFJetPtRec(0x0)
   , fh2FFJetPtGen(0x0)
@@ -594,11 +604,19 @@ void AliAnalysisTaskPID::UserCreateOutputObjects()
   // This centrality estimator deals with integers! This implies that the ranges are always [lowlim, uplim - 1]
   Double_t binsCentITSTPCTracklets[nCentBins+1] = { 0, 7, 13, 20, 29, 40, 50, 60, 72, 83, 95, 105, 115 };
   
+  // Special centrality binning for pp
+  Double_t binsCentpp[nCentBins+1] =   { 0, 0.01, 0.1, 1, 5, 10, 15, 20, 30, 40, 50, 70, 100};
+  
   if (fCentralityEstimator.CompareTo("ITSTPCtracklets", TString::kIgnoreCase) == 0 && fStoreCentralityPercentile) {
     // Special binning for this centrality estimator; but keep number of bins!
     for (Int_t i = 0; i < nCentBins+1; i++)
       binsCent[i] = binsCentITSTPCTracklets[i];
   }
+  else if (fCentralityEstimator.Contains("ppMult", TString::kIgnoreCase) && fStoreCentralityPercentile) {
+    // Special binning for this pp centrality estimator; but keep number of bins!
+    for (Int_t i = 0; i < nCentBins+1; i++)
+      binsCent[i] = binsCentpp[i];
+  }
 
   const Int_t nJetPtBins = 11;
   Double_t binsJetPt[nJetPtBins+1] = {0, 2, 5, 10, 15, 20, 30, 40, 60, 80, 120, 200};
@@ -753,26 +771,26 @@ void AliAnalysisTaskPID::UserCreateOutputObjects()
   
   
   fhEventsProcessed = new TH1D("fhEventsProcessed",
-                               "Number of events passing trigger selection, vtx and zvtx cuts and pile-up rejection;Centrality percentile", 
+                               "Number of events passing trigger selection, vtx and zvtx cuts and pile-up rejection;Centrality Percentile", 
                                nCentBins, binsCent);
   fhEventsProcessed->Sumw2();
   fOutputContainer->Add(fhEventsProcessed);
   
   fhEventsTriggerSelVtxCut = new TH1D("fhEventsTriggerSelVtxCut",
-                                      "Number of events passing trigger selection and vtx cut;Centrality percentile", 
+                                      "Number of events passing trigger selection and vtx cut;Centrality Percentile", 
                                       nCentBins, binsCent);
   fhEventsTriggerSelVtxCut->Sumw2();
   fOutputContainer->Add(fhEventsTriggerSelVtxCut);
   
   fhEventsTriggerSel = new TH1D("fhEventsTriggerSel",
-                                "Number of events passing trigger selection;Centrality percentile", 
+                                "Number of events passing trigger selection;Centrality Percentile", 
                                 nCentBins, binsCent);
   fOutputContainer->Add(fhEventsTriggerSel);
   fhEventsTriggerSel->Sumw2();
   
   
   fhEventsProcessedNoPileUpRejection = new TH1D("fhEventsProcessedNoPileUpRejection",
-                                                "Number of events passing trigger selection, vtx and zvtx cuts;Centrality percentile", 
+                                                "Number of events passing trigger selection, vtx and zvtx cuts;Centrality Percentile", 
                                                 nCentBins, binsCent);
   fOutputContainer->Add(fhEventsProcessedNoPileUpRejection);
   fhEventsProcessedNoPileUpRejection->Sumw2();
@@ -960,6 +978,37 @@ void AliAnalysisTaskPID::UserCreateOutputObjects()
     fQAContainer->Add(fDeDxCheck);
   }
   
+  if (fDoBinZeroStudy) {
+    const Double_t etaLow = -0.9;
+    const Double_t etaUp = 0.9;
+    const Int_t nEtaBins = 18;
+    
+    const Int_t nBinsBinZeroStudy = kBinZeroStudyNumAxes;
+    Int_t binZeroStudyBins[nBinsBinZeroStudy]    = { nCentBins,                   nPtBins, nEtaBins };
+    Double_t binZeroStudyXmin[nBinsBinZeroStudy] = { binsCent[0],               binsPt[0], etaLow   };
+    Double_t binZeroStudyXmax[nBinsBinZeroStudy] = { binsCent[nCentBins], binsPt[nPtBins], etaUp    };
+    
+    fChargedGenPrimariesTriggerSel = new THnSparseD("fChargedGenPrimariesTriggerSel", "Trigger sel.", nBinsBinZeroStudy, binZeroStudyBins,
+                                                    binZeroStudyXmin, binZeroStudyXmax);
+    SetUpBinZeroStudyHist(fChargedGenPrimariesTriggerSel, binsCent, binsPt);
+    fOutputContainer->Add(fChargedGenPrimariesTriggerSel);
+    
+    fChargedGenPrimariesTriggerSelVtxCut = new THnSparseD("fChargedGenPrimariesTriggerSelVtxCut", "Vertex cut", nBinsBinZeroStudy,
+                                                          binZeroStudyBins, binZeroStudyXmin, binZeroStudyXmax);
+    SetUpBinZeroStudyHist(fChargedGenPrimariesTriggerSelVtxCut, binsCent, binsPt);
+    fOutputContainer->Add(fChargedGenPrimariesTriggerSelVtxCut);
+    
+    fChargedGenPrimariesTriggerSelVtxCutZ = new THnSparseD("fChargedGenPrimariesTriggerSelVtxCutZ", "Vertex #it{z} cut", nBinsBinZeroStudy,
+                                                          binZeroStudyBins, binZeroStudyXmin, binZeroStudyXmax);
+    SetUpBinZeroStudyHist(fChargedGenPrimariesTriggerSelVtxCutZ, binsCent, binsPt);
+    fOutputContainer->Add(fChargedGenPrimariesTriggerSelVtxCutZ);
+    
+    fChargedGenPrimariesTriggerSelVtxCutZPileUpRej = new THnSparseD("fChargedGenPrimariesTriggerSelVtxCutZPileUpRej", "Vertex #it{z} cut", 
+                                                                    nBinsBinZeroStudy, binZeroStudyBins, binZeroStudyXmin, binZeroStudyXmax);
+    SetUpBinZeroStudyHist(fChargedGenPrimariesTriggerSelVtxCutZPileUpRej, binsCent, binsPt);
+    fOutputContainer->Add(fChargedGenPrimariesTriggerSelVtxCutZPileUpRej);
+  }
+  
   if(fDebug > 2)
     printf("File: %s, Line: %d: UserCreateOutputObjects -> Posting output data\n", (char*)__FILE__, __LINE__);
   
@@ -1012,40 +1061,81 @@ void AliAnalysisTaskPID::UserExec(Option_t *)
         AliError("Not esd event -> Cannot use tracklet multiplicity estimator!");
         centralityPercentile = -1;
       }
-      else {
+      else
         centralityPercentile = AliESDtrackCuts::GetReferenceMultiplicity(esdEvent, AliESDtrackCuts::kTrackletsITSTPC, fEtaAbsCutUp);
-      }
     }
-    else
+    else if (fCentralityEstimator.Contains("ppMult", TString::kIgnoreCase)) {
+      // Another special pp centrality estimator
+      centralityPercentile = fAnaUtils->GetMultiplicityPercentile(fEvent, GetPPCentralityEstimator().Data());
+    }
+    else {
+      // Ordinary centrality estimator
       centralityPercentile = fEvent->GetCentrality()->GetCentralityPercentile(fCentralityEstimator.Data());
+    }
   }
   
-  IncrementEventCounter(centralityPercentile, kTriggerSel);
-  
   // Check if vertex is ok, but don't apply cut on z position
-  if (!GetVertexIsOk(fEvent, kFALSE))
-    return;
+  const Bool_t passedVertexSelection = GetVertexIsOk(fEvent, kFALSE);
+  // Now check again, but also require z position to be in desired range
+  const Bool_t passedVertexZSelection = GetVertexIsOk(fEvent, kTRUE);
+  // Check pile-up
+  const Bool_t isPileUp = GetIsPileUp(fEvent, fPileUpRejectionType);
+  
   
-  fESD = dynamic_cast<AliESDEvent*>(fEvent);
-  const AliVVertex* primaryVertex = fESD ? fESD->GetPrimaryVertexTracks() : fEvent->GetPrimaryVertex(); 
-  if (!primaryVertex)
-    return;
   
-  if(primaryVertex->GetNContributors() <= 0) 
+  if (fDoBinZeroStudy && fMC) {
+    for (Int_t iPart = 0; iPart < fMC->GetNumberOfTracks(); iPart++) { 
+      AliMCParticle *mcPart  = dynamic_cast<AliMCParticle*>(fMC->GetTrack(iPart));
+      
+      if (!mcPart)
+          continue;
+      
+      if (!fMC->IsPhysicalPrimary(iPart)) 
+          continue;
+      
+      const Double_t etaGen = mcPart->Eta();
+      const Double_t ptGen = mcPart->Pt();
+      
+      Double_t values[kBinZeroStudyNumAxes] = { 0. };
+      values[kBinZeroStudyCentrality] = centralityPercentile;
+      values[kBinZeroStudyGenPt] = ptGen;
+      values[kBinZeroStudyGenEta] = etaGen;
+      
+      fChargedGenPrimariesTriggerSel->Fill(values);
+      if (passedVertexSelection) {
+          fChargedGenPrimariesTriggerSelVtxCut->Fill(values);
+        if (passedVertexZSelection) {
+            fChargedGenPrimariesTriggerSelVtxCutZ->Fill(values);
+          if (!isPileUp) {
+            fChargedGenPrimariesTriggerSelVtxCutZPileUpRej->Fill(values);
+          }
+        }
+      }
+    }
+  }
+  
+  
+  
+  // Event counters for trigger selection, vertex cuts and pile-up rejection
+  IncrementEventCounter(centralityPercentile, kTriggerSel);
+  
+  if (!passedVertexSelection)
     return;
   
   IncrementEventCounter(centralityPercentile, kTriggerSelAndVtxCut);
   
-  // Now check again, but also require z position to be in desired range
-  if (!GetVertexIsOk(fEvent, kTRUE))
+  if (!passedVertexZSelection)
     return;
   
   IncrementEventCounter(centralityPercentile, kTriggerSelAndVtxCutAndZvtxCutNoPileUpRejection);
-  //TODO ATTENTION: Is this the right place for the pile-up rejection? Important to have still the proper bin-0 correction,
+  // ATTENTION: Is this the right place for the pile-up rejection? Important to have still the proper bin-0 correction,
   // which is done solely with sel and selVtx, since the zvtx selection does ~not change the spectra. The question is whether the pile-up
   // rejection changes the spectra. If not, then it is perfectly fine to put it here and keep the usual histo for the normalisation to number
   // of events. But if it does change the spectra, this must somehow be corrected for.
-  if (GetIsPileUp(fEvent, fPileUpRejectionType))
+  // NOTE: multiplicity >= 0 usually implies a properly reconstructed vertex. Hence, the bin-0 correction cannot be done in multiplicity bins.
+  // Furthermore, there seems to be no MC simulation with pile-up rejection, so the bin-0 correction cannot be extracted with it. Pile-up
+  // rejection has only a minor impact, so maybe there is no need to dig further.
+  if (isPileUp)
     return;
   
   IncrementEventCounter(centralityPercentile, kTriggerSelAndVtxCutAndZvtxCut);
@@ -1787,6 +1877,7 @@ Double_t AliAnalysisTaskPID::GetMCStrangenessFactorCMS(Int_t motherPDG, Double_t
   }
   
   if (absMotherPDG == 3122) { // Lambda
+  //if (absMotherPDG == 3122 || absMotherPDG == 3112 || absMotherPDG == 3222) { // Lambda / Sigma- / Sigma+
     if (0.00 <= motherGenPt && motherGenPt < 0.20) fac = 0.645162;
     else if(0.20 <= motherGenPt && motherGenPt < 0.40) fac = 0.627431;
     else if(0.40 <= motherGenPt && motherGenPt < 0.60) fac = 0.457136;
@@ -2275,6 +2366,7 @@ void AliAnalysisTaskPID::PrintSettings(Bool_t printSystematicsSettings) const
   printf("Do Efficiency: %d\n", fDoEfficiency);
   printf("Do PtResolution: %d\n", fDoPtResolution);
   printf("Do dEdxCheck: %d\n", fDoDeDxCheck);
+  printf("Do binZeroStudy: %d\n", fDoBinZeroStudy);
   
   printf("\n");
 
@@ -3524,7 +3616,7 @@ void AliAnalysisTaskPID::SetUpGenHist(THnSparse* hist, Double_t* binsPt, Double_
   
   hist->GetAxis(kGenDeltaPrimeSpecies)->SetTitle("TPC #Delta'_{species} (arb. unit)");
   
-  hist->GetAxis(kGenCentrality)->SetTitle(Form("Centrality Percentile (%s)", fCentralityEstimator.Data()));
+  hist->GetAxis(kGenCentrality)->SetTitle(Form("Centrality Percentile (%s)", GetPPCentralityEstimator().Data()));
   
   if (fStoreAdditionalJetInformation) {
     hist->GetAxis(kGenJetPt)->SetTitle("p_{T}^{jet} (GeV/c)");
@@ -3562,7 +3654,7 @@ void AliAnalysisTaskPID::SetUpGenYieldHist(THnSparse* hist, Double_t* binsPt, Do
   // Set axes titles
   hist->GetAxis(kGenYieldMCID)->SetTitle("MC PID");
   hist->GetAxis(kGenYieldPt)->SetTitle("p_{T}^{gen} (GeV/c)");
-  hist->GetAxis(kGenYieldCentrality)->SetTitle(Form("Centrality Percentile (%s)", fCentralityEstimator.Data()));
+  hist->GetAxis(kGenYieldCentrality)->SetTitle(Form("Centrality Percentile (%s)", GetPPCentralityEstimator().Data()));
   
   if (fStoreAdditionalJetInformation) {
     hist->GetAxis(kGenYieldJetPt)->SetTitle("p_{T}^{jet, gen} (GeV/c)");
@@ -3606,7 +3698,7 @@ void AliAnalysisTaskPID::SetUpHist(THnSparse* hist, Double_t* binsPt, Double_t*
     
   hist->GetAxis(kDataDeltaPrimeSpecies)->SetTitle("TPC #Delta'_{species} (arb. unit)");
   
-  hist->GetAxis(kDataCentrality)->SetTitle(Form("Centrality Percentile (%s)", fCentralityEstimator.Data()));
+  hist->GetAxis(kDataCentrality)->SetTitle(Form("Centrality Percentile (%s)", GetPPCentralityEstimator().Data()));
   
   if (fStoreAdditionalJetInformation) {
     hist->GetAxis(kDataJetPt)->SetTitle("p_{T}^{jet} (GeV/c)");
@@ -3644,7 +3736,7 @@ void AliAnalysisTaskPID::SetUpPtResHist(THnSparse* hist, Double_t* binsPt, Doubl
   hist->GetAxis(kPtResRecPt)->SetTitle("p_{T}^{rec} (GeV/c)");  
   
   hist->GetAxis(kPtResCharge)->SetTitle("Charge (e_{0})");
-  hist->GetAxis(kPtResCentrality)->SetTitle(Form("Centrality Percentile (%s)", fCentralityEstimator.Data()));
+  hist->GetAxis(kPtResCentrality)->SetTitle(Form("Centrality Percentile (%s)", GetPPCentralityEstimator().Data()));
 }
 
 
@@ -3686,4 +3778,18 @@ void AliAnalysisTaskPID::SetUpDeDxCheckHist(THnSparse* hist, const Double_t* bin
   hist->GetAxis(kDeDxCheckEtaAbs)->SetTitle("|#eta|");  
   hist->GetAxis(kDeDxCheckP)->SetTitle("p_{TPC} (GeV/c)"); 
   hist->GetAxis(kDeDxCheckDeDx)->SetTitle("TPC dE/dx (arb. unit)");
+}
+
+
+//________________________________________________________________________
+void AliAnalysisTaskPID::SetUpBinZeroStudyHist(THnSparse* hist, const Double_t* binsCent, const Double_t* binsPt) const
+{
+  // Sets bin limits for axes which are not standard binned and the axes titles.
+  hist->SetBinEdges(kBinZeroStudyCentrality, binsCent);
+  hist->SetBinEdges(kBinZeroStudyGenPt, binsPt);
+  
+  // Set axes titles
+  hist->GetAxis(kBinZeroStudyCentrality)->SetTitle(Form("Centrality Percentile (%s)", GetPPCentralityEstimator().Data()));
+  hist->GetAxis(kBinZeroStudyGenEta)->SetTitle("#it{#eta}^{gen}");  
+  hist->GetAxis(kBinZeroStudyGenPt)->SetTitle("#it{p}_{T}^{gen} (GeV/#it{c})"); 
 }
\ No newline at end of file