added histograms for stats about the vertex, and delta phi of tracklets
authorjgrosseo <jgrosseo@f7af4fe6-9843-0410-8265-dc069ae4e863>
Tue, 29 Jan 2008 17:42:27 +0000 (17:42 +0000)
committerjgrosseo <jgrosseo@f7af4fe6-9843-0410-8265-dc069ae4e863>
Tue, 29 Jan 2008 17:42:27 +0000 (17:42 +0000)
PWG0/dNdEta/AlidNdEtaCorrectionTask.cxx
PWG0/dNdEta/AlidNdEtaCorrectionTask.h
PWG0/dNdEta/AlidNdEtaTask.cxx
PWG0/dNdEta/AlidNdEtaTask.h

index ae21d29..4ba3692 100644 (file)
@@ -2,10 +2,13 @@
 
 #include "AlidNdEtaCorrectionTask.h"
 
+#include <TROOT.h>
 #include <TCanvas.h>
 #include <TChain.h>
 #include <TFile.h>
 #include <TH1F.h>
+#include <TH2F.h>
+#include <TProfile.h>
 #include <TParticle.h>
 
 #include <AliLog.h>
@@ -26,6 +29,7 @@
 //#include "AliCorrectionMatrix3D.h"
 #include "dNdEta/dNdEtaAnalysis.h"
 #include "dNdEta/AlidNdEtaCorrection.h"
+#include "AliVertexerTracks.h"
 
 ClassImp(AlidNdEtaCorrectionTask)
 
@@ -42,6 +46,9 @@ AlidNdEtaCorrectionTask::AlidNdEtaCorrectionTask(const char* opt) :
   fdNdEtaAnalysisESD(0),
   fPIDParticles(0),
   fPIDTracks(0),
+  fVertexCorrelation(0),
+  fVertexProfile(0),
+  fVertexShiftNorm(0),
   fSigmaVertexTracks(0),
   fSigmaVertexPrim(0)
 {
@@ -145,9 +152,9 @@ void AlidNdEtaCorrectionTask::CreateOutputObjects()
   fOutput->Add(fdNdEtaAnalysisESD);
 
   if (fOption.Contains("process-types")) {
-    fdNdEtaCorrectionProcessType[0] = new AlidNdEtaCorrection("dndeta_correction_ND", "dndeta_correction_ND");
-    fdNdEtaCorrectionProcessType[1] = new AlidNdEtaCorrection("dndeta_correction_SD", "dndeta_correction_SD");
-    fdNdEtaCorrectionProcessType[2] = new AlidNdEtaCorrection("dndeta_correction_DD", "dndeta_correction_DD");
+    fdNdEtaCorrectionProcessType[0] = new AlidNdEtaCorrection("dndeta_correction_ND", "dndeta_correction_ND", fAnalysisMode);
+    fdNdEtaCorrectionProcessType[1] = new AlidNdEtaCorrection("dndeta_correction_SD", "dndeta_correction_SD", fAnalysisMode);
+    fdNdEtaCorrectionProcessType[2] = new AlidNdEtaCorrection("dndeta_correction_DD", "dndeta_correction_DD", fAnalysisMode);
 
     fOutput->Add(fdNdEtaCorrectionProcessType[0]);
     fOutput->Add(fdNdEtaCorrectionProcessType[1]);
@@ -162,6 +169,10 @@ void AlidNdEtaCorrectionTask::CreateOutputObjects()
     fOutput->Add(fSigmaVertexPrim);
     Printf("WARNING: sigma-vertex analysis enabled. This will produce weird results in the AliESDtrackCuts histograms");
   }
+
+  fVertexCorrelation = new TH2F("fVertexCorrelation", "fVertexCorrelation;MC z-vtx;ESD z-vtx", 80, -20, 20, 80, -20, 20);
+  fVertexProfile = new TProfile("fVertexProfile", "fVertexProfile;MC z-vtx;MC z-vtx - ESD z-vtx", 40, -20, 20);
+  fVertexShiftNorm = new TH1F("fVertexShiftNorm", "fVertexShiftNorm;(MC z-vtx - ESD z-vtx) / #sigma_{ESD z-vtx};Entries", 200, -100, 100);
 }
 
 void AlidNdEtaCorrectionTask::Exec(Option_t*)
