- Added monitoring output to check jetPt dependence of dEdx - Final bug fix: Introduc...

diff --git a/PWGJE/UserTasks/AliAnalysisTaskPID.cxx b/PWGJE/UserTasks/AliAnalysisTaskPID.cxx
index ec6218f12dc5594276ccbd291a7bcc1d49f845d2..b60f08a6091326ead10ef6adef495d9a45afe862 100644 (file)
--- a/PWGJE/UserTasks/AliAnalysisTaskPID.cxx
+++ b/PWGJE/UserTasks/AliAnalysisTaskPID.cxx
@@ -54,6 +54,7 @@ AliAnalysisTaskPID::AliAnalysisTaskPID()
   , fDoPID(kTRUE)
   , fDoEfficiency(kTRUE)
   , fDoPtResolution(kTRUE)
+  , fDoDeDxCheck(kTRUE)
   , fStoreCentralityPercentile(kFALSE)
   , fStoreAdditionalJetInformation(kFALSE)
   , fTakeIntoAccountMuons(kFALSE)
@@ -131,6 +132,7 @@ AliAnalysisTaskPID::AliAnalysisTaskPID()
   , fGenRespPrDeltaPi(new Double_t[fgkMaxNumGenEntries])
   , fGenRespPrDeltaPr(new Double_t[fgkMaxNumGenEntries])
   */
+  , fDeltaPrimeAxis(0x0)
   , fhMaxEtaVariation(0x0)
   , fhEventsProcessed(0x0)
   , fhEventsTriggerSel(0x0)
@@ -142,8 +144,9 @@ AliAnalysisTaskPID::AliAnalysisTaskPID()
   , fh1Xsec(0x0)
   , fh1Trials(0x0)
   , fContainerEff(0x0)
+  , fDeDxCheck(0x0)
   , fOutputContainer(0x0)
