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