@@ -176,14 +187,67 @@ void AlidNdEtaCorrectionTask::Exec(Option_t*)
   }
 
   // trigger definition
-  Bool_t eventTriggered = AliPWG0Helper::IsEventTriggered(fESD->GetTriggerMask(), AliPWG0Helper::kMB1);
+  Bool_t eventTriggered = AliPWG0Helper::IsEventTriggered(fESD->GetTriggerMask(), AliPWG0Helper::kMB2);
+
+  // post the data already here
+  PostData(0, fOutput);
+
+  // MC info
+  AliMCEventHandler* eventHandler = dynamic_cast<AliMCEventHandler*> (AliAnalysisManager::GetAnalysisManager()->GetMCtruthEventHandler());
+  if (!eventHandler) {
+    Printf("ERROR: Could not retrieve MC event handler");
+    return;
+  }
+
+  AliMCEvent* mcEvent = eventHandler->MCEvent();
+  if (!mcEvent) {
+    Printf("ERROR: Could not retrieve MC event");
+    return;
+  }
+
+  AliStack* stack = mcEvent->Stack();
+  if (!stack)
+  {
+    AliDebug(AliLog::kError, "Stack not available");
+    return;
+  }
+
+  AliHeader* header = mcEvent->Header();
+  if (!header)
+  {
+    AliDebug(AliLog::kError, "Header not available");
+    return;
+  }
+
+  // get process type; NB: this only works for Pythia
+  Int_t processType = AliPWG0Helper::GetPythiaEventProcessType(header);
+  AliDebug(AliLog::kDebug+1, Form("Found pythia process type %d", processType));
+
+  if (processType<0)
+    AliDebug(AliLog::kError, Form("Unknown Pythia process type %d.", processType));
+
+  // get the MC vertex
+  AliGenEventHeader* genHeader = header->GenEventHeader();
+  TArrayF vtxMC(3);
+  genHeader->PrimaryVertex(vtxMC);
+
+  // get the ESD vertex
+  const AliESDVertex* vtxESD = AliPWG0Helper::GetVertex(fESD, fAnalysisMode);
 
   Bool_t eventVertex = kFALSE;
-  if (AliPWG0Helper::GetVertex(fESD, fAnalysisMode))
+  Double_t vtx[3];
+  if (vtxESD) 
+  {
+    vtxESD->GetXYZ(vtx);
     eventVertex = kTRUE;
+    
+    Double_t diff = vtxMC[2] - vtx[2];
+    if (vtxESD->GetZRes() > 0) 
+      fVertexShiftNorm->Fill(diff / vtxESD->GetZRes());
+  } 
+  else
+    Printf("No vertex found");
 
-  // post the data already here
-  PostData(0, fOutput);
 
   // create list of (label, eta, pt) tuples
   Int_t inputCount = 0;
@@ -231,6 +295,8 @@ void AlidNdEtaCorrectionTask::Exec(Option_t*)
     TObjArray* list = fEsdTrackCuts->GetAcceptedTracks(fESD);
     Int_t nGoodTracks = list->GetEntries();
 
+    Printf("Accepted %d tracks", nGoodTracks);
+
     labelArr = new Int_t[nGoodTracks];
     etaArr = new Float_t[nGoodTracks];
     ptArr = new Float_t[nGoodTracks];
@@ -250,47 +316,18 @@ void AlidNdEtaCorrectionTask::Exec(Option_t*)
       ptArr[inputCount] = esdTrack->Pt();
       ++inputCount;
     }
