PWG0/dNdEta/AliMultiplicityCorrectionSelector.cxx

   1 /* $Id$ */
   2
   3 #include "AliMultiplicityCorrectionSelector.h"
   4
   5 #include <TStyle.h>
   6 #include <TSystem.h>
   7 #include <TCanvas.h>
   8 #include <TParticle.h>
   9 #include <TParticlePDG.h>
  10 #include <TH1F.h>
  11
  12 #include <TChain.h>
  13 #include <TSelector.h>
  14 #include <TFile.h>
  15
  16 #include <AliLog.h>
  17 #include <AliTracker.h>
  18 #include <AliESDVertex.h>
  19 #include <AliESD.h>
  20 #include <AliESDtrack.h>
  21 #include <AliRunLoader.h>
  22 #include <AliStack.h>
  23
  24 #include <AliMultiplicity.h>
  25
  26 #include <AliHeader.h>
  27 #include <AliGenEventHeader.h>
  28 #include <../PYTHIA6/AliGenPythiaEventHeader.h>
  29 #include <../EVGEN/AliGenCocktailEventHeader.h>
  30
  31 #include "esdTrackCuts/AliESDtrackCuts.h"
  32 #include "dNdEta/AliMultiplicityCorrection.h"
  33 #include "AliPWG0Helper.h"
  34 #include "AliPWG0depHelper.h"
  35
  36 ClassImp(AliMultiplicityCorrectionSelector)
  37
  38 AliMultiplicityCorrectionSelector::AliMultiplicityCorrectionSelector() :
  39   AliSelectorRL()
  40 {
  41   //
  42   // Constructor. Initialization of pointers
  43   //
  44 }
  45
  46 AliMultiplicityCorrectionSelector::~AliMultiplicityCorrectionSelector()
  47 {
  48   //
  49   // Destructor
  50   //
  51
  52   // histograms are in the output list and deleted when the output
  53   // list is deleted by the TSelector dtor
  54 }
  55
  56
  57 void AliMultiplicityCorrectionSelector::ReadUserObjects(TTree* tree)
  58 {
  59   // read the user objects, called from slavebegin and begin
  60
  61   if (!fEsdTrackCuts && fInput)
  62     fEsdTrackCuts = dynamic_cast<AliESDtrackCuts*> (fInput->FindObject("AliESDtrackCuts"));
  63
  64   if (!fEsdTrackCuts && tree)
  65     fEsdTrackCuts = dynamic_cast<AliESDtrackCuts*> (tree->GetUserInfo()->FindObject("AliESDtrackCuts"));
  66
  67   if (!fEsdTrackCuts)
  68      AliDebug(AliLog::kError, "ERROR: Could not read EsdTrackCuts from input list.");
  69 }
  70
  71 void AliMultiplicityCorrectionSelector::Begin(TTree * tree)
  72 {
  73   // The Begin() function is called at the start of the query.
  74   // When running with PROOF Begin() is only called on the client.
  75   // The tree argument is deprecated (on PROOF 0 is passed).
  76
  77   AliSelectorRL::Begin(tree);
  78
  79   ReadUserObjects(tree);
  80
  81   TString option = GetOption();
  82   AliInfo(Form("Called with option %s.", option.Data()));
  83
  84 }
  85
  86 void AliMultiplicityCorrectionSelector::SlaveBegin(TTree * tree)
  87 {
  88   // The SlaveBegin() function is called after the Begin() function.
  89   // When running with PROOF SlaveBegin() is called on each slave server.
  90   // The tree argument is deprecated (on PROOF 0 is passed).
  91
  92   AliSelectorRL::SlaveBegin(tree);
  93
  94   ReadUserObjects(tree);
  95
  96   fMultiplicityCorrection = new AliMultiplicityCorrection("mult_correction", "mult_correction");
  97
  98 }
  99
 100 void AliMultiplicityCorrectionSelector::Init(TTree* tree)
 101 {
 102   // read the user objects
 103
 104   AliSelectorRL::Init(tree);
 105
 106   // Enable only the needed branches
 107   if (tree)
 108   {
 109     tree->SetBranchStatus("*", 0);
 110     tree->SetBranchStatus("fTriggerMask", 1);
 111     tree->SetBranchStatus("fSPDVertex*", 1);
 112     tree->SetBranchStatus("fTracks.fLabel", 1);
 113     tree->SetBranchStatus("fTracks.fITSncls", 1);
 114     tree->SetBranchStatus("fTracks.fTPCncls", 1);
 115     tree->SetBranchStatus("fSPDMult*", 1);
 116
 117     AliESDtrackCuts::EnableNeededBranches(tree);
 118   }
 119 }
 120
 121 Bool_t AliMultiplicityCorrectionSelector::Process(Long64_t entry)
 122 {
 123   // The Process() function is called for each entry in the tree (or possibly
 124   // keyed object in the case of PROOF) to be processed. The entry argument
 125   // specifies which entry in the currently loaded tree is to be processed.
 126   // It can be passed to either TTree::GetEntry() or TBranch::GetEntry()
 127   // to read either all or the required parts of the data. When processing
 128   // keyed objects with PROOF, the object is already loaded and is available
 129   // via the fObject pointer.
 130   //
 131   // This function should contain the "body" of the analysis. It can contain
 132   // simple or elaborate selection criteria, run algorithms on the data
 133   // of the event and typically fill histograms.
 134
 135   // WARNING when a selector is used with a TChain, you must use
 136   //  the pointer to the current TTree to call GetEntry(entry).
 137   //  The entry is always the local entry number in the current tree.
 138   //  Assuming that fTree is the pointer to the TChain being processed,
 139   //  use fTree->GetTree()->GetEntry(entry).
 140
 141   AliDebug(AliLog::kDebug+1,"Processing event ...\n");
 142
 143   if (AliSelectorRL::Process(entry) == kFALSE)
 144     return kFALSE;
 145
 146   // check prerequesites
 147   if (!fESD)
 148   {
 149     AliDebug(AliLog::kError, "ESD branch not available");
 150     return kFALSE;
 151   }
 152
 153   AliHeader* header = GetHeader();
 154   if (!header)
 155   {
 156     AliDebug(AliLog::kError, "Header not available");
 157     return kFALSE;
 158   }
 159
 160   // getting the stack
 161   AliStack* stack = GetStack();
 162   if (!stack)
 163   {
 164     AliDebug(AliLog::kError, "Stack not available");
 165     return kFALSE;
 166   }
 167
 168   // getting the its multiplicity data
 169   AliMultiplicity* itsMult = (AliMultiplicity*)fESD->GetMultiplicity();
 170   if (!itsMult) {
 171     AliDebug(AliLog::kError, "Multiplicity object not found in ESD");
 172     return kFALSE;
 173   }
 174
 175   //  if (!fEsdTrackCuts)
 176   // {
 177   //  AliDebug(AliLog::kError, "fESDTrackCuts not available");
 178   //  return kFALSE;
 179   // }
 180
 181   Bool_t vertexReconstructed = AliPWG0Helper::IsVertexReconstructed(fESD);
 182
 183   Bool_t eventTriggered = AliPWG0Helper::IsEventTriggered(fESD);
 184
 185   // only look at triggered events with a vertex
 186   if (!(vertexReconstructed && eventTriggered))
 187     return kTRUE;
 188
 189   // get the MC vertex
 190   AliGenEventHeader* genHeader = header->GenEventHeader();
 191
 192   TArrayF vtxMC(3);
 193   genHeader->PrimaryVertex(vtxMC);
 194
 195   // vertex cut - should not be hard coded!!!
 196   if (TMath::Abs(vtxMC[2]>10))
 197     return kTRUE;
 198
 199   // loop over mc particles
 200   Int_t nPrim  = stack->GetNprimary();
 201
 202   for (Int_t iMc = 0; iMc < nPrim; ++iMc) {
 203
 204     AliDebug(AliLog::kDebug+1, Form("MC Loop: Processing particle %d.", iMc));
 205
 206     TParticle* particle = stack->Particle(iMc);
 207
 208     if (!particle)
 209     {
 210       AliDebug(AliLog::kError, Form("UNEXPECTED: particle with label %d not found in stack (mc loop).", iMc));
 211       continue;
 212     }
 213
 214     if (AliPWG0Helper::IsPrimaryCharged(particle, nPrim) == kFALSE)
 215       continue;
 216
 217     if (particle->GetPDG()->Charge()==0)
 218       continue;
 219
 220     Float_t eta = particle->Eta();
 221
 222     fMultiplicityCorrection->FillMeasuredMultHit(eta);
 223
 224   }// end of mc particle
 225
 226   // ########################################################
 227   // loop over spd tracklets
 228
 229
 230   for(Int_t i=0; i<itsMult->GetNumberOfTracklets(); i++) {
 231     Float_t theta  = itsMult->GetTheta(i);
 232     Float_t eta    = -TMath::Log(TMath::Tan(0.5*theta));
 233
 234     fMultiplicityCorrection->FillTrueMultHit(eta);
 235   }
 236
 237   // very important!
 238   fMultiplicityCorrection->NewEvent();
 239
 240   return kTRUE;
 241 }
 242
 243 void AliMultiplicityCorrectionSelector::SlaveTerminate()
 244 {
 245   // The SlaveTerminate() function is called after all entries or objects
 246   // have been processed. When running with PROOF SlaveTerminate() is called
 247   // on each slave server.
 248
 249   AliSelectorRL::SlaveTerminate();
 250
 251   // Add the histograms to the output on each slave server
 252   if (!fOutput)
 253   {
 254     AliDebug(AliLog::kError, "ERROR: Output list not initialized");
 255     return;
 256   }
 257
 258   fOutput->Add(fMultiplicityCorrection);
 259
 260   AliDebug(AliLog::kDebug+1,"Slave terminate ...\n");
 261
 262 }
 263
 264 void AliMultiplicityCorrectionSelector::Terminate()
 265 {
 266   // The Terminate() function is the last function to be called during
 267   // a query. It always runs on the client, it can be used to present
 268   // the results graphically or save the results to file.
 269
 270   AliDebug(AliLog::kDebug+1,"Terminate ...\n");
 271
 272
 273   AliSelectorRL::Terminate();
 274
 275   fMultiplicityCorrection = dynamic_cast<AliMultiplicityCorrection*> (fOutput->FindObject("mult_correction"));
 276   if (!fMultiplicityCorrection)
 277   {
 278     AliDebug(AliLog::kError, "Could not read object from output list");
 279     return;
 280   }
 281
 282   fMultiplicityCorrection->Finish();
 283
 284   TFile* fout = new TFile(Form("correction_mult%s.root", GetOption()), "RECREATE");
 285
 286   //  if (fEsdTrackCuts)
 287   //  fEsdTrackCuts->SaveHistograms("esd_track_cuts");
 288   fMultiplicityCorrection->SaveHistograms();
 289
 290   fout->Write();
 291   fout->Close();
 292
 293   fMultiplicityCorrection->DrawHistograms();
 294
 295 }