Added central dN/deta task based on tracks
authorcholm <cholm@f7af4fe6-9843-0410-8265-dc069ae4e863>
Fri, 28 Jan 2011 13:27:27 +0000 (13:27 +0000)
committercholm <cholm@f7af4fe6-9843-0410-8265-dc069ae4e863>
Fri, 28 Jan 2011 13:27:27 +0000 (13:27 +0000)
PWG2/FORWARD/analysis2/AddTaskCentralMult.C [new file with mode: 0644]
PWG2/FORWARD/analysis2/AnalyseAOD.C
PWG2/FORWARD/analysis2/MakeAOD.C
PWG2/FORWARD/analysis2/scripts/AliESDCentrality.C [new file with mode: 0644]

diff --git a/PWG2/FORWARD/analysis2/AddTaskCentralMult.C b/PWG2/FORWARD/analysis2/AddTaskCentralMult.C
new file mode 100644 (file)
index 0000000..d06320a
--- /dev/null
@@ -0,0 +1,621 @@
+/**
+ * @file   AddCentralMultTask.C
+ * @author Christian Holm Christensen <cholm@dalsgaard.hehi.nbi.dk>
+ * @date   Fri Jan 28 10:22:26 2011
+ * 
+ * @brief Class and script to add a multiplicity task for the central
+ *        @f$\eta@f$ region
+ * 
+ * 
+ */
+#include <AliAnalysisTaskSE.h>
+#include <AliESDtrackCuts.h>
+class TList;
+class TH2D;
+class TH1D;
+
+/**
+ * Task to determine the 
+ * @f[
+ *   \left.\frac{d^2N_{ch}}{d\eta d\phi}\right|_{central}
+ * @f] 
+ * from tracks and tracklets. 
+ *
+ * First, global tracks are investigated. The requirements on global
+ * tracks are
+ * - @f$ n_{TPC clusters} \ge 70@f$ 
+ * - @f$ \chi^2_{TPC cluster} \le 4@f$ 
+ * - No daughter kinks 
+ * - Re-fit of TPC tracks 
+ * - Re-fit of ITS tracks 
+ * - No requirement on SPD clusters 
+ *
+ * Secondly, ITS stand-alone tracks are investigated.  The
+ * requirements on ITS standalone tracks are
+ * - Re-fit of ITS tracks 
+ * - No requirement on SPD clusters 
+ * 
+ * Tracks that does not meet these quality conditions are flagged as
+ * rejected.
+ *
+ * Both kinds of tracks (global and ITS standalone) have requirements
+ * on the distance of closest approach (DCA) to the interaction
+ * vertex @f$(v_x,v_y,v_z)@f$. 
+ *
+ * For tracks with SPD clusters: 
+ * - @f$ DCA_{xy} < 0.0182+0.0350/p_t^{1.01}@f$ 
+ * - @f$ DCA_{z} < 0.5@f$ 
+ *
+ * For tracks without SPD clusters 
+ * - @f$ DCA_{xy} < 1.5(0.0182+0.0350/p_t^{1.01})@f$ 
+ * - @f$ DCA_{z} < 0.5@f$ 
+ *
+ * Tracks that does not meet these DCA requirements are flagged as
+ * secondaries.
+ *
+ * Thirdly, the number of SPD tracklets are investigated.  If the
+ * tracklet is associated with a track, and that track has already
+ * been used, then that tracklet is ignored
+ * 
+ * An @f$(\eta,\phi)@f$ per-event histogram is then filled from these
+ * tracks and tracklets, and that histogram is stored in the output
+ * AOD event.
+ *
+ * Furthermore, a number of diagnostics @f$\eta@f$ histograms are filled: 
+ * - @f$\eta@f$ of all accepted tracks and tracklets 
+ * - @f$\eta@f$ of accepted global tracks
+ * - @f$\eta@f$ of accepted ITS tracks
+ * - @f$\eta@f$ of accepted SPD tracklets
+ *
+ * At the end of the job, these histograms are normalized to the
+ * number of accepted events and bin width to provide an estimate of
+ * the @f$dN_{ch}/d\eta@f$
+ *
+ * Only minimum bias events with a @f$v_z@f$ within the defined cut
+ * are analysed.
+ */
+class CentralMultTask : public AliAnalysisTaskSE
+{
+public:
+  enum { 
+    kValidTrigger = 0x1, 
+    kValidVertex  = 0x2
+  };
+  /** 
+   * Constructor 
+   * 
+   */
+  CentralMultTask();
+  /** 
+   * Constructor
+   * 
+   * @param name    Name of task 
+   * @param maxEta  Set @f$\eta@f$ range
+   * @param maxVtx  Set @f$v_z@f$ range
+   */
+  CentralMultTask(const char* name, 
+                 Double_t maxEta=2, 
+                 Double_t maxVtx=10);
+  /**
+   * Destructor
+   * 
+   */
+  virtual ~CentralMultTask();
+  /** 
+   * Initialise on master - does nothing
+   * 
+   */
+  virtual void   Init() {}
+  /** 
+   * Create output objects.  
+   *
+   * This is called once per slave process 
+   */
+  virtual void UserCreateOutputObjects();
+  /** 
+   * Process a single event 
+   * 
+   * @param option Not used
+   */
+  virtual void UserExec(Option_t* option);
+  /** 
+   * Called at end of event processing.. 
+   *
+   * This is called once in the master 
+   * 
+   * @param option Not used 
+   */
+  virtual void Terminate(Option_t* option);
+protected:
+  /** 
+   * Check if this event is within cuts
+   * 
+   * @param esd Event
+   * 
+   * @return true if this event is to be considered
+   */
+  UShort_t CheckEvent(const AliESDEvent& esd, Double_t& vz);
+
+  TH2D*           fHist;      // Per-event d2n/deta/dphi
+  TH1D*           fAll;       // Summed d2n/deta/dphi - all track(let)s
+  TH1D*           fGlobal;    // Summed d2n/deta/dphi - global tracks
+  TH1D*           fITS;       // Summed d2n/deta/dphi - ITS tracks
+  TH1D*           fTracklets; // Summed d2n/deta/dphi - SPD tracklets
+  TH1D*           fNEventsTr; // Number of triggered events per vertex bin
+  TH1D*           fNEventsVtx;// Number of triggered+vertex events per vertex
+  TList*          fOutput;    // Output list
+  AliESDtrackCuts fQGlo;      // Quality cut on ITS+TPC
+  AliESDtrackCuts fQITS;      // Quality cut on ITS standalone (not ITS+TPC)
+  AliESDtrackCuts fDCAwSPD;   // DCA for traks with SPD hits
+  AliESDtrackCuts fDCAwoSPD;  // DCA for traks without SPD hits
+  AliESDtrackCuts fIsPrimary; // Primary particle 
+  Bool_t          fDebug;     // Debug flag
+
+  ClassDef(CentralMultTask,1); // Determine multiplicity in central area
+};
+
+//====================================================================
+#include <TMath.h>
+#include <TH2D.h>
+#include <TH1D.h>
+#include <THStack.h>
+#include <TList.h>
+#include <AliAnalysisManager.h>
+#include <AliESDEvent.h>
+#include <AliAODHandler.h>
+#include <AliAODEvent.h>
+#include <AliESDInputHandler.h>
+#include <AliESDtrack.h>
+#include <AliMultiplicity.h>
+
+//____________________________________________________________________
+inline CentralMultTask::CentralMultTask()
+  : AliAnalysisTaskSE(), 
+    fHist(0),
+    fAll(0),
+    fGlobal(0),
+    fITS(0),
+    fTracklets(0),
+    fNEventsTr(0),
+    fNEventsVtx(0),
+    fOutput(0),
+    fDebug(false)
+{}
+
+//____________________________________________________________________
+inline CentralMultTask::CentralMultTask(const char* name,
+                                Double_t maxEta,
+                                Double_t maxVtx)
+  : AliAnalysisTaskSE(name), 
+    fHist(0),
+    fAll(0),
+    fGlobal(0),
+    fITS(0),
+    fTracklets(0),
+    fNEventsTr(0),
+    fNEventsVtx(0),
+    fOutput(0),
+    fDebug(false)
+{
+  fHist = new TH2D("Central", "d^{2}N_{ch}/d#etad#phi in central region",
+                  80, -maxEta, maxEta, 20, 0, 2*TMath::Pi());
+  fHist->SetDirectory(0);
+  fHist->SetXTitle("#eta");
+  fHist->SetYTitle("#phi [radians]");
+  fHist->SetZTitle("d^{2}N_{ch}/d#etad#phi");
+  fHist->SetStats(0);
+  fHist->Sumw2();
+
+  fAll = new TH1D("all", "Central region", 80, -maxEta, maxEta);
+  fAll->SetDirectory(0);
+  fAll->SetXTitle("#eta");
+  fAll->SetYTitle("dN_{ch}/d#eta");
+  fAll->Sumw2();
+  fAll->SetFillColor(kGray);
+  fAll->SetFillStyle(3001);
+  fAll->SetMarkerStyle(28);
+  fAll->SetMarkerColor(kGray);
+  fAll->SetStats(0);
+  
+  fGlobal = static_cast<TH1D*>(fAll->Clone("global"));
+  fGlobal->SetTitle("Global tracks");
+  fGlobal->SetDirectory(0);
+  fGlobal->Sumw2();
+  fGlobal->SetFillColor(kRed+1);
+  fGlobal->SetFillStyle(3001);
+  fGlobal->SetMarkerStyle(28);
+  fGlobal->SetMarkerColor(kRed+1);
+  fGlobal->SetStats(0);
+
+  fITS = static_cast<TH1D*>(fAll->Clone("its"));
+  fITS->SetTitle("ITS tracks");
+  fITS->SetDirectory(0);
+  fITS->Sumw2();
+  fITS->SetFillColor(kGreen+1);
+  fITS->SetFillStyle(3001);
+  fITS->SetMarkerStyle(28);
+  fITS->SetMarkerColor(kGreen+1);
+  fITS->SetStats(0);
+
+  fTracklets = static_cast<TH1D*>(fAll->Clone("tracklets"));
+  fTracklets->SetTitle("SPD tracklets");
+  fTracklets->SetDirectory(0);
+  fTracklets->Sumw2();
+  fTracklets->SetFillColor(kBlue+1);
+  fTracklets->SetFillStyle(3001);
+  fTracklets->SetMarkerStyle(28);
+  fTracklets->SetMarkerColor(kBlue+1);
+  fTracklets->SetStats(0);
+
+  fNEventsTr = new TH1D("nEventsTr", 
+                       "Events per vertex bin", 10, -maxVtx, maxVtx);
+  fNEventsTr->SetXTitle("v_{z}");
+  fNEventsTr->SetYTitle("Events");
+  fNEventsTr->SetDirectory(0);
+
+  fNEventsVtx = new TH1D("nEventsVtx", 
+                        "Events per vertex bin", 10, -maxVtx, maxVtx);
+  fNEventsVtx->SetXTitle("v_{z}");
+  fNEventsVtx->SetYTitle("Events");
+  fNEventsVtx->SetDirectory(0);
+
+  // --- Global (ITS+TPC) track quality cuts 
+  // TPC  
+  fQGlo.SetMinNClustersTPC(70);
+  fQGlo.SetMaxChi2PerClusterTPC(4);
+  fQGlo.SetAcceptKinkDaughters(kFALSE);
+  fQGlo.SetRequireTPCRefit(kTRUE);
+  // ITS
+  fQGlo.SetRequireITSRefit(kTRUE);
+  fQGlo.SetClusterRequirementITS(AliESDtrackCuts::kSPD,AliESDtrackCuts::kOff);
+  fQGlo.SetEtaRange(-maxEta, maxEta);
+
+  // --- ITS standalone track quality cuts 
+  fQITS.SetRequireITSRefit(kTRUE);
+  fQITS.SetClusterRequirementITS(AliESDtrackCuts::kSPD,AliESDtrackCuts::kOff);
+  fQITS.SetEtaRange(-maxEta, maxEta); 
+
+  // -- Distance-of-Closest-Approach cuts for tracks w/ITS hits 
+  fDCAwSPD.SetClusterRequirementITS(AliESDtrackCuts::kSPD,
+                                   AliESDtrackCuts::kAny);
+  fDCAwSPD.SetMaxDCAToVertexXYPtDep("0.0182+0.0350/pt^1.01");
+  fDCAwSPD.SetMaxDCAToVertexZ(0.5);
+  fDCAwSPD.SetEtaRange(-maxEta, maxEta);
+
+  // -- Distance-of-Closest-Approach cuts for tracks w/o ITS hits 
+  fDCAwoSPD.SetClusterRequirementITS(AliESDtrackCuts::kSPD,
+                                    AliESDtrackCuts::kNone);
+  fDCAwoSPD.SetMaxDCAToVertexXYPtDep("1.5*(0.0182+0.0350/pt^1.01)");
+  fDCAwoSPD.SetMaxDCAToVertexZ(0.5);
+  fDCAwoSPD.SetEtaRange(-maxEta, maxEta);  
+
+  // -- Primary track cut 
+  // https://twiki.cern.ch/twiki/bin/view/ALICE/SelectionOfPrimaryTracksForPpDataAnalysis
+  // Quality
+  fIsPrimary.SetMinNClustersTPC(70);
+  fIsPrimary.SetMaxChi2PerClusterTPC(4);
+  fIsPrimary.SetAcceptKinkDaughters(kFALSE);
+  fIsPrimary.SetRequireTPCRefit(kTRUE);
+  fIsPrimary.SetRequireITSRefit(kTRUE);
+  fIsPrimary.SetClusterRequirementITS(AliESDtrackCuts::kSPD,
+                                     AliESDtrackCuts::kAny);
+  //  Dca:
+  fIsPrimary.SetMaxDCAToVertexXYPtDep("0.0182+0.0350/pt^1.01");
+  fIsPrimary.SetMaxDCAToVertexZ(2);
+  fIsPrimary.SetDCAToVertex2D(kFALSE);
+  fIsPrimary.SetRequireSigmaToVertex(kFALSE);  
+
+  // Output slot #1 writes into a TH1 container
+  DefineOutput(1, TList::Class()); 
+}
+
+//____________________________________________________________________
+inline CentralMultTask::~CentralMultTask()
+{
+  if (fHist)      delete fHist;
+  if (fAll)       delete fAll;
+  if (fGlobal)    delete fGlobal;
+  if (fITS)       delete fITS;
+  if (fTracklets) delete fTracklets;
+}
+
+//________________________________________________________________________
+inline void 
+CentralMultTask::UserCreateOutputObjects()
+{
+  // Create histograms
+  // Called once (on the worker node)
+
+  fOutput = new TList;
+  fOutput->SetName(GetName());
+  fOutput->SetOwner();
+
+  fOutput->Add(fAll);
+  fOutput->Add(fGlobal);
+  fOutput->Add(fITS);
+  fOutput->Add(fTracklets);
+  fOutput->Add(fNEventsTr);
+  fOutput->Add(fNEventsVtx);
+
+  AliAnalysisManager* am = AliAnalysisManager::GetAnalysisManager();
+  AliAODHandler*      ah = 
+    dynamic_cast<AliAODHandler*>(am->GetOutputEventHandler());
+  if (!ah) AliFatal("No AOD output handler set in analysis manager");
+
+  ah->AddBranch("TH2D", &fHist);
+  
+  // Post data for ALL output slots >0 here, to get at least an empty histogram
+  PostData(1, fOutput); 
+}
+
+//____________________________________________________________________
+inline UShort_t
+CentralMultTask::CheckEvent(const AliESDEvent& esd, Double_t& vz)
+{
+  // Do some fast cuts first
+  vz = 0;
+
+  // Get the analysis manager - should always be there 
+  AliAnalysisManager* am = AliAnalysisManager::GetAnalysisManager();
+  if (!am) { 
+    AliWarning("No analysis manager defined!");
+    return 0;
+  }
+
+  // Get the input handler - should always be there 
+  AliInputEventHandler* ih = 
+    static_cast<AliInputEventHandler*>(am->GetInputEventHandler());
+  if (!ih) { 
+    AliWarning("No input handler");
+    return 0;
+  }
+
+  // Trigger mask 
+  UInt_t    mask     = ih->IsEventSelected();   
+  Bool_t   isMinBias = (mask == AliVEvent::kMB) ? 1 : 0;
+  UShort_t ret       = 0;
+  if (isMinBias) ret |= kValidTrigger;
+  
+  // check for good reconstructed vertex
+  if (!(esd.GetPrimaryVertex()->GetStatus())) {
+    if (fDebug)
+      AliWarning("Primary vertex has bad status");
+    return ret;
+  }
+  if (!(esd.GetPrimaryVertexSPD()->GetStatus())) { 
+    if (fDebug)
+      AliWarning("Primary SPD vertex has bad status");
+    return ret;
+  }
+
+  // if vertex is from spd vertexZ, require more stringent cut
+  if (esd.GetPrimaryVertex()->IsFromVertexerZ()) {
+    if (esd.GetPrimaryVertex()->GetDispersion() > 0.02 ||  
+       esd.GetPrimaryVertex()->GetZRes()       > 0.25) {
+      if (fDebug)
+       AliWarning(Form("Primary vertex dispersion=%f (0.02) zres=%f (0.05)", 
+                       esd.GetPrimaryVertex()->GetDispersion(),
+                       esd.GetPrimaryVertex()->GetZRes()));
+      return ret;
+    }
+  }
+  // One can add a cut on the vertex z position here
+  vz = esd.GetPrimaryVertex()->GetZ();
+  Double_t vl = fNEventsVtx->GetXaxis()->GetXmin();
+  Double_t vh = fNEventsVtx->GetXaxis()->GetXmax();
+  if (vz <  vl || vz > vh) {
+    if (fDebug)
+      AliWarning(Form("Primary vertex vz=%f out of range [%f,%f]", vz, vl, vh));
+    return ret; 
+  }
+  ret |= kValidVertex;
+
+  return ret;
+}
+
+//____________________________________________________________________
+inline void 
+CentralMultTask::UserExec(Option_t *) 
+{
+  // Main loop
+  AliESDEvent* esd = dynamic_cast<AliESDEvent*>(InputEvent());
+  if (!esd) {
+    AliError("Cannot get the ESD event");
+    return;
+  }  
+
+  fHist->Reset();
+
+  // Check event 
+  Double_t vz         = 0;
+  UShort_t eventFlags = CheckEvent(*esd, vz);
+  if (!(eventFlags & kValidTrigger)) 
+    return; // No trigger
+  else {
+    fNEventsTr->Fill(vz);
+    if (eventFlags & kValidVertex) 
+      fNEventsVtx->Fill(vz);
+    else 
+      return; // No trigger or no vertex
+  }
+  
+
+  // flags for secondary and rejected tracks
+  // set this bit in ESD tracks if it is rejected by a cut
+  const int kRejBit = BIT(15); 
+  // set this bit in ESD tracks if it is secondary according to a cut
+  const int kSecBit = BIT(16); 
+
+  Int_t total      = 0;
+  Int_t nESDTracks = esd->GetNumberOfTracks();
+  // Loop over ESD tracks 
+  for(Int_t i = 0; i < nESDTracks; i++){ // flag the tracks
+
+    AliESDtrack* track = esd->GetTrack(i);
+    
+    // if track is a secondary from a V0, flag as a secondary
+    if (track->IsOn(AliESDtrack::kMultInV0)) {
+      track->SetBit(kSecBit); 
+      continue;
+    } 
+
+    Double_t eta = track->Eta();
+    Double_t phi = track->Phi();
+    // check tracks with ITS part
+    if (track->IsOn(AliESDtrack::kITSin) && 
+       !track->IsOn(AliESDtrack::kITSpureSA)) { 
+      // track has ITS part but is not an ITS_SA -> TPC+ITS
+      if (track->IsOn(AliESDtrack::kTPCin)) {  
+       // Global track, has ITS and TPC contributions
+       if (fQGlo.AcceptTrack(track)) { // good TPCITS track
+         if (fDCAwSPD.AcceptTrack(track) || 
+             fDCAwoSPD.AcceptTrack(track)) {
+           fGlobal->Fill(eta);
+           fAll->Fill(eta);
+           total++;
+         }
+         else 
+           // large DCA -> secondary, don't count either track not
+           // associated tracklet
+           track->SetBit(kSecBit); 
+       }
+       else 
+         // bad quality, don't count the track, but may count
+         // tracklet if associated
+         track->SetBit(kRejBit); 
+      } // if (track->IsOn(AliESDtrack::kTPCin))
+      // ITS complimentary
+      else if (fQITS.AcceptTrack(track)) { 
+       // good ITS complimentary track
+       if (fDCAwSPD.AcceptTrack(track) || 
+           fDCAwoSPD.AcceptTrack(track)) {
+         fITS->Fill(eta);
+         fAll->Fill(eta);
+         total++;
+       }
+       else 
+         // large DCA -> secondary, don't count either track not
+         // associated tracklet
+         track->SetBit(kSecBit); 
+      } // else if (fQITS.AcceptTrack(track))
+      else 
+       // bad quality, don't count the track, but may count tracklet
+       // if associated
+       track->SetBit(kRejBit); 
+      if (track->IsOn(kSecBit) || track->IsOn(kRejBit)) continue;
+
+      // fill d2n/detadphi
+      fHist->Fill(eta, phi);
+    } /* track->IsOn(AliESDtrack::kITSin) && 
+       * !track->IsOn(AliESDtrack::kITSpureSA) */
+  }
+  
+  // Get multiplicity from SPD tracklets 
+  const AliMultiplicity* spdmult = esd->GetMultiplicity();
+  for (Int_t i=0; i<spdmult->GetNumberOfTracklets(); ++i){
+    // if counting tracks+tracklets, 
+    // check if clusters were already used in tracks
+    Int_t id1,id2;
+    spdmult->GetTrackletTrackIDs(i,0,id1,id2);
+    AliESDtrack* tr1 = id1>=0 ? esd->GetTrack(id1) : 0;
+    AliESDtrack* tr2 = id2>=0 ? esd->GetTrack(id2) : 0;
+    //
+    if ((tr1 && tr1->TestBit(kSecBit)) || 
+       (tr2 && tr2->TestBit(kSecBit))) 
+      // reject as secondary
+      continue; 
+    if ((tr1 && !tr1->TestBit(kRejBit)) || 
+       (tr2 && !tr2->TestBit(kRejBit))) 
+      // already accounted for
+      continue; 
+    ++total;
+    Double_t eta = spdmult->GetEta(i);
+    Double_t phi = spdmult->GetPhi(i);
+    
+    // Incremetn d2n/detadphi from an SPD tracklet 
+    fHist->Fill(eta, phi);
+    fTracklets->Fill(eta);
+    fAll->Fill(eta);
+  }
+  if (fDebug)
+    AliInfo(Form("A total of %d tracks", total));
+
+  // NEW HISTO should be filled before this point, as PostData puts the
+  // information for this iteration of the UserExec in the container
+  PostData(1, fOutput);
+}
+
+  
+//________________________________________________________________________
+inline void 
+CentralMultTask::Terminate(Option_t *) 
+{
+  // Draw result to screen, or perform fitting, normalizations Called
+  // once at the end of the query
+        
+  fOutput = dynamic_cast<TList*> (GetOutputData(1));
+  if(!fOutput) {
+    AliError("Could not retrieve TList fOutput"); 
+    return; 
+  }
+
+  TH1D* all       = static_cast<TH1D*>(fOutput->FindObject("all"));
+  TH1D* global    = static_cast<TH1D*>(fOutput->FindObject("global"));
+  TH1D* its       = static_cast<TH1D*>(fOutput->FindObject("its"));
+  TH1D* tracklets = static_cast<TH1D*>(fOutput->FindObject("tracklets"));
+  TH1D* eventsTr  = static_cast<TH1D*>(fOutput->FindObject("nEventsTr"));
+  TH1D* eventsVtx = static_cast<TH1D*>(fOutput->FindObject("nEventsVtx"));
+
+  Int_t nTriggers  = eventsTr->GetEntries();
+  Int_t nVertex    = eventsVtx->GetEntries();
+  if (nTriggers <= 0 || nVertex <= 0) {
+    AliWarning("No data in the events histogram!");
+    nTriggers = 1;
+  }
+
+  all      ->Scale(1. / nTriggers, "width");
+  global   ->Scale(1. / nTriggers, "width");
+  its      ->Scale(1. / nTriggers, "width");
+  tracklets->Scale(1. / nTriggers, "width");
+
+  THStack* stack = new THStack("components", "Components");
+  stack->Add(tracklets);
+  stack->Add(global);
+  stack->Add(its);
+
+  fOutput->Add(stack);
+}
+
+//========================================================================
+inline AliAnalysisTask*
+AddTaskCentralMult()
+{
+  // analysis manager
+  AliAnalysisManager* mgr = AliAnalysisManager::GetAnalysisManager();
+  
+  // Make our object.  2nd argumenent is absolute max Eta 
+  // 3rd argument is absolute max Vz
+  CentralMultTask* task = new CentralMultTask("Global", 2, 10);
+  // if physics selection performed in UserExec(),
+  // this line should be commented 
+  // task->SelectCollisionCandidates(AliVEvent::kMB); 
+  mgr->AddTask(task);
+
+  // create containers for input/output
+  AliAnalysisDataContainer *output = 
+    mgr->CreateContainer("Central", TList::Class(), 
+                        AliAnalysisManager::kOutputContainer, 
+                        AliAnalysisManager::GetCommonFileName());
+  
+  // connect input/output
+  mgr->ConnectInput(task, 0, mgr->GetCommonInputContainer());
+  mgr->ConnectOutput(task, 1, output);
+
+  return task;
+}
+
+  
+//________________________________________________________________________
+//
+// EOF
+// 
index 4a066f365c0e1fe98d8ebca48554faf87e1de12c..fb40b177aea012799ad2705e79d364d260c3f35a 100644 (file)
@@ -62,6 +62,10 @@ public:
   TH2D*              fPrimary;
   /** Sum primary information */
   TH2D*              fSumPrimary;