-  }
-  else
-    return;
-
-  AliMCEventHandler* eventHandler = dynamic_cast<AliMCEventHandler*> (AliAnalysisManager::GetAnalysisManager()->GetMCtruthEventHandler());
-  if (!eventHandler) {
-    Printf("ERROR: Could not retrieve MC event handler");
-    return;
-  }
 
-  AliMCEvent* mcEvent = eventHandler->MCEvent();
-  if (!mcEvent) {
-    Printf("ERROR: Could not retrieve MC event");
-    return;
+    delete list;
   }
-
-  AliStack* stack = mcEvent->Stack();
-  if (!stack)
-  {
-    AliDebug(AliLog::kError, "Stack not available");
+  else
     return;
-  }
 
-  AliHeader* header = mcEvent->Header();
-  if (!header)
+  if (eventTriggered && eventVertex)
   {
-    AliDebug(AliLog::kError, "Header not available");
-    return;
+    fVertexCorrelation->Fill(vtxMC[2], vtx[2]);
+    fVertexProfile->Fill(vtxMC[2], vtxMC[2] - vtx[2]);
   }
 
-  // get process type; NB: this only works for Pythia
-  Int_t processType = AliPWG0Helper::GetPythiaEventProcessType(header);
-  AliDebug(AliLog::kDebug+1, Form("Found pythia process type %d", processType));
-
-  if (processType<0)
-    AliDebug(AliLog::kError, Form("Unknown Pythia process type %d.", processType));
-
-  // get the MC vertex
-  AliGenEventHeader* genHeader = header->GenEventHeader();
-  TArrayF vtxMC(3);
-  genHeader->PrimaryVertex(vtxMC);
 
   // loop over mc particles
   Int_t nPrim  = stack->GetNprimary();
@@ -340,7 +377,7 @@ void AlidNdEtaCorrectionTask::Exec(Option_t*)
 
     if (label < 0)
     {
-      AliDebug(AliLog::kWarning, Form("WARNING: cannot find corresponding mc part for track %d.", label));
+      Printf("WARNING: cannot find corresponding mc part for track(let) %d with label %d.", i, label);
       continue;
     }
 
@@ -504,10 +541,17 @@ void AlidNdEtaCorrectionTask::Terminate(Option_t *)
   if (fSigmaVertexPrim)
     fSigmaVertexPrim->Write();
 
+  if (fVertexCorrelation)
+    fVertexCorrelation->Write();
+  if (fVertexProfile)
+    fVertexProfile->Write();
+  if (fVertexShiftNorm)
+    fVertexShiftNorm->Write();
+
   fout->Write();
   fout->Close();
 
-  fdNdEtaCorrection->DrawHistograms();
+  //fdNdEtaCorrection->DrawHistograms();
 
   Printf("Writting result to %s", fileName.Data());
 
index c9419e4..a204ff3 100644 (file)
@@ -13,6 +13,8 @@ class AlidNdEtaCorrection;
 class TH1F;
 class AliESDEvent;
 class TParticlePDG;