-  , fPtResolutionContainer(0x0)
+  , fQAContainer(0x0)
 {
   // default Constructor
   
@@ -185,6 +188,7 @@ AliAnalysisTaskPID::AliAnalysisTaskPID(const char *name)
   , fDoPID(kTRUE)
   , fDoEfficiency(kTRUE)
   , fDoPtResolution(kTRUE)
+  , fDoDeDxCheck(kTRUE)
   , fStoreCentralityPercentile(kFALSE)
   , fStoreAdditionalJetInformation(kFALSE)
   , fTakeIntoAccountMuons(kFALSE)
@@ -262,6 +266,7 @@ AliAnalysisTaskPID::AliAnalysisTaskPID(const char *name)
   , fGenRespPrDeltaPi(new Double_t[fgkMaxNumGenEntries])
   , fGenRespPrDeltaPr(new Double_t[fgkMaxNumGenEntries])
   */
+  , fDeltaPrimeAxis(0x0)
   , fhMaxEtaVariation(0x0)
   , fhEventsProcessed(0x0)
   , fhEventsTriggerSel(0x0)
@@ -273,8 +278,9 @@ AliAnalysisTaskPID::AliAnalysisTaskPID(const char *name)
   , fh1Xsec(0x0)
   , fh1Trials(0x0)
   , fContainerEff(0x0)
+  , fDeDxCheck(0x0)
   , fOutputContainer(0x0)
-  , fPtResolutionContainer(0x0)
+  , fQAContainer(0x0)
 {
   // Constructor
   
@@ -328,12 +334,15 @@ AliAnalysisTaskPID::~AliAnalysisTaskPID()
   delete fOutputContainer;
   fOutputContainer = 0x0;
   
-  delete fPtResolutionContainer;
-  fPtResolutionContainer = 0x0;
+  delete fQAContainer;
+  fQAContainer = 0x0;
 
   delete fConvolutedGausDeltaPrime;
   fConvolutedGausDeltaPrime = 0x0;
   
+  delete fDeltaPrimeAxis;
+  fDeltaPrimeAxis = 0x0;
+  
   delete [] fGenRespElDeltaPrimeEl;
   delete [] fGenRespElDeltaPrimeKa;
   delete [] fGenRespElDeltaPrimePi;
@@ -446,7 +455,7 @@ void AliAnalysisTaskPID::SetUpPIDcombined()
 {
   // Initialise the PIDcombined object
   
-  if (!fDoPID)
+  if (!fDoPID && !fDoDeDxCheck)
     return;
   
   if(fDebug > 1)
@@ -613,6 +622,8 @@ void AliAnalysisTaskPID::UserCreateOutputObjects()
     deltaPrimeBins[i] = factor * deltaPrimeBins[i - 1];
   }
   
+  fDeltaPrimeAxis = new TAxis(deltaPrimeNBins, deltaPrimeBins);
+  
   const Int_t nMCPIDbins = 5;
   const Double_t mcPIDmin = 0.;
   const Double_t mcPIDmax = 5.;
@@ -874,16 +885,17 @@ void AliAnalysisTaskPID::UserCreateOutputObjects()
   fOutputContainer->Add(fh1Xsec);
   fOutputContainer->Add(fh1Trials);
   
-  if (fDoPtResolution) {
+  if (fDoDeDxCheck || fDoPtResolution) {
     OpenFile(3);
+    fQAContainer = new TObjArray(1);
+    fQAContainer->SetName(Form("%s_QA", GetName()));
+    fQAContainer->SetOwner(kTRUE);
     
     if(fDebug > 2)
-    printf("File: %s, Line: %d: UserCreateOutputObjects -> OpenFile(3) successful\n", (char*)__FILE__, __LINE__);
-    
-    fPtResolutionContainer = new TObjArray(1);
-    fPtResolutionContainer->SetName(Form("%s_PtResolution", GetName()));
-    fPtResolutionContainer->SetOwner(kTRUE);
-    
+      printf("File: %s, Line: %d: UserCreateOutputObjects -> OpenFile(3) successful\n", (char*)__FILE__, __LINE__);
+  }
+  
+  if (fDoPtResolution) {
     const Int_t nPtBinsRes = 100;
     Double_t pTbinsRes[nPtBinsRes + 1];
 
@@ -909,16 +921,35 @@ void AliAnalysisTaskPID::UserCreateOutputObjects()
                                         Form("Pt resolution for primaries, %s", AliPID::ParticleLatexName(i)),
                                         nBinsPtResolution, ptResolutionBins, ptResolutionXmin, ptResolutionXmax);
       SetUpPtResHist(fPtResolution[i], pTbinsRes, binsJetPt, binsCent);
-      fPtResolutionContainer->Add(fPtResolution[i]);
+      fQAContainer->Add(fPtResolution[i]);
     }
   }
   
+  
+  
+  if (fDoDeDxCheck) {
+    const Int_t nEtaAbsBins = 9;
+    const Double_t binsEtaAbs[nEtaAbsBins+1] = { 0.0, 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9 };
+    
+    const Double_t dEdxMin = 20;
+    const Double_t dEdxMax = 110;
+    const Int_t nDeDxBins = (Int_t) ((dEdxMax - dEdxMin) / 0.02);
+    const Int_t nBinsDeDxCheck= kDeDxCheckNumAxes;
+    Int_t dEdxCheckBins[kDeDxCheckNumAxes]    = { nSelSpeciesBins, nPtBins, nJetPtBins, nEtaAbsBins, nDeDxBins };
+    Double_t dEdxCheckXmin[kDeDxCheckNumAxes] = { selSpeciesMin, binsPt[0], binsJetPt[0], binsEtaAbs[0], dEdxMin };
+    Double_t dEdxCheckXmax[kDeDxCheckNumAxes] = { selSpeciesMax, binsPt[nPtBins], binsJetPt[nJetPtBins], binsEtaAbs[nEtaAbsBins], dEdxMax };
+    
+    fDeDxCheck = new THnSparseD("fDeDxCheck", "dEdx check", nBinsDeDxCheck, dEdxCheckBins, dEdxCheckXmin, dEdxCheckXmax);
+    SetUpDeDxCheckHist(fDeDxCheck, binsPt, binsJetPt, binsEtaAbs);
+    fQAContainer->Add(fDeDxCheck);
+  }
+  
   if(fDebug > 2)
     printf("File: %s, Line: %d: UserCreateOutputObjects -> Posting output data\n", (char*)__FILE__, __LINE__);
   
   PostData(1, fOutputContainer);
   PostData(2, fContainerEff);
-  PostData(3, fPtResolutionContainer);
+  PostData(3, fQAContainer);
   
   if(fDebug > 2)
     printf("File: %s, Line: %d: UserCreateOutputObjects -> Done\n", (char*)__FILE__, __LINE__);
@@ -1112,7 +1143,7 @@ void AliAnalysisTaskPID::UserExec(Option_t *)
     // Only process tracks inside the desired eta window    
     if (!IsInAcceptedEtaRange(TMath::Abs(track->Eta()))) continue;
    
-    if (fDoPID) 
+    if (fDoPID || fDoDeDxCheck) 
       ProcessTrack(track, pdg, centralityPercentile, -1); // No jet information in this case -> Set jet pT to -1
     
     if (fDoPtResolution) {
@@ -1410,7 +1441,7 @@ Bool_t AliAnalysisTaskPID::GetParticleFraction(Double_t trackPt, Double_t jetPt,
   else {
     Double_t x0 = 0., x1 = 0., y0 = 0., y1 = 0.;
     Double_t y0errStat = 0., y1errStat = 0., y0errSys = 0., y1errSys = 0.;
-    Int_t trackPtBin = xAxis->FindBin(trackPt);
+    Int_t trackPtBin = xAxis->FindFixBin(trackPt);
     
     // Linear interpolation between nearest neighbours in trackPt
     if (trackPt <= xAxis->GetBinCenter(trackPtBin)) {
@@ -2142,6 +2173,7 @@ void AliAnalysisTaskPID::PrintSettings(Bool_t printSystematicsSettings) const
   printf("Do PID: %d\n", fDoPID);
   printf("Do Efficiency: %d\n", fDoEfficiency);
   printf("Do PtResolution: %d\n", fDoPtResolution);
+  printf("Do dEdxCheck: %d\n", fDoDeDxCheck);
   
   printf("\n");
 
@@ -2191,7 +2223,7 @@ Bool_t AliAnalysisTaskPID::ProcessTrack(const AliVTrack* track, Int_t particlePD
   if(fDebug > 1)
     printf("File: %s, Line: %d: ProcessTrack\n", (char*)__FILE__, __LINE__);
   
-  if (!fDoPID)
+  if (!fDoPID && !fDoDeDxCheck)
     return kFALSE;
   
   if(fDebug > 2)
@@ -2504,7 +2536,7 @@ Bool_t AliAnalysisTaskPID::ProcessTrack(const AliVTrack* track, Int_t particlePD
     
   }
   else {
-    // No systematic studies on expected signal - just take it as it comve from the TPCPIDResponse
+    // No systematic studies on expected signal - just take it as it comes from the TPCPIDResponse
     dEdxEl = fPIDResponse->GetTPCResponse().GetExpectedSignal(track, AliPID::kElectron, AliTPCPIDResponse::kdEdxDefault, 
                                                               fPIDResponse->UseTPCEtaCorrection(),
                                                               fPIDResponse->UseTPCMultiplicityCorrection());
@@ -2640,6 +2672,69 @@ Bool_t AliAnalysisTaskPID::ProcessTrack(const AliVTrack* track, Int_t particlePD
     // In all other cases, the maximum remains untouched from the scaling (and is < AliPID::kSPECIES by construction)
   }
   
+  if (fDoDeDxCheck) {
+    // Generate the expected responses in DeltaPrime and translate to real dEdx. Only take responses for exactly that species
+    // (i.e. with pre-PID)
+    
+    Int_t numGenEntries = fgkMaxNumGenEntries; // fgkMaxNumGenEntries = 500
+    
+    ErrorCode errCode = kNoErrors;
+  
+    if(fDebug > 2)
+      printf("File: %s, Line: %d: ProcessTrack -> Generate Responses for dEdx check\n", (char*)__FILE__, __LINE__);
+    
+    // Electrons
+    errCode = GenerateDetectorResponse(errCode, 1., sigmaEl / dEdxEl, fGenRespElDeltaPrimeEl, numGenEntries);
+
+    // Kaons
+    errCode = GenerateDetectorResponse(errCode, 1., sigmaKa / dEdxKa, fGenRespKaDeltaPrimeKa, numGenEntries);
+
+    // Pions
+    errCode = GenerateDetectorResponse(errCode, 1., sigmaPi / dEdxPi, fGenRespPiDeltaPrimePi, numGenEntries);
+
+    // Protons
+    errCode = GenerateDetectorResponse(errCode, 1., sigmaPr / dEdxPr, fGenRespPrDeltaPrimePr, numGenEntries);
+    
+    if (errCode == kNoErrors) {
+      Double_t genEntry[kDeDxCheckNumAxes];
+      genEntry[kDeDxCheckJetPt] = jetPt;
+      genEntry[kDeDxCheckEtaAbs] = TMath::Abs(track->Eta());
+      genEntry[kDeDxCheckP] = pTPC;
+      
+        
+      for (Int_t n = 0; n < numGenEntries; n++)  {
+        // If no MC info is available or shall not be used, use weighting with priors to generate entries for the different species
+        Double_t rnd = fRandom->Rndm(); // Produce uniformly distributed floating point in ]0, 1]
+        
+        // Consider generated response as originating from...
+        if (rnd <= prob[AliPID::kElectron]) {
+          genEntry[kDeDxCheckPID] = 0; // ... an electron
+          genEntry[kDeDxCheckDeDx] = fGenRespElDeltaPrimeEl[n] * dEdxEl;
+        }
+        else if (rnd <= prob[AliPID::kElectron] + prob[AliPID::kKaon]) {
+          genEntry[kDeDxCheckPID] = 1;  // ... a kaon
+          genEntry[kDeDxCheckDeDx] = fGenRespKaDeltaPrimeKa[n] * dEdxKa;
+        }
+        else if (rnd <= prob[AliPID::kElectron] + prob[AliPID::kKaon] + prob[AliPID::kMuon] + prob[AliPID::kPion]) {
+          genEntry[kDeDxCheckPID] = 2;  // ... a pion (or a muon)
+          genEntry[kDeDxCheckDeDx] = fGenRespPiDeltaPrimePi[n] * dEdxPi;
+        }
+        else {
+          genEntry[kDeDxCheckPID] = 3;  // ... a proton
+          genEntry[kDeDxCheckDeDx] = fGenRespPrDeltaPrimePr[n] * dEdxPr;
+        }
+     
+        fDeDxCheck->Fill(genEntry);
+      }
+    }
+    
+    if(fDebug > 2)
+      printf("File: %s, Line: %d: ProcessTrack -> Generate Responses for dEdx check done\n", (char*)__FILE__, __LINE__);
+  }
+  
+  if (!fDoPID)
+    return kTRUE;
+  
   if (!isMC) {
     // Clean up the most probable PID for low momenta (not an issue for the probabilities, since fractions small,
     // but to get clean (and nice looking) templates (if most probable PID is used instead of generated responses), it should be done).
@@ -3259,8 +3354,11 @@ AliAnalysisTaskPID::ErrorCode AliAnalysisTaskPID::SetParamsForConvolutedGaus(Dou
     rangeEnd = fConvolutedGausDeltaPrime->GetX(maximumFraction, upBoundSearchBoundLow, upBoundSearchBoundUp, (maxX < 0.4) ? 1e-5 : 0.001);
     
     fConvolutedGausDeltaPrime->SetRange(rangeStart,rangeEnd);
-    fConvolutedGausDeltaPrime->SetNpx(fhPIDdataAll->GetAxis(kDataDeltaPrimeSpecies)->FindBin(rangeEnd)
-                                      - fhPIDdataAll->GetAxis(kDataDeltaPrimeSpecies)->FindBin(rangeStart) + 1);
+    fConvolutedGausDeltaPrime->SetNpx(fDeltaPrimeAxis->FindFixBin(rangeEnd) - fDeltaPrimeAxis->FindFixBin(rangeStart) + 1);
+    
+    fConvolutedGausDeltaPrime->SetNpx(fDeltaPrimeAxis->FindFixBin(rangeEnd) - fDeltaPrimeAxis->FindFixBin(rangeStart) + 1);
+    //fConvolutedGausDeltaPrime->SetNpx(fhPIDdataAll->GetAxis(kDataDeltaPrimeSpecies)->FindFixBin(rangeEnd)
+    //                                  - fhPIDdataAll->GetAxis(kDataDeltaPrimeSpecies)->FindFixBin(rangeStart) + 1);
     //fConvolutedGausDeltaPrime->SetNpx((rangeEnd - rangeStart) / fDeltaPrimeBinWidth + 1);
   }
   
@@ -3424,3 +3522,26 @@ void AliAnalysisTaskPID::SetUpPtResHist(THnSparse* hist, Double_t* binsPt, Doubl
   hist->GetAxis(kPtResCharge)->SetTitle("Charge (e_{0})");
   hist->GetAxis(kPtResCentrality)->SetTitle(Form("Centrality Percentile (%s)", fCentralityEstimator.Data()));
 }
+
+//________________________________________________________________________
+void AliAnalysisTaskPID::SetUpDeDxCheckHist(THnSparse* hist, const Double_t* binsPt, const Double_t* binsJetPt, const Double_t* binsEtaAbs) const
+{
+  // Sets bin limits for axes which are not standard binned and the axes titles.
+  hist->SetBinEdges(kDeDxCheckP, binsPt);
+  hist->SetBinEdges(kDeDxCheckJetPt, binsJetPt);
+  hist->SetBinEdges(kDeDxCheckEtaAbs, binsEtaAbs);
+  
+  
+  // Set axes titles
+  hist->GetAxis(kDeDxCheckPID)->SetTitle("PID");
+  hist->GetAxis(kDeDxCheckPID)->SetBinLabel(1, "e");
+  hist->GetAxis(kDeDxCheckPID)->SetBinLabel(2, "K");
+  hist->GetAxis(kDeDxCheckPID)->SetBinLabel(3, "#pi");
+  hist->GetAxis(kDeDxCheckPID)->SetBinLabel(4, "p");
+  
+  
+  hist->GetAxis(kDeDxCheckJetPt)->SetTitle("p_{T}^{jet} (GeV/c)");
+  hist->GetAxis(kDeDxCheckEtaAbs)->SetTitle("|#eta|");  
+  hist->GetAxis(kDeDxCheckP)->SetTitle("p_{TPC} (GeV/c)"); 
+  hist->GetAxis(kDeDxCheckDeDx)->SetTitle("TPC dE/dx (arb. unit)");
+}
\ No newline at end of file

diff --git a/PWGJE/UserTasks/AliAnalysisTaskPID.h b/PWGJE/UserTasks/AliAnalysisTaskPID.h
index a063ff9472fc50507a9d047cf826d09ba28060ab..7ec7640e8846b1aded23964ef3328e2c62b729a6 100644 (file)
--- a/PWGJE/UserTasks/AliAnalysisTaskPID.h
+++ b/PWGJE/UserTasks/AliAnalysisTaskPID.h
@@ -25,6 +25,7 @@ class AliTOFPIDResponse;
 class AliVEvent;
 class AliVTrack;
 
+#include "TAxis.h"
 #include "TH1D.h"
 #include "TH2D.h"
 #include "TH3D.h"
@@ -61,6 +62,9 @@ class AliAnalysisTaskPID : public AliAnalysisTaskPIDV0base {
   
   enum ptResolutionAxes { kPtResJetPt = 0, kPtResGenPt = 1, kPtResRecPt = 2, kPtResCharge = 3, kPtResCentrality = 4, kPtResNumAxes = 5 };
   
+  enum dEdxCheckAxes { kDeDxCheckPID = 0, kDeDxCheckP = 1, kDeDxCheckJetPt = 2, kDeDxCheckEtaAbs = 3 , kDeDxCheckDeDx = 4,
+                       kDeDxCheckNumAxes = 5 };
+  
   enum efficiencyAxes { kEffMCID = 0, kEffTrackPt = 1, kEffTrackEta = 2, kEffTrackCharge = 3, kEffCentrality = 4, kEffJetPt = 5,
                         kEffZ = 6, kEffXi = 7, kEffNumAxes = 8 };
   
@@ -142,6 +146,9 @@ class AliAnalysisTaskPID : public AliAnalysisTaskPIDV0base {
   Bool_t GetDoPtResolution() const { return fDoPtResolution; };
   void SetDoPtResolution(Bool_t flag) { fDoPtResolution = flag; };
   
+  Bool_t GetDoDeDxCheck() const { return fDoDeDxCheck; };
+  void SetDoDeDxCheck(Bool_t flag) { fDoDeDxCheck = flag; };
+  
   Bool_t GetStoreCentralityPercentile() const { return fStoreCentralityPercentile; };
   void SetStoreCentralityPercentile(Bool_t flag) { fStoreCentralityPercentile = flag; };
   
@@ -253,6 +260,7 @@ class AliAnalysisTaskPID : public AliAnalysisTaskPIDV0base {
   virtual void SetUpGenYieldHist(THnSparse* hist, Double_t* binsPt, Double_t* binsCent, Double_t* binsJetPt) const;
   virtual void SetUpHist(THnSparse* hist, Double_t* binsPt, Double_t* binsDeltaPrime, Double_t* binsCent, Double_t* binsJetPt) const;
   virtual void SetUpPtResHist(THnSparse* hist, Double_t* binsPt, Double_t* binsJetPt, Double_t* binsCent) const;
+  virtual void SetUpDeDxCheckHist(THnSparse* hist, const Double_t* binsPt, const Double_t* binsJetPt, const Double_t* binsEtaAbs) const;
   virtual void SetUpPIDcombined();
   
   static const Int_t fgkNumJetAxes; // Number of additional axes for jets
@@ -270,6 +278,7 @@ class AliAnalysisTaskPID : public AliAnalysisTaskPIDV0base {
   Bool_t fDoPID; // Only do PID processing (and post the output), if flag is set to kTRUE
   Bool_t fDoEfficiency; // Only do efficiency processing (and post the output), if flag is set to kTRUE
   Bool_t fDoPtResolution; // Only do pT resolution processing (and post the output), if flag is set to kTRUE
+  Bool_t fDoDeDxCheck; // Only check dEdx, if flag set to kTRUE
   
   Bool_t fStoreCentralityPercentile; // If set to kTRUE, store centrality percentile for each event. In case of kFALSE (appropriate for pp), centrality percentile will be set to -1 for every event
   Bool_t fStoreAdditionalJetInformation; // If set to kTRUE, additional jet information like jetPt, z, xi will be stored in the THnSparses
@@ -368,6 +377,7 @@ class AliAnalysisTaskPID : public AliAnalysisTaskPIDV0base {
   Double_t* fGenRespPrDeltaPr; //! Generated responses for a single track
   */
   
+  TAxis* fDeltaPrimeAxis; //! Axis holding the deltaPrime binning
   TH1D* fhMaxEtaVariation; //! Histo holding the maximum deviation of the eta correction factor from unity vs. 1/dEdx(splines)
   
   TH1D* fhEventsProcessed; //! Histo holding the number of processed events (i.e. passing trigger selection, vtx and zvtx cuts
@@ -387,14 +397,16 @@ class AliAnalysisTaskPID : public AliAnalysisTaskPIDV0base {
   
   THnSparseD* fPtResolution[AliPID::kSPECIES]; //! Pt Resolution for the individual species
   
+  THnSparseD* fDeDxCheck; //! dEdx check
+  
   TObjArray* fOutputContainer;  //! output data container
   
-  TObjArray* fPtResolutionContainer;  //! output data container for pt resolution
+  TObjArray* fQAContainer; //! output data container for QA
   
   AliAnalysisTaskPID(const AliAnalysisTaskPID&); // not implemented
   AliAnalysisTaskPID& operator=(const AliAnalysisTaskPID&); // not implemented
   
-  ClassDef(AliAnalysisTaskPID, 18);
+  ClassDef(AliAnalysisTaskPID, 19);
 };
 
 
@@ -606,8 +618,8 @@ inline void AliAnalysisTaskPID::PostOutputData()
   if (fDoEfficiency)
     PostData(2, fContainerEff);
   
-  if (fDoPtResolution)
-    PostData(3, fPtResolutionContainer);
+  if (fDoPtResolution || fDoDeDxCheck)
+    PostData(3, fQAContainer);
 }
 
 #endif