+  /** Central region */
+  TH2D*              fCentral;
+  /** Central region */
+  TH2D*              fSumCentral;
   /** Vertex efficiency */
   Double_t           fVtxEff;
   /** Title to put on the plot */
@@ -96,6 +100,8 @@ public:
       fMCSum(0),
       fPrimary(0), 
       fSumPrimary(0),
+      fCentral(0),
+      fSumCentral(0),
       fVtxEff(0),
       fTitle(""),
       fDoHHD(kTRUE)
@@ -118,11 +124,16 @@ public:
     if (fSum)        delete fSum;
     if (fMCSum)      delete fMCSum;
     if (fSumPrimary) delete fSumPrimary;
+    if (fCentral)    delete fCentral;
+    if (fSumCentral) delete fSumCentral;
     fTree       = 0;
     fOut        = 0;
     fSum        = 0;
     fMCSum      = 0;
+    fPrimary    = 0;
     fSumPrimary = 0;
+    fCentral    = 0;
+    fSumCentral = 0;
     fVtxEff     = 0;
     
   }
@@ -282,6 +293,10 @@ public:
     // Set the branch pointer 
     if (fTree->GetBranch("primary"))
       fTree->SetBranchAddress("primary", &fPrimary);
+    
+    // Set branch pointer 
+    if (fTree->GetBranch("Central")) 
+      fTree->SetBranchAddress("Central", &fCentral);
 
     fOut = TFile::Open(outname, "RECREATE");
     if (!fOut) { 
@@ -391,12 +406,23 @@ public:
        fSumPrimary = 
          static_cast<TH2D*>(fPrimary->Clone("primarySum"));
        Info("Process", "Created MC truth histogram (%d,%f,%f)x(%d,%f,%f)", 
-            fMCSum->GetNbinsX(), 
-            fMCSum->GetXaxis()->GetXmin(),
-            fMCSum->GetXaxis()->GetXmax(),
-            fMCSum->GetNbinsY(), 
-            fMCSum->GetYaxis()->GetXmin(),
-            fMCSum->GetYaxis()->GetXmax());
+            fSumPrimary->GetNbinsX(), 
+            fSumPrimary->GetXaxis()->GetXmin(),
+            fSumPrimary->GetXaxis()->GetXmax(),
+            fSumPrimary->GetNbinsY(), 
+            fSumPrimary->GetYaxis()->GetXmin(),
+            fSumPrimary->GetYaxis()->GetXmax());
+      }
+      if (!fSumCentral && fTree->GetBranch("Central")) { 
+       fSumCentral = 
+         static_cast<TH2D*>(fCentral->Clone("centralSum"));
+       Info("Process", "Created Cental histogram (%d,%f,%f)x(%d,%f,%f)", 
+            fSumCentral->GetNbinsX(), 
+            fSumCentral->GetXaxis()->GetXmin(),
+            fSumCentral->GetXaxis()->GetXmax(),
+            fSumCentral->GetNbinsY(), 
+            fSumCentral->GetYaxis()->GetXmin(),
+            fSumCentral->GetYaxis()->GetXmax());
       }
       
       // Add contribution from this event 