+class TH2F;
+class TProfile;
 
 class AlidNdEtaCorrectionTask : public AliAnalysisTask {
   public:
@@ -47,6 +49,10 @@ class AlidNdEtaCorrectionTask : public AliAnalysisTask {
     // control histograms
     TH1F* fPIDParticles;                         //! pid of primary particles
     TH1F* fPIDTracks;                            //! pid of reconstructed tracks
+    TH2F* fVertexCorrelation;                    //! ESD z-vtx vs MC z-vtx
+    TProfile* fVertexProfile;                    //! Profile of MC z-vtx - ESD z-vtx vs. MC z-vtx
+    TH1F* fVertexShiftNorm;                      //! (MC z-vtx - ESD z-vtx) / (sigma_ESD-z-vtx) histogrammed
 
     // histograms for systematic studies (must be enabled with option)
 
index 917ea74..0480c61 100644 (file)
@@ -44,17 +44,20 @@ AlidNdEtaTask::AlidNdEtaTask(const char* opt) :
   fOption(opt),
   fAnalysisMode(AliPWG0Helper::kTPC),
   fReadMC(kFALSE),
+  fUseMCVertex(kFALSE),
   fEsdTrackCuts(0),
   fdNdEtaAnalysisESD(0),
   fMult(0),
   fMultVtx(0),
   fEvents(0),
+  fVertexResolution(0),
   fdNdEtaAnalysis(0),
   fdNdEtaAnalysisTr(0),
   fdNdEtaAnalysisTrVtx(0),
   fdNdEtaAnalysisTracks(0),
   fVertex(0),
-  fPartPt(0)
+  fPartPt(0),
+  fDeltaPhi(0)
 {
   //
   // Constructor. Initialization of pointers
@@ -99,12 +102,15 @@ void AlidNdEtaTask::ConnectInputData(Option_t *)
     tree->SetBranchStatus("fSPDVertex*", 1);
     // PrimaryVertex also needed
 
-    if (fAnalysisMode == AliPWG0Helper::kSPD)
+    if (fAnalysisMode == AliPWG0Helper::kSPD) {
       tree->SetBranchStatus("fSPDMult*", 1);
-
+      //AliPWG0Helper::SetBranchStatusRecursive(tree, "AliMultiplicity", kTRUE);
+    }
+    
     if (fAnalysisMode == AliPWG0Helper::kTPC || fAnalysisMode == AliPWG0Helper::kTPCITS) {
       AliESDtrackCuts::EnableNeededBranches(tree);
       tree->SetBranchStatus("fTracks.fLabel", 1);
+      //AliPWG0Helper::SetBranchStatusRecursive(tree, "Tracks", kTRUE);
     }
 
     AliESDInputHandler *esdH = dynamic_cast<AliESDInputHandler*> (AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler());
@@ -114,12 +120,17 @@ void AlidNdEtaTask::ConnectInputData(Option_t *)
     } else
       fESD = esdH->GetEvent();
   }
+
+  // disable info messages of AliMCEvent (per event)
+  AliLog::SetClassDebugLevel("AliMCEvent", AliLog::kWarning - AliLog::kDebug + 1);
 }
 
 void AlidNdEtaTask::CreateOutputObjects()
 {
   // create result objects and add to output list
 
+  Printf("AlidNdEtaTask::CreateOutputObjects");
+
   fOutput = new TList;
   fOutput->SetOwner();
 
@@ -142,6 +153,9 @@ void AlidNdEtaTask::CreateOutputObjects()
   fEvents = new TH1F("dndeta_check_vertex", "dndeta_check_vertex", 40, -20, 20);
   fOutput->Add(fEvents);
 
+  fVertexResolution = new TH1F("dndeta_vertex_resolution_z", "dndeta_vertex_resolution_z", 1000, 0, 10);
+  fOutput->Add(fVertexResolution);
+
   if (fReadMC)
   {
     fdNdEtaAnalysis = new dNdEtaAnalysis("dndeta", "dndeta", fAnalysisMode);
@@ -163,6 +177,9 @@ void AlidNdEtaTask::CreateOutputObjects()
     fPartPt->Sumw2();
     fOutput->Add(fPartPt);
   }
+
+  if (fAnalysisMode == AliPWG0Helper::kSPD)
+    fDeltaPhi = new TH1F("fDeltaPhi", "fDeltaPhi;#Delta #phi;Entries", 1000, -3.14, 3.14);
 }
 
 void AlidNdEtaTask::Exec(Option_t*)
@@ -183,8 +200,15 @@ void AlidNdEtaTask::Exec(Option_t*)
   const AliESDVertex* vtxESD = AliPWG0Helper::GetVertex(fESD, fAnalysisMode);
   
   Double_t vtx[3];
