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