PWG0/dNdEta/AlidNdEtaCorrectionSelector.cxx

   1 /* $Id$ */
   2
   3 #include "AlidNdEtaCorrectionSelector.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 <AliHeader.h>
  25 #include <AliGenEventHeader.h>
  26 #include <../STEER/AliGenPythiaEventHeader.h>
  27 #include <../STEER/AliGenCocktailEventHeader.h>
  28
  29 #include "esdTrackCuts/AliESDtrackCuts.h"
  30 #include "dNdEta/AlidNdEtaCorrection.h"
  31 #include "AliPWG0Helper.h"
  32
  33 #include "dNdEta/dNdEtaAnalysis.h"
  34
  35 ClassImp(AlidNdEtaCorrectionSelector)
  36
  37 AlidNdEtaCorrectionSelector::AlidNdEtaCorrectionSelector() :
  38   AliSelectorRL(),
  39   fEsdTrackCuts(0),
  40   fdNdEtaCorrection(0),
  41   fdNdEtaAnalysisMC(0),
  42   fdNdEtaAnalysisESD(0),
  43   fPIDParticles(0),
  44   fPIDTracks(0),
  45   fClustersITSPos(0),
  46   fClustersTPCPos(0),
  47   fClustersITSNeg(0),
  48   fClustersTPCNeg(0),
  49   fSignMode(0)
  50 {
  51   //
  52   // Constructor. Initialization of pointers
  53   //
  54 }
  55
  56 AlidNdEtaCorrectionSelector::~AlidNdEtaCorrectionSelector()
  57 {
  58   //
  59   // Destructor
  60   //
  61
  62   // histograms are in the output list and deleted when the output
  63   // list is deleted by the TSelector dtor
  64 }
  65
  66 Bool_t AlidNdEtaCorrectionSelector::SignOK(TParticlePDG* particle)
  67 {
  68   // returns if a particle with this sign should be counted
  69   // this is determined by the value of fSignMode, which should have the same sign
  70   // as the charge
  71   // if fSignMode is 0 all particles are counted
  72
  73   if (fSignMode == 0)
  74     return kTRUE;
  75
  76   if (!particle)
  77   {
  78     printf("WARNING: not counting a particle that does not have a pdg particle\n");
  79     return kFALSE;
  80   }
  81
  82   Double_t charge = particle->Charge();
  83
  84   if (fSignMode > 0)
  85     if (charge < 0)
  86       return kFALSE;
  87
  88   if (fSignMode < 0)
  89     if (charge > 0)
  90       return kFALSE;
  91
  92   return kTRUE;
  93 }
  94
  95 void AlidNdEtaCorrectionSelector::ReadUserObjects(TTree* tree)
  96 {
  97   // read the user objects, called from slavebegin and begin
  98
  99   if (!fEsdTrackCuts && fInput)
 100     fEsdTrackCuts = dynamic_cast<AliESDtrackCuts*> (fInput->FindObject("AliESDtrackCuts"));
 101
 102   if (!fEsdTrackCuts && tree)
 103     fEsdTrackCuts = dynamic_cast<AliESDtrackCuts*> (tree->GetUserInfo()->FindObject("AliESDtrackCuts"));
 104
 105   if (!fEsdTrackCuts)
 106      AliDebug(AliLog::kError, "ERROR: Could not read EsdTrackCuts from input list.");
 107 }
 108
 109 void AlidNdEtaCorrectionSelector::Begin(TTree * tree)
 110 {
 111   // The Begin() function is called at the start of the query.
 112   // When running with PROOF Begin() is only called on the client.
 113   // The tree argument is deprecated (on PROOF 0 is passed).
 114
 115   AliSelectorRL::Begin(tree);
 116
 117   ReadUserObjects(tree);
 118
 119   TString option = GetOption();
 120   AliInfo(Form("Called with option %s.", option.Data()));
 121
 122   if (option.Contains("only-positive"))
 123   {
 124     AliInfo("Processing only positive particles.");
 125     fSignMode = 1;
 126   }
 127   else if (option.Contains("only-negative"))
 128   {
 129     AliInfo("Processing only negative particles.");
 130     fSignMode = -1;
 131   }
 132 }
 133
 134 void AlidNdEtaCorrectionSelector::SlaveBegin(TTree * tree)
 135 {
 136   // The SlaveBegin() function is called after the Begin() function.
 137   // When running with PROOF SlaveBegin() is called on each slave server.
 138   // The tree argument is deprecated (on PROOF 0 is passed).
 139
 140   AliSelectorRL::SlaveBegin(tree);
 141
 142   ReadUserObjects(tree);
 143
 144   fdNdEtaCorrection = new AlidNdEtaCorrection("dndeta_correction", "dndeta_correction");
 145
 146   fPIDParticles = new TH1F("pid_particles", "PID of generated primary particles", 10001, -5000.5, 5000.5);
 147   fPIDTracks = new TH1F("pid_tracks", "MC PID of reconstructed tracks", 10001, -5000.5, 5000.5);
 148
 149   fClustersITSPos = new TH1F("clusters_its_pos", "clusters_its_pos", 7, -0.5, 6.5);
 150   fClustersTPCPos = new TH1F("clusters_tpc_pos", "clusters_tpc_pos", 160, -0.5, 159.5);
 151
 152   fClustersITSNeg = new TH1F("clusters_its_neg", "clusters_its_neg", 7, -0.5, 6.5);
 153   fClustersTPCNeg = new TH1F("clusters_tpc_neg", "clusters_tpc_neg", 160, -0.5, 159.5);
 154
 155   fdNdEtaAnalysisMC = new dNdEtaAnalysis("dndetaMC", "dndetaMC");
 156   fdNdEtaAnalysisESD = new dNdEtaAnalysis("dndetaESD", "dndetaESD");
 157 }
 158
 159 void AlidNdEtaCorrectionSelector::Init(TTree* tree)
 160 {
 161   // read the user objects
 162
 163   AliSelectorRL::Init(tree);
 164
 165   // Enable only the needed branches
 166   if (tree)
 167   {
 168     tree->SetBranchStatus("*", 0);
 169     tree->SetBranchStatus("fTriggerMask", 1);
 170     tree->SetBranchStatus("fSPDVertex*", 1);
 171     tree->SetBranchStatus("fTracks.fLabel", 1);
 172     tree->SetBranchStatus("fTracks.fITSncls", 1);
 173     tree->SetBranchStatus("fTracks.fTPCncls", 1);
 174
 175     AliESDtrackCuts::EnableNeededBranches(tree);
 176   }
 177 }
 178
 179 Bool_t AlidNdEtaCorrectionSelector::Process(Long64_t entry)
 180 {
 181   // The Process() function is called for each entry in the tree (or possibly
 182   // keyed object in the case of PROOF) to be processed. The entry argument
 183   // specifies which entry in the currently loaded tree is to be processed.
 184   // It can be passed to either TTree::GetEntry() or TBranch::GetEntry()
 185   // to read either all or the required parts of the data. When processing
 186   // keyed objects with PROOF, the object is already loaded and is available
 187   // via the fObject pointer.
 188   //
 189   // This function should contain the "body" of the analysis. It can contain
 190   // simple or elaborate selection criteria, run algorithms on the data
 191   // of the event and typically fill histograms.
 192
 193   // WARNING when a selector is used with a TChain, you must use
 194   //  the pointer to the current TTree to call GetEntry(entry).
 195   //  The entry is always the local entry number in the current tree.
 196   //  Assuming that fTree is the pointer to the TChain being processed,
 197   //  use fTree->GetTree()->GetEntry(entry).
 198
 199   AliDebug(AliLog::kDebug+1,"Processing event ...\n");
 200
 201   if (AliSelectorRL::Process(entry) == kFALSE)
 202     return kFALSE;
 203
 204   // check prerequesites
 205   if (!fESD)
 206   {
 207     AliDebug(AliLog::kError, "ESD branch not available");
 208     return kFALSE;
 209   }
 210
 211   AliHeader* header = GetHeader();
 212   if (!header)
 213   {
 214     AliDebug(AliLog::kError, "Header not available");
 215     return kFALSE;
 216   }
 217
 218   AliStack* stack = GetStack();
 219   if (!stack)
 220   {
 221     AliDebug(AliLog::kError, "Stack not available");
 222     return kFALSE;
 223   }
 224
 225   if (!fEsdTrackCuts)
 226   {
 227     AliDebug(AliLog::kError, "fESDTrackCuts not available");
 228     return kFALSE;
 229   }
 230
 231   Bool_t vertexReconstructed = AliPWG0Helper::IsVertexReconstructed(fESD);
 232
 233   Bool_t eventTriggered = AliPWG0Helper::IsEventTriggered(fESD);
 234
 235   // get the MC vertex
 236   AliGenEventHeader* genHeader = header->GenEventHeader();
 237
 238   // primary vertex (from MC)
 239   TArrayF vtxMC(3);
 240   genHeader->PrimaryVertex(vtxMC);
 241
 242   // get process type
 243   Int_t processType = AliPWG0Helper::GetPythiaEventProcessType(header);
 244   AliDebug(AliLog::kDebug+1, Form("Found pythia procces type %d", processType));
 245
 246   if (processType<0)
 247     AliDebug(AliLog::kError, Form("Unknown Pythia process type %d.", processType));
 248
 249   // loop over mc particles
 250   Int_t nPrim  = stack->GetNprimary();
 251
 252   for (Int_t iMc = 0; iMc < nPrim; ++iMc)
 253   {
 254     AliDebug(AliLog::kDebug+1, Form("MC Loop: Processing particle %d.", iMc));
 255
 256     TParticle* particle = stack->Particle(iMc);
 257
 258     if (!particle)
 259     {
 260       AliDebug(AliLog::kError, Form("UNEXPECTED: particle with label %d not found in stack (mc loop).", iMc));
 261       continue;
 262     }
 263
 264     if (AliPWG0Helper::IsPrimaryCharged(particle, nPrim) == kFALSE)
 265       continue;
 266
 267     if (SignOK(particle->GetPDG()) == kFALSE)
 268       continue;
 269
 270     Float_t eta = particle->Eta();
 271     Float_t pt = particle->Pt();
 272
 273     fdNdEtaCorrection->FillMCParticle(vtxMC[2], eta, pt, eventTriggered, vertexReconstructed, processType);
 274
 275     if (eventTriggered)
 276     {
 277       if (vertexReconstructed)
 278       {
 279         fdNdEtaAnalysisMC->FillTrack(vtxMC[2], eta, pt);
 280       }
 281     }
 282   } // end of mc particle
 283
 284   // ########################################################
 285   // loop over esd tracks
 286   Int_t nTracks = fESD->GetNumberOfTracks();
 287
 288   // count the number of "good" tracks as parameter for vertex reconstruction efficiency
 289   Int_t nGoodTracks = 0;
 290   for (Int_t t=0; t<nTracks; t++)
 291   {
 292     AliDebug(AliLog::kDebug+1, Form("ESD Loop: Processing track %d.", t));
 293
 294     AliESDtrack* esdTrack = fESD->GetTrack(t);
 295
 296     // cut the esd track?
 297     if (!fEsdTrackCuts->AcceptTrack(esdTrack))
 298       continue;
 299
 300     nGoodTracks++;
 301
 302     // using the properties of the mc particle
 303     Int_t label = TMath::Abs(esdTrack->GetLabel());
 304     if (label == 0)
 305     {
 306       AliDebug(AliLog::kWarning, Form("WARNING: cannot find corresponding mc part for track %d.", t));
 307       continue;
 308     }
 309
 310     TParticle* particle = stack->Particle(label);
 311     if (!particle)
 312     {
 313       AliDebug(AliLog::kError, Form("UNEXPECTED: particle with label %d not found in stack (track loop).", label));
 314       continue;
 315     }
 316
 317     if (SignOK(particle->GetPDG()) == kFALSE)
 318         continue;
 319
 320     if (eventTriggered && vertexReconstructed)
 321     {
 322       fdNdEtaCorrection->FillTrackedParticle(vtxMC[2], particle->Eta(), particle->Pt());
 323       fdNdEtaAnalysisESD->FillTrack(vtxMC[2], particle->Eta(), particle->Pt());
 324       if (particle->Pt() > 0.1 && particle->Pt() < 0.2)
 325       {
 326         fPIDTracks->Fill(particle->GetPdgCode());
 327         if (particle->GetPDG()->Charge() > 0)
 328         {
 329           fClustersITSPos->Fill(esdTrack->GetITSclusters(0));
 330           fClustersTPCPos->Fill(esdTrack->GetTPCclusters(0));
 331         }
 332         else
 333         {
 334           fClustersITSNeg->Fill(esdTrack->GetITSclusters(0));
 335           fClustersTPCNeg->Fill(esdTrack->GetTPCclusters(0));
 336         }
 337       }
 338     }
 339   } // end of track loop
 340
 341   if (eventTriggered && vertexReconstructed)
 342     fdNdEtaAnalysisMC->FillEvent(vtxMC[2], nGoodTracks);
 343
 344   // stuff regarding the vertex reco correction and trigger bias correction
 345   fdNdEtaCorrection->FillEvent(vtxMC[2], nGoodTracks, eventTriggered, vertexReconstructed, processType);
 346   if (eventTriggered) {
 347     if (vertexReconstructed)
 348     {
 349       fdNdEtaAnalysisESD->FillEvent(vtxMC[2], nGoodTracks);
 350     }
 351   }
 352
 353   return kTRUE;
 354 }
 355
 356 void AlidNdEtaCorrectionSelector::SlaveTerminate()
 357 {
 358   // The SlaveTerminate() function is called after all entries or objects
 359   // have been processed. When running with PROOF SlaveTerminate() is called
 360   // on each slave server.
 361
 362   AliSelectorRL::SlaveTerminate();
 363
 364   // Add the histograms to the output on each slave server
 365   if (!fOutput)
 366   {
 367     AliDebug(AliLog::kError, "ERROR: Output list not initialized");
 368     return;
 369   }
 370
 371   fOutput->Add(fdNdEtaCorrection);
 372   fOutput->Add(fdNdEtaAnalysisMC);
 373   fOutput->Add(fdNdEtaAnalysisESD);
 374 }
 375
 376 void AlidNdEtaCorrectionSelector::Terminate()
 377 {
 378   // The Terminate() function is the last function to be called during
 379   // a query. It always runs on the client, it can be used to present
 380   // the results graphically or save the results to file.
 381
 382   AliSelectorRL::Terminate();
 383
 384   fdNdEtaCorrection = dynamic_cast<AlidNdEtaCorrection*> (fOutput->FindObject("dndeta_correction"));
 385   fdNdEtaAnalysisMC = dynamic_cast<dNdEtaAnalysis*> (fOutput->FindObject("dndetaMC"));
 386   fdNdEtaAnalysisESD = dynamic_cast<dNdEtaAnalysis*> (fOutput->FindObject("dndetaESD"));
 387   if (!fdNdEtaCorrection || !fdNdEtaAnalysisMC || !fdNdEtaAnalysisESD)
 388   {
 389     AliDebug(AliLog::kError, "Could not read object from output list");
 390     return;
 391   }
 392
 393   fdNdEtaCorrection->Finish();
 394
 395   TFile* fout = new TFile(Form("correction_map%s.root", GetOption()), "RECREATE");
 396
 397   if (fEsdTrackCuts)
 398     fEsdTrackCuts->SaveHistograms("esd_track_cuts");
 399   fdNdEtaCorrection->SaveHistograms();
 400   fdNdEtaAnalysisMC->SaveHistograms();
 401   fdNdEtaAnalysisESD->SaveHistograms();
 402
 403   fout->Write();
 404   fout->Close();
 405
 406   fdNdEtaCorrection->DrawHistograms();
 407
 408   if (fPIDParticles && fPIDTracks)
 409   {
 410     new TCanvas("pidcanvas", "pidcanvas", 500, 500);
 411
 412     fPIDParticles->Draw();
 413     fPIDTracks->SetLineColor(2);
 414     fPIDTracks->Draw("SAME");
 415
 416     TDatabasePDG* pdgDB = new TDatabasePDG;
 417
 418     for (Int_t i=0; i <= fPIDParticles->GetNbinsX()+1; ++i)
 419       if (fPIDParticles->GetBinContent(i) > 0)
 420         printf("PDG = %d (%s): generated: %d, reconstructed: %d, ratio: %f\n", (Int_t) fPIDParticles->GetBinCenter(i), pdgDB->GetParticle((Int_t) fPIDParticles->GetBinCenter(i))->GetName(), (Int_t) fPIDParticles->GetBinContent(i), (Int_t) fPIDTracks->GetBinContent(i), ((fPIDTracks->GetBinContent(i) > 0) ? fPIDParticles->GetBinContent(i) / fPIDTracks->GetBinContent(i) : -1));
 421
 422     delete pdgDB;
 423     pdgDB = 0;
 424   }
 425
 426   if (fClustersITSPos && fClustersITSNeg && fClustersTPCPos && fClustersTPCNeg)
 427   {
 428     TCanvas* canvas = new TCanvas("clusters", "clusters", 1000, 500);
 429     canvas->Divide(2, 1);
 430
 431     canvas->cd(1);
 432     fClustersITSPos->Draw();
 433     fClustersITSNeg->SetLineColor(kRed);
 434     fClustersITSNeg->Draw("SAME");
 435
 436     canvas->cd(2);
 437     fClustersTPCPos->Draw();
 438     fClustersTPCNeg->SetLineColor(kRed);
 439     fClustersTPCNeg->Draw("SAME");
 440   }
 441 }