-  if (vtxESD)
+  if (vtxESD) {
     vtxESD->GetXYZ(vtx);
+  }
+  else
+    Printf("No vertex");
+  
+  // fill z vertex resolution
+  if (vtxESD)
+    fVertexResolution->Fill(vtxESD->GetZRes());
 
   // post the data already here
   PostData(0, fOutput);
@@ -217,11 +241,15 @@ void AlidNdEtaTask::Exec(Option_t*)
       if (mult->GetDeltaPhi(i) < -1000)
         continue;
 
+      fDeltaPhi->Fill(mult->GetDeltaPhi(i));
+
       etaArr[inputCount] = mult->GetEta(i);
       labelArr[inputCount] = mult->GetLabel(i);
       ptArr[inputCount] = 0; // no pt for tracklets
       ++inputCount;
     }
+
+    Printf("Accepted %d tracklets", inputCount);
   }
   else if (fAnalysisMode == AliPWG0Helper::kTPC || fAnalysisMode == AliPWG0Helper::kTPCITS)
   {
@@ -255,11 +283,47 @@ void AlidNdEtaTask::Exec(Option_t*)
       ++inputCount;
     }
 
+    Printf("Accepted %d tracks", nGoodTracks);
+
     delete list;
   }
   else
     return;
 
+  if (fUseMCVertex) {
+    Printf("WARNING: Replacing vertex by MC vertex. This is for systematical checks only.");
+    if (!fReadMC) {
+      Printf("ERROR: fUseMCVertex set without fReadMC set!");
+      return;
+    }
+
+    AliMCEventHandler* eventHandler = dynamic_cast<AliMCEventHandler*> (AliAnalysisManager::GetAnalysisManager()->GetMCtruthEventHandler());
+    if (!eventHandler) {
+      Printf("ERROR: Could not retrieve MC event handler");
+      return;
+    }
+
+    AliMCEvent* mcEvent = eventHandler->MCEvent();
+    if (!mcEvent) {
+      Printf("ERROR: Could not retrieve MC event");
+      return;
+    }
+
+    AliHeader* header = mcEvent->Header();
+    if (!header)
+    {
+      AliDebug(AliLog::kError, "Header not available");
+      return;
+    }
+
+    // get the MC vertex
+    AliGenEventHeader* genHeader = header->GenEventHeader();
+    TArrayF vtxMC(3);
+    genHeader->PrimaryVertex(vtxMC);
+
+    vtx[2] = vtxMC[2];
+  }
+
   // Processing of ESD information (always)
   if (eventTriggered)
   {
@@ -291,7 +355,7 @@ void AlidNdEtaTask::Exec(Option_t*)
       // for event count per vertex
       fdNdEtaAnalysisESD->FillEvent(vtx[2], inputCount);
 
-      // control hists
+      // control hist
       fEvents->Fill(vtx[2]);
     }
   }
@@ -349,8 +413,8 @@ void AlidNdEtaTask::Exec(Option_t*)
       Float_t eta = particle->Eta();
       Float_t pt = particle->Pt();
 
-       // make a rough eta cut (so that nAcceptedParticles is not too far off)
-      if (TMath::Abs(eta) < 1.5)
+       // make a rough eta cut (so that nAcceptedParticles is not too far off the true measured value (NB: this histograms are only gathered for comparison))
+      if (TMath::Abs(eta) < 1.5 && pt > 0.3)
         nAcceptedParticles++;
 
       fdNdEtaAnalysis->FillTrack(vtxMC[2], eta, pt);
@@ -424,6 +488,7 @@ void AlidNdEtaTask::Terminate(Option_t *)
   fMult = dynamic_cast<TH1F*> (fOutput->FindObject("fMult"));
   fMultVtx = dynamic_cast<TH1F*> (fOutput->FindObject("fMultVtx"));
   fEvents = dynamic_cast<TH1F*> (fOutput->FindObject("dndeta_check_vertex"));