@@ -423,7 +449,8 @@ public:
       // Add contribution from this event
       if (fMCSum) fMCSum->Add(&(fMCAOD->GetHistogram()));
 
-      
+      // Add contribution from this event
+      if (fSumCentral) fSumCentral->Add(fCentral);      
     }
     printf("\n");
     fVtxEff = Double_t(fNWithVertex)/fNTriggered;
@@ -497,15 +524,33 @@ public:
       dndetaTruth->SetFillStyle(0);
       Rebin(dndetaTruth, rebin);
     }
+    TH1D* dndetaCentral = 0;
+    if (fSumCentral) { 
+      dndetaCentral = fSumCentral->ProjectionX("dndetaCentral", -1, -1, "e");
+      dndetaCentral->SetTitle("ALICE Central");
+      dndetaCentral->Scale(1./fNTriggered, "width");
+      dndetaCentral->SetMarkerColor(kGray+3);
+      dndetaCentral->SetMarkerStyle(22);
+      dndetaCentral->SetMarkerSize(1);
+      dndetaCentral->SetFillStyle(0);
+      Rebin(dndetaCentral, rebin <= 1 ? 1 : 2*(rebin/2));
+      // 1 -> 1
+      // 2 -> 2*2/2 -> 2*1 -> 2
+      // 3 -> 2*3/2 -> 2*1 -> 2
+      // 4 -> 2*4/2 -> 2*2 -> 4
+      // 5 -> 2*5/2 -> 2*2 -> 4
+      // 6 -> 2*6/2 -> 2*3 -> 6
+    }
 
