PWG0/dNdEta/AlidNdEtaAnalysisMCSelector.cxx

   1 /* $Id$ */
   2
   3 #include "AlidNdEtaAnalysisMCSelector.h"
   4
   5 #include <TStyle.h>
   6 #include <TSystem.h>
   7 #include <TCanvas.h>
   8 #include <TParticle.h>
   9 #include <TParticlePDG.h>
  10 #include <TVector3.h>
  11 #include <TH1F.h>
  12 #include <TH3F.h>
  13 #include <TTree.h>
  14 #include <TFile.h>
  15
  16 #include <AliLog.h>
  17 #include <AliGenEventHeader.h>
  18 #include <AliHeader.h>
  19 #include <AliStack.h>
  20 #include <AliESDtrack.h>
  21
  22 #include "dNdEta/dNdEtaAnalysis.h"
  23 #include <AliPWG0Helper.h>
  24 #include <AliCorrection.h>
  25 #include <AliCorrectionMatrix2D.h>
  26 #include "esdTrackCuts/AliESDtrackCuts.h"
  27
  28
  29 ClassImp(AlidNdEtaAnalysisMCSelector)
  30
  31 AlidNdEtaAnalysisMCSelector::AlidNdEtaAnalysisMCSelector() :
  32   AliSelectorRL(),
  33   fEsdTrackCuts(0),
  34   fdNdEtaAnalysis(0),
  35   fdNdEtaAnalysisTr(0),
  36   fdNdEtaAnalysisTrVtx(0),
  37   fdNdEtaAnalysisTracks(0),
  38   fVertex(0),
  39   fPartPt(0),
  40   fEvents(0)
  41 {
  42   //
  43   // Constructor. Initialization of pointers
  44   //
  45 }
  46
  47 AlidNdEtaAnalysisMCSelector::~AlidNdEtaAnalysisMCSelector()
  48 {
  49   //
  50   // Destructor
  51   //
  52 }
  53
  54 void AlidNdEtaAnalysisMCSelector::Begin(TTree* tree)
  55 {
  56   // Begin function
  57
  58   ReadUserObjects(tree);
  59 }
  60
  61 void AlidNdEtaAnalysisMCSelector::ReadUserObjects(TTree* tree)
  62 {
  63   // read the user objects, called from slavebegin and begin
  64
  65   if (!fEsdTrackCuts && fInput)
  66     fEsdTrackCuts = dynamic_cast<AliESDtrackCuts*> (fInput->FindObject("AliESDtrackCuts"));
  67
  68   if (!fEsdTrackCuts && tree)
  69     fEsdTrackCuts = dynamic_cast<AliESDtrackCuts*> (tree->GetUserInfo()->FindObject("AliESDtrackCuts"));
  70
  71   if (!fEsdTrackCuts)
  72      AliDebug(AliLog::kError, "ERROR: Could not read EsdTrackCuts from input list.");
  73 }
  74
  75 void AlidNdEtaAnalysisMCSelector::SlaveBegin(TTree * tree)
  76 {
  77   // The SlaveBegin() function is called after the Begin() function.
  78   // When running with PROOF SlaveBegin() is called on each slave server.
  79   // The tree argument is deprecated (on PROOF 0 is passed).
  80
  81   AliSelectorRL::SlaveBegin(tree);
  82
  83   ReadUserObjects(tree);
  84
  85   fdNdEtaAnalysis = new dNdEtaAnalysis("dndeta", "dndeta");
  86   fdNdEtaAnalysisTr = new dNdEtaAnalysis("dndetaTr", "dndetaTr");
  87   fdNdEtaAnalysisTrVtx = new dNdEtaAnalysis("dndetaTrVtx", "dndetaTrVtx");
  88   fdNdEtaAnalysisTracks = new dNdEtaAnalysis("dndetaTracks", "dndetaTracks");
  89   fVertex = new TH3F("vertex_check", "vertex_check", 50, -50, 50, 50, -50, 50, 50, -50, 50);
  90   for (Int_t i=0; i<3; ++i)
  91   {
  92     fPartEta[i] = new TH1F(Form("dndeta_check_%d", i), Form("dndeta_check_%d", i), 60, -6, 6);
  93     fPartEta[i]->Sumw2();
  94   }
  95   fPartPt =  new TH1F("dndeta_check_pt", "dndeta_check_pt", 1000, 0, 10);
  96   fPartPt->Sumw2();
  97   fEvents = new TH1F("dndeta_check_vertex", "dndeta_check_vertex", 40, -20, 20);
  98 }
  99
 100 void AlidNdEtaAnalysisMCSelector::Init(TTree *tree)
 101 {
 102   AliSelectorRL::Init(tree);
 103
 104   if (!tree)
 105   {
 106     AliDebug(AliLog::kError, "tree not available");
 107     return;
 108   }
 109
 110   tree->SetBranchStatus("*", 0);
 111   tree->SetBranchStatus("fTriggerMask", 1);
 112   tree->SetBranchStatus("fSPDVertex*", 1);
 113   tree->SetBranchStatus("fTracks.fLabel", 1);
 114
 115   AliESDtrackCuts::EnableNeededBranches(tree);
 116 }
 117
 118 Bool_t AlidNdEtaAnalysisMCSelector::Process(Long64_t entry)
 119 {
 120   // fill the dNdEtaAnalysis class from the monte carlo
 121
 122   if (AliSelectorRL::Process(entry) == kFALSE)
 123     return kFALSE;
 124
 125   // Check prerequisites
 126   if (!fESD)
 127   {
 128     AliDebug(AliLog::kError, "ESD branch not available");
 129     return kFALSE;
 130   }
 131
 132   AliStack* stack = GetStack();
 133   if (!stack)
 134   {
 135     AliDebug(AliLog::kError, "Stack not available");
 136     return kFALSE;
 137   }
 138
 139   AliHeader* header = GetHeader();
 140   if (!header)
 141   {
 142     AliDebug(AliLog::kError, "Header not available");
 143     return kFALSE;
 144   }
 145
 146   Bool_t vertexReconstructed = AliPWG0Helper::IsVertexReconstructed(fESD);
 147   Bool_t eventTriggered = AliPWG0Helper::IsEventTriggered(fESD);
 148
 149   // get the MC vertex
 150   AliGenEventHeader* genHeader = header->GenEventHeader();
 151
 152   TArrayF vtxMC(3);
 153   genHeader->PrimaryVertex(vtxMC);
 154
 155   // loop over mc particles
 156   Int_t nPrim  = stack->GetNprimary();
 157
 158   Int_t nAcceptedParticles = 0;
 159
 160   for (Int_t iMc = 0; iMc < nPrim; ++iMc)
 161   {
 162     TParticle* particle = stack->Particle(iMc);
 163
 164     if (!particle)
 165       continue;
 166
 167     if (AliPWG0Helper::IsPrimaryCharged(particle, nPrim) == kFALSE)
 168       continue;
 169
 170     AliDebug(AliLog::kDebug+1, Form("Accepted primary %d, unique ID: %d", iMc, particle->GetUniqueID()));
 171     ++nAcceptedParticles;
 172
 173     fdNdEtaAnalysis->FillTrack(vtxMC[2], particle->Eta(), particle->Pt());
 174     fVertex->Fill(particle->Vx(), particle->Vy(), particle->Vz());
 175
 176     if (eventTriggered)
 177     {
 178       fdNdEtaAnalysisTr->FillTrack(vtxMC[2], particle->Eta(), particle->Pt());
 179       if (vertexReconstructed)
 180         fdNdEtaAnalysisTrVtx->FillTrack(vtxMC[2], particle->Eta(), particle->Pt());
 181     }
 182
 183     if (TMath::Abs(vtxMC[2]) < 20)
 184     {
 185       fPartEta[0]->Fill(particle->Eta());
 186
 187       if (vtxMC[2] < 0)
 188         fPartEta[1]->Fill(particle->Eta());
 189       else
 190         fPartEta[2]->Fill(particle->Eta());
 191     }
 192
 193     if (TMath::Abs(particle->Eta()) < 0.8)
 194       fPartPt->Fill(particle->Pt());
 195   }
 196
 197   fEvents->Fill(vtxMC[2]);
 198
 199   fdNdEtaAnalysis->FillEvent(vtxMC[2], nAcceptedParticles);
 200   if (eventTriggered)
 201   {
 202     fdNdEtaAnalysisTr->FillEvent(vtxMC[2], nAcceptedParticles);
 203     if (vertexReconstructed)
 204       fdNdEtaAnalysisTrVtx->FillEvent(vtxMC[2], nAcceptedParticles);
 205   }
 206
 207   if (!eventTriggered || !vertexReconstructed)
 208     return kTRUE;
 209
 210   // from tracks is only done for triggered and vertex reconstructed events
 211
 212   // get number of "good" tracks
 213   TObjArray* list = fEsdTrackCuts->GetAcceptedTracks(fESD);
 214   Int_t nGoodTracks = list->GetEntries();
 215
 216   // loop over esd tracks
 217   for (Int_t t=0; t<nGoodTracks; t++)
 218   {
 219     AliESDtrack* esdTrack = dynamic_cast<AliESDtrack*> (list->At(t));
 220     if (!esdTrack)
 221     {
 222       AliDebug(AliLog::kError, Form("ERROR: Could not retrieve track %d.", t));
 223       continue;
 224     }
 225
 226     Int_t label = TMath::Abs(esdTrack->GetLabel());
 227
 228     TParticle* particle = stack->Particle(label);
 229     if (!particle)
 230     {
 231       AliDebug(AliLog::kError, Form("ERROR: Could not retrieve particle %d.", esdTrack->GetLabel()));
 232       continue;
 233     }
 234
 235     fdNdEtaAnalysisTracks->FillTrack(vtxMC[2], particle->Eta(), particle->Pt());
 236   } // end of track loop
 237
 238   delete list;
 239   list = 0;
 240
 241   // for event count per vertex
 242   fdNdEtaAnalysisTracks->FillEvent(vtxMC[2], nGoodTracks);
 243
 244   return kTRUE;
 245 }
 246
 247 void AlidNdEtaAnalysisMCSelector::SlaveTerminate()
 248 {
 249   // The SlaveTerminate() function is called after all entries or objects
 250   // have been processed. When running with PROOF SlaveTerminate() is called
 251   // on each slave server.
 252
 253   AliSelectorRL::SlaveTerminate();
 254
 255   // Add the histograms to the output on each slave server
 256   if (!fOutput)
 257   {
 258     AliDebug(AliLog::kError, Form("ERROR: Output list not initialized."));
 259     return;
 260   }
 261
 262   fOutput->Add(fdNdEtaAnalysis);
 263   fOutput->Add(fdNdEtaAnalysisTr);
 264   fOutput->Add(fdNdEtaAnalysisTrVtx);
 265   fOutput->Add(fdNdEtaAnalysisTracks);
 266
 267   fOutput->Add(fPartPt);
 268   fOutput->Add(fEvents);
 269   for (Int_t i=0; i<3; ++i)
 270     fOutput->Add(fPartEta[i]);
 271 }
 272
 273 void AlidNdEtaAnalysisMCSelector::Terminate()
 274 {
 275   //
 276
 277   AliSelectorRL::Terminate();
 278
 279   fdNdEtaAnalysis = dynamic_cast<dNdEtaAnalysis*> (fOutput->FindObject("dndeta"));
 280   fdNdEtaAnalysisTr = dynamic_cast<dNdEtaAnalysis*> (fOutput->FindObject("dndetaTr"));
 281   fdNdEtaAnalysisTrVtx = dynamic_cast<dNdEtaAnalysis*> (fOutput->FindObject("dndetaTrVtx"));
 282   fdNdEtaAnalysisTracks = dynamic_cast<dNdEtaAnalysis*> (fOutput->FindObject("dndetaTracks"));
 283   fPartPt = dynamic_cast<TH1F*> (fOutput->FindObject("dndeta_check_pt"));
 284   fEvents = dynamic_cast<TH1F*> (fOutput->FindObject("dndeta_check_vertex"));
 285   for (Int_t i=0; i<3; ++i)
 286     fPartEta[i] = dynamic_cast<TH1F*> (fOutput->FindObject(Form("dndeta_check_%d", i)));
 287
 288   if (!fdNdEtaAnalysis || !fdNdEtaAnalysisTr || !fdNdEtaAnalysisTrVtx || !fPartPt || !fEvents)
 289   {
 290     AliDebug(AliLog::kError, Form("ERROR: Histograms not available %p %p", (void*) fdNdEtaAnalysis, (void*) fPartPt));
 291     return;
 292   }
 293
 294   fdNdEtaAnalysis->Finish(0, -1, AlidNdEtaCorrection::kNone);
 295   fdNdEtaAnalysisTr->Finish(0, -1, AlidNdEtaCorrection::kNone);
 296   fdNdEtaAnalysisTrVtx->Finish(0, -1, AlidNdEtaCorrection::kNone);
 297   fdNdEtaAnalysisTracks->Finish(0, -1, AlidNdEtaCorrection::kNone);
 298
 299   Int_t events = (Int_t) fdNdEtaAnalysis->GetData()->GetEventCorrection()->GetMeasuredHistogram()->Integral();
 300   fPartPt->Scale(1.0/events);
 301   fPartPt->Scale(1.0/fPartPt->GetBinWidth(1));
 302
 303   if (fPartEta[0])
 304   {
 305     Int_t events1 = (Int_t) fEvents->Integral(fEvents->GetXaxis()->FindBin(-19.9), fEvents->GetXaxis()->FindBin(-0.001));
 306     Int_t events2 = (Int_t) fEvents->Integral(fEvents->GetXaxis()->FindBin(0.001), fEvents->GetXaxis()->FindBin(19.9));
 307
 308     fPartEta[0]->Scale(1.0 / (events1 + events2));
 309     fPartEta[1]->Scale(1.0 / events1);
 310     fPartEta[2]->Scale(1.0 / events2);
 311
 312     for (Int_t i=0; i<3; ++i)
 313       fPartEta[i]->Scale(1.0 / fPartEta[i]->GetBinWidth(1));
 314
 315     new TCanvas("control", "control", 500, 500);
 316     for (Int_t i=0; i<3; ++i)
 317     {
 318       fPartEta[i]->SetLineColor(i+1);
 319       fPartEta[i]->Draw((i==0) ? "" : "SAME");
 320     }
 321   }
 322
 323   TFile* fout = new TFile("analysis_mc.root","RECREATE");
 324
 325   fdNdEtaAnalysis->SaveHistograms();
 326   fdNdEtaAnalysisTr->SaveHistograms();
 327   fdNdEtaAnalysisTrVtx->SaveHistograms();
 328   fdNdEtaAnalysisTracks->SaveHistograms();
 329   fPartPt->Write();
 330
 331   fout->Write();
 332   fout->Close();
 333
 334   if (fPartPt)
 335   {
 336     new TCanvas("control2", "control2", 500, 500);
 337     fPartPt->DrawCopy();
 338   }
 339
 340   if (fEvents)
 341   {
 342     new TCanvas("control3", "control3", 500, 500);
 343     fEvents->DrawCopy();
 344   }
 345 }