+  fVertexResolution = dynamic_cast<TH1F*> (fOutput->FindObject("dndeta_vertex_resolution_z"));
   for (Int_t i=0; i<3; ++i)
     fPartEta[i] = dynamic_cast<TH1F*> (fOutput->FindObject(Form("dndeta_check_%d", i)));
 
@@ -446,6 +511,10 @@ void AlidNdEtaTask::Terminate(Option_t *)
     Int_t events1 = (Int_t) fEvents->Integral(fEvents->GetXaxis()->FindBin(-19.9), fEvents->GetXaxis()->FindBin(-0.001));
     Int_t events2 = (Int_t) fEvents->Integral(fEvents->GetXaxis()->FindBin(0.001), fEvents->GetXaxis()->FindBin(19.9));
 
+    Printf("%d events with vertex used", events1 + events2);
+
+    if (events1 > 0 && events2 > 0)
+    {
     fPartEta[0]->Scale(1.0 / (events1 + events2));
     fPartEta[1]->Scale(1.0 / events1);
     fPartEta[2]->Scale(1.0 / events2);
@@ -459,6 +528,7 @@ void AlidNdEtaTask::Terminate(Option_t *)
       fPartEta[i]->SetLineColor(i+1);
       fPartEta[i]->Draw((i==0) ? "" : "SAME");
     }
+    }
   }
 
   if (fEvents)
@@ -488,6 +558,12 @@ void AlidNdEtaTask::Terminate(Option_t *)
   if (fEvents)
     fEvents->Write();
 
+  if (fVertexResolution)
+    fVertexResolution->Write();
+
+  if (fDeltaPhi)
+    fDeltaPhi->Write();
+
   fout->Write();
   fout->Close();
 
index 4fb6bef..66a4073 100644 (file)
@@ -26,6 +26,7 @@ class AlidNdEtaTask : public AliAnalysisTask {
     void SetTrackCuts(AliESDtrackCuts* cuts) { fEsdTrackCuts = cuts; }
     void SetAnalysisMode(AliPWG0Helper::AnalysisMode mode) { fAnalysisMode = mode; }
     void SetReadMC(Bool_t flag = kTRUE) { fReadMC = flag; }
+    void SetUseMCVertex(Bool_t flag = kTRUE) { fUseMCVertex = flag; }
 
  protected:
     AliESDEvent *fESD;    //! ESD object
@@ -35,6 +36,7 @@ class AlidNdEtaTask : public AliAnalysisTask {
     AliPWG0Helper::AnalysisMode fAnalysisMode; // detector that is used for analysis
 
     Bool_t  fReadMC;       // if true reads MC data (to build correlation maps)
+    Bool_t  fUseMCVertex;  // the MC vtx is used instead of the ESD vertex (for syst. check)
 
     AliESDtrackCuts* fEsdTrackCuts;         // Object containing the parameters of the esd track cuts
 
@@ -45,6 +47,7 @@ class AlidNdEtaTask : public AliAnalysisTask {
     TH1F* fMultVtx;                            //! raw multiplicity histogram of evts with vtx (control histogram)
     TH1F* fPartEta[3];            //! counted particles as function of eta (full vertex range, below 0 range, above 0 range)
     TH1F* fEvents;                //! events counted as function of vtx
+    TH1F* fVertexResolution;      //! z resolution of the vertex
 
     // Gathered from MC (when fReadMC is set)
     dNdEtaAnalysis* fdNdEtaAnalysis;        //! contains the dndeta from the full sample
@@ -54,6 +57,7 @@ class AlidNdEtaTask : public AliAnalysisTask {
     // the following are control histograms to check the dNdEtaAnalysis class
     TH3F* fVertex;                //! vertex of counted particles
     TH1F* fPartPt;                //! counted particles as function of pt
+    TH1F* fDeltaPhi;              //! histogram of delta_phi values for tracklets (only for SPD analysis)
 
  private:
     AlidNdEtaTask(const AlidNdEtaTask&);