PWG0/dNdEta/AlidNdEtaSystematicsSelector.cxx

   1 /* $Id$ */
   2
   3 #include "AlidNdEtaSystematicsSelector.h"
   4
   5 #include <TStyle.h>
   6 #include <TSystem.h>
   7 #include <TCanvas.h>
   8 #include <TVector3.h>
   9 #include <TChain.h>
  10 #include <TFile.h>
  11 #include <TH2F.h>
  12 #include <TH1F.h>
  13 #include <TParticle.h>
  14
  15 #include <AliLog.h>
  16 #include <AliESD.h>
  17 #include <AliStack.h>
  18 #include <AliHeader.h>
  19 #include <AliGenEventHeader.h>
  20 #include <AliGenPythiaEventHeader.h>
  21 #include <AliGenCocktailEventHeader.h>
  22
  23
  24 #include "esdTrackCuts/AliESDtrackCuts.h"
  25 #include "AliPWG0Helper.h"
  26 #include "dNdEta/AlidNdEtaCorrection.h"
  27
  28 ClassImp(AlidNdEtaSystematicsSelector)
  29
  30 AlidNdEtaSystematicsSelector::AlidNdEtaSystematicsSelector() :
  31   AliSelectorRL(),
  32   fSecondaries(0),
  33   fSigmaVertex(0),
  34   fEsdTrackCuts(0),
  35   fPIDParticles(0),
  36   fPIDTracks(0)
  37 {
  38   //
  39   // Constructor. Initialization of pointers
  40   //
  41
  42   for (Int_t i=0; i<4; ++i)
  43     fdNdEtaCorrectionSpecies[i] = 0;
  44
  45   for (Int_t i=0; i<3; ++i) {
  46     fdNdEtaCorrectionVertexReco[i] = 0;
  47     fdNdEtaCorrectionTriggerBias[i] = 0;
  48   }
  49 }
  50
  51 AlidNdEtaSystematicsSelector::~AlidNdEtaSystematicsSelector()
  52 {
  53   //
  54   // Destructor
  55   //
  56 }
  57
  58 void AlidNdEtaSystematicsSelector::Begin(TTree* tree)
  59 {
  60   // Begin function
  61
  62   AliSelectorRL::Begin(tree);
  63
  64   ReadUserObjects(tree);
  65 }
  66
  67 void AlidNdEtaSystematicsSelector::ReadUserObjects(TTree* tree)
  68 {
  69   // read the user objects, called from slavebegin and begin
  70
  71   if (!fEsdTrackCuts && fInput)
  72     fEsdTrackCuts = dynamic_cast<AliESDtrackCuts*> (fInput->FindObject("AliESDtrackCuts"));
  73
  74   if (!fEsdTrackCuts && tree)
  75     fEsdTrackCuts = dynamic_cast<AliESDtrackCuts*> (tree->GetUserInfo()->FindObject("AliESDtrackCuts"));
  76
  77   if (!fEsdTrackCuts)
  78      AliDebug(AliLog::kError, "ERROR: Could not read EsdTrackCuts from input list.");
  79 }
  80
  81 void AlidNdEtaSystematicsSelector::SlaveBegin(TTree* tree)
  82 {
  83   // The SlaveBegin() function is called after the Begin() function.
  84   // When running with PROOF SlaveBegin() is called on each slave server.
  85   // The tree argument is deprecated (on PROOF 0 is passed).
  86
  87   AliSelector::SlaveBegin(tree);
  88
  89   ReadUserObjects(tree);
  90
  91   TString option(GetOption());
  92
  93   printf("Running AlidNdEtaSystematicsSelector with options %s\n", option.Data());
  94
  95   if (option.Contains("secondaries"))
  96   {
  97     fSecondaries = new TH2F("fSecondaries", "fSecondaries;Case;N", 8, 0.5, 8.5, 2001, -0.5, 2000.5);
  98     fSecondaries->GetXaxis()->SetBinLabel(1, "All Primaries");
  99     fSecondaries->GetXaxis()->SetBinLabel(2, "All Secondaries");
 100     fSecondaries->GetXaxis()->SetBinLabel(3, "Primaries pT > 0.3 GeV/c");
 101     fSecondaries->GetXaxis()->SetBinLabel(4, "Secondaries pT > 0.3 GeV/c");
 102     fSecondaries->GetXaxis()->SetBinLabel(5, "Tracks from Primaries");
 103     fSecondaries->GetXaxis()->SetBinLabel(6, "Tracks from Secondaries");
 104     fSecondaries->GetXaxis()->SetBinLabel(7, "Accepted Tracks from Primaries");
 105     fSecondaries->GetXaxis()->SetBinLabel(8, "Accepted Tracks from Secondaries");
 106   }
 107
 108   if (option.Contains("particle-composition"))
 109   {
 110     for (Int_t i=0; i<4; ++i)
 111     {
 112       TString name;
 113       name.Form("correction_%d", i);
 114       fdNdEtaCorrectionSpecies[i] = new AlidNdEtaCorrection(name, name);
 115     }
 116   }
 117
 118   if (option.Contains("sigma-vertex"))
 119   {
 120     fSigmaVertex = new TH1F("fSigmaVertex", "fSigmaVertex;Nsigma2vertex;NacceptedTracks", 50, 0.05, 5.05);
 121     printf("WARNING: sigma-vertex analysis enabled. This will produce weird results in the AliESDtrackCuts histograms\n");
 122   }
 123
 124   if (option.Contains("vertexreco")) {
 125     fdNdEtaCorrectionVertexReco[0] = new AlidNdEtaCorrection("vertexRecoND", "vertexRecoND");
 126     fdNdEtaCorrectionVertexReco[1] = new AlidNdEtaCorrection("vertexRecoSD", "vertexRecoSD");
 127     fdNdEtaCorrectionVertexReco[2] = new AlidNdEtaCorrection("vertexRecoDD", "vertexRecoDD");
 128   }
 129
 130   if (option.Contains("triggerbias")) {
 131     fdNdEtaCorrectionTriggerBias[0] = new AlidNdEtaCorrection("triggerBiasND", "triggerBiasND");
 132     fdNdEtaCorrectionTriggerBias[1] = new AlidNdEtaCorrection("triggerBiasSD", "triggerBiasSD");
 133     fdNdEtaCorrectionTriggerBias[2] = new AlidNdEtaCorrection("triggerBiasDD", "triggerBiasDD");
 134   }
 135
 136
 137   fPIDParticles = new TH1F("pid_particles", "PID of generated primary particles", 10001, -5000.5, 5000.5);
 138   fPIDTracks = new TH1F("pid_tracks", "MC PID of reconstructed tracks", 10001, -5000.5, 5000.5);
 139 }
 140
 141 Bool_t AlidNdEtaSystematicsSelector::Process(Long64_t entry)
 142 {
 143   // The Process() function is called for each entry in the tree (or possibly
 144   // keyed object in the case of PROOF) to be processed. The entry argument
 145   // specifies which entry in the currently loaded tree is to be processed.
 146   // It can be passed to either TTree::GetEntry() or TBranch::GetEntry()
 147   // to read either all or the required parts of the data. When processing
 148   // keyed objects with PROOF, the object is already loaded and is available
 149   // via the fObject pointer.
 150   //
 151   // This function should contain the "body" of the analysis. It can contain
 152   // simple or elaborate selection criteria, run algorithms on the data
 153   // of the event and typically fill histograms.
 154
 155   // WARNING when a selector is used with a TChain, you must use
 156   //  the pointer to the current TTree to call GetEntry(entry).
 157   //  The entry is always the local entry number in the current tree.
 158   //  Assuming that fTree is the pointer to the TChain being processed,
 159   //  use fTree->GetTree()->GetEntry(entry).
 160
 161   if (AliSelectorRL::Process(entry) == kFALSE)
 162     return kFALSE;
 163
 164   // Check prerequisites
 165   if (!fESD)
 166   {
 167     AliDebug(AliLog::kError, "ESD branch not available");
 168     return kFALSE;
 169   }
 170
 171   if (!fEsdTrackCuts)
 172   {
 173     AliDebug(AliLog::kError, "fESDTrackCuts not available");
 174     return kFALSE;
 175   }
 176
 177   AliStack* stack = GetStack();
 178   if (!stack)
 179   {
 180     AliDebug(AliLog::kError, "Stack not available");
 181     return kFALSE;
 182   }
 183
 184   TObjArray* list = fEsdTrackCuts->GetAcceptedTracks(fESD);
 185
 186   if (fdNdEtaCorrectionSpecies[0])
 187     FillCorrectionMaps(list);
 188
 189   if (fSecondaries)
 190     FillSecondaries();
 191
 192   if (fSigmaVertex)
 193     FillSigmaVertex();
 194
 195   // stuff for vertex reconstruction correction systematics
 196   Bool_t vertexRecoStudy  = kFALSE;
 197   Bool_t triggerBiasStudy = kFALSE;
 198   if (fdNdEtaCorrectionVertexReco[0]) vertexRecoStudy = kTRUE;
 199   if (fdNdEtaCorrectionTriggerBias[0]) triggerBiasStudy = kTRUE;
 200
 201   if (vertexRecoStudy || triggerBiasStudy) {
 202     AliHeader* header = GetHeader();
 203     if (!header) {
 204       AliDebug(AliLog::kError, "Header not available");
 205       return kFALSE;
 206     }
 207
 208     // getting process information
 209     Int_t processtype = AliPWG0Helper::GetPythiaEventProcessType(header);
 210
 211     AliDebug(AliLog::kInfo, Form("Pythia process type %d.",processtype));
 212
 213     // can only read pythia headers, either directly or from cocktalil header
 214     AliGenEventHeader* genHeader = (AliGenEventHeader*)(header->GenEventHeader());
 215
 216     if (!genHeader) {
 217       AliDebug(AliLog::kError, "Gen header not available");
 218       return kFALSE;
 219     }
 220
 221     // get the MC vertex
 222     TArrayF vtxMC(3);
 223     genHeader->PrimaryVertex(vtxMC);
 224
 225     Int_t nGoodTracks = list->GetEntries();
 226
 227     Bool_t eventTriggered = AliPWG0Helper::IsEventTriggered(fESD);
 228     Bool_t vertexReconstructed = AliPWG0Helper::IsVertexReconstructed(fESD);
 229
 230     // non diffractive
 231     if (processtype!=92 && processtype!=93 && processtype!=94) {
 232       if (triggerBiasStudy) fdNdEtaCorrectionTriggerBias[0]->FillEventAll(vtxMC[2], nGoodTracks);
 233
 234       if (eventTriggered) {
 235         if (triggerBiasStudy) fdNdEtaCorrectionTriggerBias[0]->FillEventWithTrigger(vtxMC[2], nGoodTracks);
 236         if (vertexRecoStudy) fdNdEtaCorrectionVertexReco[0]->FillEventWithTrigger(vtxMC[2], nGoodTracks);
 237
 238         if (vertexReconstructed)
 239           if (vertexRecoStudy) fdNdEtaCorrectionVertexReco[0]->FillEventWithTriggerWithReconstructedVertex(vtxMC[2], nGoodTracks);
 240       }
 241     }
 242
 243     // single diffractive
 244     if (processtype==92 || processtype==93) {
 245       if (triggerBiasStudy) fdNdEtaCorrectionTriggerBias[1]->FillEventAll(vtxMC[2], nGoodTracks);
 246
 247       if (eventTriggered) {
 248         if (triggerBiasStudy) fdNdEtaCorrectionTriggerBias[1]->FillEventWithTrigger(vtxMC[2], nGoodTracks);
 249         if (vertexRecoStudy) fdNdEtaCorrectionVertexReco[1]->FillEventWithTrigger(vtxMC[2], nGoodTracks);
 250
 251         if (vertexReconstructed)
 252           if (vertexRecoStudy) fdNdEtaCorrectionVertexReco[1]->FillEventWithTriggerWithReconstructedVertex(vtxMC[2], nGoodTracks);
 253       }
 254     }
 255
 256     // double diffractive
 257     if (processtype==94) {
 258       if (triggerBiasStudy) fdNdEtaCorrectionTriggerBias[2]->FillEventAll(vtxMC[2], nGoodTracks);
 259
 260       if (eventTriggered) {
 261         if (triggerBiasStudy) fdNdEtaCorrectionTriggerBias[2]->FillEventWithTrigger(vtxMC[2], nGoodTracks);
 262         if (vertexRecoStudy) fdNdEtaCorrectionVertexReco[2]->FillEventWithTrigger(vtxMC[2], nGoodTracks);
 263
 264         if (vertexReconstructed)
 265           if (vertexRecoStudy) fdNdEtaCorrectionVertexReco[2]->FillEventWithTriggerWithReconstructedVertex(vtxMC[2], nGoodTracks);
 266       }
 267     }
 268   }
 269
 270   delete list;
 271   list = 0;
 272
 273   return kTRUE;
 274 }
 275
 276 void AlidNdEtaSystematicsSelector::FillCorrectionMaps(TObjArray* listOfTracks)
 277 {
 278   // fills the correction maps for different particle species
 279
 280   AliStack* stack = GetStack();
 281   AliHeader* header = GetHeader();
 282
 283   Bool_t vertexReconstructed = AliPWG0Helper::IsVertexReconstructed(fESD);
 284   Bool_t eventTriggered = AliPWG0Helper::IsEventTriggered(fESD);
 285
 286   // get the MC vertex
 287   AliGenEventHeader* genHeader = header->GenEventHeader();
 288
 289   TArrayF vtxMC(3);
 290   genHeader->PrimaryVertex(vtxMC);
 291
 292   // loop over mc particles
 293   Int_t nPrim  = stack->GetNprimary();
 294
 295   for (Int_t iMc = 0; iMc < nPrim; ++iMc)
 296   {
 297     TParticle* particle = stack->Particle(iMc);
 298
 299     if (!particle)
 300     {
 301       AliDebug(AliLog::kError, Form("UNEXPECTED: particle with label %d not found in stack (mc loop).", iMc));
 302       continue;
 303     }
 304
 305     if (AliPWG0Helper::IsPrimaryCharged(particle, nPrim) == kFALSE)
 306       continue;
 307
 308     Float_t eta = particle->Eta();
 309     Float_t pt = particle->Pt();
 310
 311     Int_t id = -1;
 312     switch (TMath::Abs(particle->GetPdgCode()))
 313     {
 314       case 211: id = 0; break;
 315       case 321: id = 1; break;
 316       case 2212: id = 2; break;
 317       case 11:
 318         {
 319           /*if (pt < 0.1 && particle->GetMother(0) > -1)
 320           {
 321             TParticle* mother = stack->Particle(particle->GetMother(0));
 322             printf("Mother pdg code is %d\n", mother->GetPdgCode());
 323           } */
 324           //particle->Dump();
 325           //if (particle->GetUniqueID() == 1)
 326
 327           //printf("Found illegal particle mc id = %d file = %s event = %d\n", iMc,  fTree->GetCurrentFile()->GetName(), fTree->GetTree()->GetReadEntry());
 328         }
 329       default: id = 3; break;
 330     }
 331
 332     if (vertexReconstructed)
 333     {
 334       fdNdEtaCorrectionSpecies[id]->FillParticle(vtxMC[2], eta, pt);
 335       //if (pt < 0.1)
 336         fPIDParticles->Fill(particle->GetPdgCode());
 337     }
 338
 339     //fdNdEtaCorrectionSpecies[id]->FillParticleAllEvents(eta, pt);
 340     //if (eventTriggered)
 341     // fdNdEtaCorrectionSpecies[id]->FillParticleWhenEventTriggered(eta, pt);
 342   }// end of mc particle
 343
 344   // loop over esd tracks
 345   TIterator* iter = listOfTracks->MakeIterator();
 346   TObject* obj = 0;
 347   while ((obj = iter->Next()) != 0)
 348   {
 349     AliESDtrack* esdTrack = dynamic_cast<AliESDtrack*> (obj);
 350     if (!esdTrack)
 351       continue;
 352
 353     // using the properties of the mc particle
 354     Int_t label = TMath::Abs(esdTrack->GetLabel());
 355     TParticle* particle = stack->Particle(label);
 356     if (!particle)
 357     {
 358       AliDebug(AliLog::kError, Form("UNEXPECTED: particle with label %d not found in stack (track loop).", label));
 359       continue;
 360     }
 361
 362     TParticle* mother = particle;
 363     // find primary particle that created this particle
 364     while (label >= nPrim)
 365     {
 366       //printf("Particle %d (pdg %d) is not a primary. Let's check its mother %d\n", label, mother->GetPdgCode(), mother->GetMother(0));
 367
 368       if (mother->GetMother(0) == -1)
 369       {
 370         AliDebug(AliLog::kError, Form("UNEXPECTED: Could not find mother of secondary particle %d.", label));
 371         mother = 0;
 372         break;
 373       }
 374
 375       label = mother->GetMother(0);
 376
 377       mother = stack->Particle(label);
 378       if (!mother)
 379       {
 380         AliDebug(AliLog::kError, Form("UNEXPECTED: particle with label %d not found in stack (find mother loop).", label));
 381         break;
 382       }
 383     }
 384
 385     if (!mother)
 386       continue;
 387
 388     //printf("We continue with particle %d (pdg %d)\n", label, mother->GetPdgCode());
 389
 390     Int_t id = -1;
 391     switch (TMath::Abs(mother->GetPdgCode()))
 392     {
 393       case 211:  id = 0; break;
 394       case 321:  id = 1; break;
 395       case 2212: id = 2; break;
 396       default:   id = 3; break;
 397     }
 398
 399     if (vertexReconstructed)
 400     {
 401       fdNdEtaCorrectionSpecies[id]->FillParticleWhenMeasuredTrack(vtxMC[2], particle->Eta(), particle->Pt());
 402       //if (particle->Pt() < 0.1)
 403         fPIDTracks->Fill(particle->GetPdgCode());
 404     }
 405   } // end of track loop
 406
 407   Int_t nGoodTracks = listOfTracks->GetEntries();
 408
 409   for (Int_t i=0; i<4; ++i)
 410   {
 411     fdNdEtaCorrectionSpecies[i]->FillEventAll(vtxMC[2], nGoodTracks);
 412     if (eventTriggered)
 413     {
 414       fdNdEtaCorrectionSpecies[i]->FillEventWithTrigger(vtxMC[2], nGoodTracks);
 415       if (vertexReconstructed)
 416         fdNdEtaCorrectionSpecies[i]->FillEventWithTriggerWithReconstructedVertex(vtxMC[2], nGoodTracks);
 417     }
 418   }
 419
 420   delete iter;
 421   iter = 0;
 422 }
 423
 424 void AlidNdEtaSystematicsSelector::FillSecondaries()
 425 {
 426   // fills the secondary histograms
 427
 428   AliStack* stack = GetStack();
 429
 430   Int_t particlesPrimaries = 0;
 431   Int_t particlesSecondaries = 0;
 432   Int_t particlesPrimariesPtCut = 0;
 433   Int_t particlesSecondariesPtCut = 0;
 434
 435   for (Int_t iParticle = 0; iParticle < stack->GetNtrack(); iParticle++)
 436   {
 437     TParticle* particle = stack->Particle(iParticle);
 438     if (!particle)
 439     {
 440       AliDebug(AliLog::kError, Form("UNEXPECTED: particle with label %d not found in stack (particle loop).", iParticle));
 441       continue;
 442     }
 443
 444     if (TMath::Abs(particle->Eta()) > 0.9)
 445       continue;
 446
 447     Bool_t isPrimary = kFALSE;
 448     if (iParticle < stack->GetNprimary())
 449     {
 450       if (AliPWG0Helper::IsPrimaryCharged(particle, stack->GetNprimary()) == kFALSE)
 451         continue;
 452
 453       isPrimary = kTRUE;
 454     }
 455
 456     if (isPrimary)
 457       particlesPrimaries++;
 458     else
 459       particlesSecondaries++;
 460
 461     if (particle->Pt() < 0.3)
 462       continue;
 463
 464     if (isPrimary)
 465       particlesPrimariesPtCut++;
 466     else
 467       particlesSecondariesPtCut++;
 468   }
 469
 470   fSecondaries->Fill(1, particlesPrimaries);
 471   fSecondaries->Fill(2, particlesSecondaries);
 472   fSecondaries->Fill(3, particlesPrimariesPtCut);
 473   fSecondaries->Fill(4, particlesSecondariesPtCut);
 474
 475   Int_t allTracksPrimaries = 0;
 476   Int_t allTracksSecondaries = 0;
 477   Int_t acceptedTracksPrimaries = 0;
 478   Int_t acceptedTracksSecondaries = 0;
 479
 480   for (Int_t iTrack = 0; iTrack < fESD->GetNumberOfTracks(); iTrack++)
 481   {
 482     AliESDtrack* esdTrack = fESD->GetTrack(iTrack);
 483     if (!esdTrack)
 484     {
 485       AliDebug(AliLog::kError, Form("UNEXPECTED: Could not get track %d.", iTrack));
 486       continue;
 487     }
 488
 489     Int_t label = TMath::Abs(esdTrack->GetLabel());
 490     TParticle* particle = stack->Particle(label);
 491     if (!particle)
 492     {
 493       AliDebug(AliLog::kError, Form("UNEXPECTED: particle with label %d not found in stack (track loop).", label));
 494       continue;
 495     }
 496
 497    if (label < stack->GetNprimary())
 498       allTracksPrimaries++;
 499     else
 500       allTracksSecondaries++;
 501
 502     if (!fEsdTrackCuts->AcceptTrack(esdTrack))
 503       continue;
 504
 505     if (label < stack->GetNprimary())
 506       acceptedTracksPrimaries++;
 507     else
 508       acceptedTracksSecondaries++;
 509   }
 510
 511   fSecondaries->Fill(5, allTracksPrimaries);
 512   fSecondaries->Fill(6, allTracksSecondaries);
 513   fSecondaries->Fill(7, acceptedTracksPrimaries);
 514   fSecondaries->Fill(8, acceptedTracksSecondaries);
 515
 516   //printf("P = %d, S = %d, P_pt = %d, S_pt = %d, T_P = %d, T_S = %d, T_P_acc = %d, T_S_acc = %d\n", particlesPrimaries, particlesSecondaries, particlesPrimariesPtCut, particlesSecondariesPtCut, allTracksPrimaries, allTracksSecondaries, acceptedTracksPrimaries, acceptedTracksSecondaries);
 517 }
 518
 519 void AlidNdEtaSystematicsSelector::FillSigmaVertex()
 520 {
 521   // fills the fSigmaVertex histogram
 522
 523   // save the old value
 524   Float_t oldSigmaVertex = fEsdTrackCuts->GetMinNsigmaToVertex();
 525
 526   // set to maximum
 527   fEsdTrackCuts->SetMinNsigmaToVertex(5);
 528
 529   TObjArray* list = fEsdTrackCuts->GetAcceptedTracks(fESD);
 530
 531   TIterator* iter = list->MakeIterator();
 532   TObject* obj = 0;
 533   while ((obj = iter->Next()) != 0)
 534   {
 535     AliESDtrack* esdTrack = dynamic_cast<AliESDtrack*> (obj);
 536     if (!esdTrack)
 537       continue;
 538
 539     Float_t sigma2Vertex = fEsdTrackCuts->GetSigmaToVertex(esdTrack);
 540
 541     for (Double_t nSigma = 0.1; nSigma < 5.05; nSigma += 0.1)
 542     {
 543       if (sigma2Vertex < nSigma)
 544         fSigmaVertex->Fill(nSigma);
 545     }
 546   }
 547
 548   delete iter;
 549   iter = 0;
 550
 551   delete list;
 552   list = 0;
 553
 554   // set back the old value
 555   fEsdTrackCuts->SetMinNsigmaToVertex(oldSigmaVertex);
 556 }
 557
 558 void AlidNdEtaSystematicsSelector::SlaveTerminate()
 559 {
 560   // The SlaveTerminate() function is called after all entries or objects
 561   // have been processed. When running with PROOF SlaveTerminate() is called
 562   // on each slave server.
 563
 564   AliSelector::SlaveTerminate();
 565
 566   // Add the histograms to the output on each slave server
 567   if (!fOutput)
 568   {
 569     AliDebug(AliLog::kError, Form("ERROR: Output list not initialized."));
 570     return;
 571   }
 572
 573   if (fSecondaries)
 574     fOutput->Add(fSecondaries);
 575
 576   for (Int_t i=0; i<4; ++i)
 577     if (fdNdEtaCorrectionSpecies[i])
 578       fOutput->Add(fdNdEtaCorrectionSpecies[i]);
 579
 580   if (fSigmaVertex)
 581     fOutput->Add(fSigmaVertex);
 582
 583   for (Int_t i=0; i<3; ++i) {
 584     if (fdNdEtaCorrectionVertexReco[i])
 585       fOutput->Add(fdNdEtaCorrectionVertexReco[i]);
 586
 587     if (fdNdEtaCorrectionTriggerBias[i])
 588       fOutput->Add(fdNdEtaCorrectionTriggerBias[i]);
 589   }
 590
 591 }
 592
 593 void AlidNdEtaSystematicsSelector::Terminate()
 594 {
 595   // The Terminate() function is the last function to be called during
 596   // a query. It always runs on the client, it can be used to present
 597   // the results graphically or save the results to file.
 598
 599   AliSelector::Terminate();
 600
 601   fSecondaries = dynamic_cast<TH2F*> (fOutput->FindObject("fSecondaries"));
 602   for (Int_t i=0; i<4; ++i)
 603     fdNdEtaCorrectionSpecies[i] = dynamic_cast<AlidNdEtaCorrection*> (fOutput->FindObject(Form("correction_%d", i)));
 604   fSigmaVertex = dynamic_cast<TH1F*> (fOutput->FindObject("fSigmaVertex"));
 605
 606   if (fPIDParticles)
 607   {
 608     TDatabasePDG* pdgDB = new TDatabasePDG;
 609
 610     for (Int_t i=0; i <= fPIDParticles->GetNbinsX()+1; ++i)
 611       if (fPIDParticles->GetBinContent(i) > 0)
 612         printf("PDG = %d (%s): generated: %d, reconstructed: %d, ratio: %f\n",
 613                (Int_t) fPIDParticles->GetBinCenter(i),
 614                pdgDB->GetParticle((Int_t) fPIDParticles->GetBinCenter(i))->GetName(),
 615                (Int_t) fPIDParticles->GetBinContent(i),
 616                (Int_t) fPIDTracks->GetBinContent(i),
 617                ((fPIDTracks->GetBinContent(i) > 0) ? fPIDParticles->GetBinContent(i) / fPIDTracks->GetBinContent(i) : -1));
 618
 619     delete pdgDB;
 620     pdgDB = 0;
 621   }
 622
 623   for (Int_t i=0; i<3; ++i) {
 624     if (fdNdEtaCorrectionVertexReco[i])
 625       fdNdEtaCorrectionVertexReco[i]->Finish();
 626
 627     if (fdNdEtaCorrectionTriggerBias[i])
 628       fdNdEtaCorrectionTriggerBias[i]->Finish();
 629   }
 630
 631   TFile* fout = TFile::Open("systematics.root", "RECREATE");
 632
 633   if (fEsdTrackCuts)
 634     fEsdTrackCuts->SaveHistograms("esd_track_cuts");
 635
 636   if (fSecondaries)
 637     fSecondaries->Write();
 638
 639   if (fSigmaVertex)
 640     fSigmaVertex->Write();
 641
 642   for (Int_t i=0; i<4; ++i)
 643     if (fdNdEtaCorrectionSpecies[i])
 644       fdNdEtaCorrectionSpecies[i]->SaveHistograms();
 645
 646   for (Int_t i=0; i<3; ++i) {
 647     if (fdNdEtaCorrectionVertexReco[i])
 648       fdNdEtaCorrectionVertexReco[i]->SaveHistograms();
 649
 650     if (fdNdEtaCorrectionTriggerBias[i])
 651       fdNdEtaCorrectionTriggerBias[i]->SaveHistograms();
 652   }
 653
 654   fout->Write();
 655   fout->Close();
 656
 657   if (fSecondaries)
 658   {
 659     new TCanvas;
 660     fSecondaries->Draw("COLZ");
 661   }
 662
 663   if (fSigmaVertex)
 664   {
 665     new TCanvas;
 666     fSigmaVertex->Draw();
 667   }
 668 }