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