-    DrawIt(dndeta, dndetaMC, dndetaTruth, mask, energy, doHHD, doComp);
+    DrawIt(dndeta, dndetaMC, dndetaTruth, dndetaCentral, 
+          mask, energy, doHHD, doComp);
 
     return kTRUE;
   }
   //__________________________________________________________________
   /** 
    */
-  void DrawIt(TH1* dndeta, TH1* dndetaMC, TH1* dndetaTruth, 
+  void DrawIt(TH1* dndeta, TH1* dndetaMC, TH1* dndetaTruth, TH1* dndetaCentral,
              Int_t mask, Int_t energy, Bool_t doHHD, 
              Bool_t doComp)
   {
@@ -563,6 +608,13 @@ public:
     stack->Add(dndetaSym);
     stack->Add(dndeta);
 
+    // If we have central region data, add that 
+    if (dndetaCentral) { 
+      stack->Add(dndetaCentral);
+      max = TMath::Max(dndetaCentral->GetMaximum(),max);
+    }
+
+
     // Get graph of 'other' data - e.g., UA5, CMS, ... - and check if
     // there's any graphs.  Set the pad division based on that.
     // Info("DrawIt", "Will %sdraw other results", (doComp ? "" : "not "));
index 4dadd409a539920d97d91d724f83c41ea6997ae6..13cf83df410d91b01e86611ecb80ef947a6da49e 100644 (file)
@@ -60,7 +60,7 @@ void MakeAOD(const char* esddir,
                                  "EMCALTrigger "
                                  "SPDPileupVertices " 
                                  "TrkPileupVertices " 
-                                 "Tracks "
+                                 // "Tracks "
                                  "MuonTracks " 
                                  "PmdTracks "
                                  "TrdTracks "
@@ -101,6 +101,11 @@ void MakeAOD(const char* esddir,
   Compile("$ALICE_ROOT/PWG2/FORWARD/analysis2/scripts/AliESDCentrality.C","+g");
   // gDebug = 0;
   AddTaskCopyHeader();
+
+
+  // Central multiplicity
+  Compile("$ALICE_ROOT/PWG2/FORWARD/analysis2/AddTaskCentralMult.C","+");
+  AddTaskCentralMult();
 #endif
 
   // FMD 
diff --git a/PWG2/FORWARD/analysis2/scripts/AliESDCentrality.C b/PWG2/FORWARD/analysis2/scripts/AliESDCentrality.C
new file mode 100644 (file)
index 0000000..9e4299c
--- /dev/null
@@ -0,0 +1,169 @@
+//-*- Mode: C++ -*-
+#ifndef ALIESDCentrality_H
+#define ALIESDCentrality_H
+/* This file is property of and copyright by the ALICE HLT Project        *
+ * ALICE Experiment at CERN, All rights reserved.                         *
+ * See cxx source for full Copyright notice                               */
+
+//*****************************************************
+//   Class AliCentralitySelectionTask
+//   author: Alberica Toia
+//*****************************************************
+
+#include "TNamed.h"
+
+class AliESDCentrality : public TNamed
+{
+ public:
+
+  AliESDCentrality();  /// constructor
+  ~AliESDCentrality();  /// destructor
+  AliESDCentrality(const AliESDCentrality& cnt); /// copy constructor
+  AliESDCentrality& operator=(const AliESDCentrality& cnt);   /// assignment operator
+
+  /// set centrality result
+  void SetCentralityV0M(Float_t cent) {fCentralityV0M = cent;} 
+  void SetCentralityFMD(Float_t cent) {fCentralityFMD = cent;}
+  void SetCentralityTRK(Float_t cent) {fCentralityTRK = cent;}
+  void SetCentralityTKL(Float_t cent) {fCentralityTKL = cent;}
+  void SetCentralityCL0(Float_t cent) {fCentralityCL0 = cent;}
+  void SetCentralityV0MvsFMD(Float_t cent) {fCentralityV0MvsFMD = cent;}
+  void SetCentralityTKLvsV0M(Float_t cent) {fCentralityTKLvsV0M = cent;}
+  void SetCentralityZEMvsZDC(Float_t cent) {fCentralityZEMvsZDC = cent;}
+
+  /// get centrality result
+  Float_t GetCentralityPercentile(const char *method);
+  Int_t   GetCentralityClass10(const char *method);
+  Int_t   GetCentralityClass5(const char *method);
+  Bool_t  IsEventInCentralityClass(Float_t a, Float_t b, const char *method);
+
+ private:
+  Float_t fCentralityV0M;   // Centrality from V0
+  Float_t fCentralityFMD;   // Centrality from FMD
+  Float_t fCentralityTRK;   // Centrality from tracks
+  Float_t fCentralityTKL;   // Centrality from tracklets
+  Float_t fCentralityCL0;   // Centrality from Clusters in layer 0
+  Float_t fCentralityV0MvsFMD;   // Centrality from V0 vs FMD
+  Float_t fCentralityTKLvsV0M;   // Centrality from tracklets vs V0
+  Float_t fCentralityZEMvsZDC;   // Centrality from ZEM vs ZDC
+
+  ClassDef(AliESDCentrality, 2)
+};
+
+inline
+AliESDCentrality::AliESDCentrality() : TNamed("ESDCentrality", "Centrality"),
+  fCentralityV0M(0),
+  fCentralityFMD(0),
+  fCentralityTRK(0),
+  fCentralityTKL(0),
+  fCentralityCL0(0),
+  fCentralityV0MvsFMD(0),
+  fCentralityTKLvsV0M(0),
+  fCentralityZEMvsZDC(0)
+{
+  /// constructor
+}
+
+inline
+AliESDCentrality::AliESDCentrality(const AliESDCentrality& cnt) : 
+  TNamed(cnt), 
+  fCentralityV0M(cnt.fCentralityV0M),
+  fCentralityFMD(cnt.fCentralityFMD),
+  fCentralityTRK(cnt.fCentralityTRK),
+  fCentralityTKL(cnt.fCentralityTKL),
+  fCentralityCL0(cnt.fCentralityCL0),
+  fCentralityV0MvsFMD(cnt.fCentralityV0MvsFMD),
+  fCentralityTKLvsV0M(cnt.fCentralityTKLvsV0M),
+  fCentralityZEMvsZDC(cnt.fCentralityZEMvsZDC)
+{
+  /// Copy constructor
+}
+
+inline
+AliESDCentrality& AliESDCentrality::operator=(const AliESDCentrality& c)
+{
+  /// Assignment operator
+  if (this!=&c) {
+    TNamed::operator=(c);
+    fCentralityV0M = c.fCentralityV0M;
+    fCentralityFMD = c.fCentralityFMD;
+    fCentralityTRK = c.fCentralityTRK;
+    fCentralityTKL = c.fCentralityTKL;
+    fCentralityCL0 = c.fCentralityCL0;
+    fCentralityV0MvsFMD = c.fCentralityV0MvsFMD;
+    fCentralityTKLvsV0M = c.fCentralityTKLvsV0M;
+    fCentralityZEMvsZDC = c.fCentralityZEMvsZDC;
+  }
+
+  return *this;
+}
+
+inline
+AliESDCentrality::~AliESDCentrality()
+{
+  /// destructor
+}
+
+inline
+Float_t AliESDCentrality::GetCentralityPercentile(const char *x)
+{
+  TString method = x;
+  if(method.CompareTo("V0M")==0)      return fCentralityV0M;
+  if(method.CompareTo("FMD")==0)      return fCentralityFMD;
+  if(method.CompareTo("TRK")==0)      return fCentralityTRK;
+  if(method.CompareTo("TKL")==0)      return fCentralityTKL;
+  if(method.CompareTo("CL0")==0)      return fCentralityCL0;
+  if(method.CompareTo("V0MvsFMD")==0) return fCentralityV0MvsFMD;
+  if(method.CompareTo("TKLvsV0M")==0) return fCentralityTKLvsV0M;
+  if(method.CompareTo("ZENvsZDC")==0) return fCentralityZEMvsZDC;
+  return -1;
+}
+
+inline
+Int_t AliESDCentrality::GetCentralityClass10(const char *x)
+{
+  TString method = x;
+  if(method.CompareTo("V0M")==0)      return (Int_t) (fCentralityV0M / 10.0);
+  if(method.CompareTo("FMD")==0)      return (Int_t) (fCentralityFMD / 10.0);
+  if(method.CompareTo("TRK")==0)      return (Int_t) (fCentralityTRK / 10.0);
+  if(method.CompareTo("TKL")==0)      return (Int_t) (fCentralityTKL / 10.0);
+  if(method.CompareTo("CL0")==0)      return (Int_t) (fCentralityCL0 / 10.0);
+  if(method.CompareTo("V0MvsFMD")==0) return (Int_t) (fCentralityV0MvsFMD / 10.0);
+  if(method.CompareTo("TKLvsV0M")==0) return (Int_t) (fCentralityTKLvsV0M / 10.0);
+  if(method.CompareTo("ZENvsZDC")==0) return (Int_t) (fCentralityZEMvsZDC / 10.0);
+  return -1;
+}
+
+Int_t AliESDCentrality::GetCentralityClass5(const char *x)
+{
+ TString method = x;
+  if(method.CompareTo("V0M")==0)      return (Int_t) (fCentralityV0M / 5.0);
+  if(method.CompareTo("FMD")==0)      return (Int_t) (fCentralityFMD / 5.0);
+  if(method.CompareTo("TRK")==0)      return (Int_t) (fCentralityTRK / 5.0);
+  if(method.CompareTo("TKL")==0)      return (Int_t) (fCentralityTKL / 5.0);
+  if(method.CompareTo("CL0")==0)      return (Int_t) (fCentralityCL0 / 5.0);
+  if(method.CompareTo("V0MvsFMD")==0) return (Int_t) (fCentralityV0MvsFMD / 5.0);
+  if(method.CompareTo("TKLvsV0M")==0) return (Int_t) (fCentralityTKLvsV0M / 5.0);
+  if(method.CompareTo("ZENvsZDC")==0) return (Int_t) (fCentralityZEMvsZDC / 5.0);
+  return -1;
+}
+
+inline
+Bool_t AliESDCentrality::IsEventInCentralityClass(Float_t a, Float_t b, const char *x)
+{
+  TString method = x;
+  if ((method.CompareTo("V0M")==0) && (fCentralityV0M >=a && fCentralityV0M < b)) return kTRUE;
+  if ((method.CompareTo("FMD")==0) && (fCentralityFMD >=a && fCentralityFMD < b)) return kTRUE;
+  if ((method.CompareTo("TRK")==0) && (fCentralityTRK >=a && fCentralityTRK < b)) return kTRUE;
+  if ((method.CompareTo("TKL")==0) && (fCentralityTKL >=a && fCentralityTKL < b)) return kTRUE;
+  if ((method.CompareTo("CL0")==0) && (fCentralityCL0 >=a && fCentralityCL0 < b)) return kTRUE;
+  if ((method.CompareTo("V0MvsFMD")==0) && (fCentralityV0MvsFMD >=a && fCentralityV0MvsFMD < b)) return kTRUE;
+  if ((method.CompareTo("TKLvsV0M")==0) && (fCentralityTKLvsV0M >=a && fCentralityTKLvsV0M < b)) return kTRUE;
+  if ((method.CompareTo("ZEMvsZDC")==0) && (fCentralityZEMvsZDC >=a && fCentralityZEMvsZDC < b)) return kTRUE;
+  else return kFALSE;
+}
+
+
+
+
+#endif //ALIESDCENTRALITY_H