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