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
 395     TParticle* mother = particle;
 396     // find primary particle that created this particle
 397     while (label >= nPrim)
 398     {
 399       //printf("Particle %d (pdg %d) is not a primary. Let's check its mother %d\n", label, mother->GetPdgCode(), mother->GetMother(0));
 400
 401       if (mother->GetMother(0) == -1)
 402       {
 403         AliDebug(AliLog::kError, Form("UNEXPECTED: Could not find mother of secondary particle %d.", label));
 404         mother = 0;
 405         break;
 406       }
 407
 408       label = mother->GetMother(0);
 409
 410       mother = stack->Particle(label);
 411       if (!mother)
 412       {
 413         AliDebug(AliLog::kError, Form("UNEXPECTED: particle with label %d not found in stack (find mother loop).", label));
 414         break;
 415       }
 416     }
 417
 418     if (!mother)
 419       continue;
 420
 421     if (SignOK(mother->GetPDG()) == kFALSE)
 422       continue;
 423
 424     //printf("We continue with particle %d (pdg %d)\n", label, mother->GetPdgCode());
 425
 426     Int_t id = -1;
 427     switch (TMath::Abs(mother->GetPdgCode()))
 428     {
 429       case 211:  id = 0; break;
 430       case 321:  id = 1; break;
 431       case 2212: id = 2; break;
 432       default:   id = 3; break;
 433     }
 434
 435     if (vertexReconstructed && eventTriggered) {
 436       if (fdNdEtaCorrectionSpecies[id])
 437         fdNdEtaCorrectionSpecies[id]->FillTrackedParticle(vtxMC[2], eta, pt);
 438       //if (pt < 0.1)
 439       if (fPIDTracks)
 440         fPIDTracks->Fill(particle->GetPdgCode());
 441     }
 442   } // end of track loop
 443
 444
 445   delete iter;
 446   iter = 0;
 447
 448   if (fSecondaries)
 449     FillSecondaries();
 450
 451   if (fSigmaVertex)
 452     FillSigmaVertex();
 453
 454
 455   delete listOfTracks;
 456   listOfTracks = 0;
 457
 458   return kTRUE;
 459 }
 460
 461 void AlidNdEtaSystematicsSelector::FillSecondaries()
 462 {
 463   // fills the secondary histograms
 464
 465   AliStack* stack = GetStack();
 466
 467   Int_t particlesPrimaries = 0;
 468   Int_t particlesSecondaries = 0;
 469   Int_t particlesPrimariesPtCut = 0;
 470   Int_t particlesSecondariesPtCut = 0;
 471
 472   for (Int_t iParticle = 0; iParticle < stack->GetNtrack(); iParticle++)
 473   {
 474     TParticle* particle = stack->Particle(iParticle);
 475     if (!particle)
 476     {
 477       AliDebug(AliLog::kError, Form("UNEXPECTED: particle with label %d not found in stack (particle loop).", iParticle));
 478       continue;
 479     }
 480
 481     if (TMath::Abs(particle->Eta()) > 0.9)
 482       continue;
 483
 484     Bool_t isPrimary = kFALSE;
 485     if (iParticle < stack->GetNprimary())
 486     {
 487       if (AliPWG0Helper::IsPrimaryCharged(particle, stack->GetNprimary()) == kFALSE)
 488         continue;
 489
 490       isPrimary = kTRUE;
 491     }
 492
 493     if (isPrimary)
 494       particlesPrimaries++;
 495     else
 496       particlesSecondaries++;
 497
 498     if (particle->Pt() < 0.3)
 499       continue;
 500
 501     if (isPrimary)
 502       particlesPrimariesPtCut++;
 503     else
 504       particlesSecondariesPtCut++;
 505   }
 506
 507   fSecondaries->Fill(1, particlesPrimaries);
 508   fSecondaries->Fill(2, particlesSecondaries);
 509   fSecondaries->Fill(3, particlesPrimariesPtCut);
 510   fSecondaries->Fill(4, particlesSecondariesPtCut);
 511
 512   Int_t allTracksPrimaries = 0;
 513   Int_t allTracksSecondaries = 0;
 514   Int_t acceptedTracksPrimaries = 0;
 515   Int_t acceptedTracksSecondaries = 0;
 516
 517   for (Int_t iTrack = 0; iTrack < fESD->GetNumberOfTracks(); iTrack++)
 518   {
 519     AliESDtrack* esdTrack = fESD->GetTrack(iTrack);
 520     if (!esdTrack)
 521     {
 522       AliDebug(AliLog::kError, Form("UNEXPECTED: Could not get track %d.", iTrack));
 523       continue;
 524     }
 525
 526     Int_t label = TMath::Abs(esdTrack->GetLabel());
 527     TParticle* particle = stack->Particle(label);
 528     if (!particle)
 529     {
 530       AliDebug(AliLog::kError, Form("UNEXPECTED: particle with label %d not found in stack (track loop).", label));
 531       continue;
 532     }
 533
 534    if (label < stack->GetNprimary())
 535       allTracksPrimaries++;
 536     else
 537       allTracksSecondaries++;
 538
 539     if (!fEsdTrackCuts->AcceptTrack(esdTrack))
 540       continue;
 541
 542     if (label < stack->GetNprimary())
 543       acceptedTracksPrimaries++;
 544     else
 545       acceptedTracksSecondaries++;
 546   }
 547
 548   fSecondaries->Fill(5, allTracksPrimaries);
 549   fSecondaries->Fill(6, allTracksSecondaries);
 550   fSecondaries->Fill(7, acceptedTracksPrimaries);
 551   fSecondaries->Fill(8, acceptedTracksSecondaries);
 552
 553   //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);
 554 }
 555
 556 void AlidNdEtaSystematicsSelector::FillSigmaVertex()
 557 {
 558   // fills the fSigmaVertex histogram
 559
 560   // save the old value
 561   Float_t oldSigmaVertex = fEsdTrackCuts->GetMinNsigmaToVertex();
 562
 563   // set to maximum
 564   fEsdTrackCuts->SetMinNsigmaToVertex(5);
 565
 566   TObjArray* list = fEsdTrackCuts->GetAcceptedTracks(fESD);
 567
 568   TIterator* iter = list->MakeIterator();
 569   TObject* obj = 0;
 570   while ((obj = iter->Next()) != 0)
 571   {
 572     AliESDtrack* esdTrack = dynamic_cast<AliESDtrack*> (obj);
 573     if (!esdTrack)
 574       continue;
 575
 576     Float_t sigma2Vertex = fEsdTrackCuts->GetSigmaToVertex(esdTrack);
 577
 578     for (Double_t nSigma = 0.1; nSigma < 5.05; nSigma += 0.1)
 579     {
 580       if (sigma2Vertex < nSigma)
 581         fSigmaVertex->Fill(nSigma);
 582     }
 583   }
 584
 585   delete iter;
 586   iter = 0;
 587
 588   delete list;
 589   list = 0;
 590
 591   // set back the old value
 592   fEsdTrackCuts->SetMinNsigmaToVertex(oldSigmaVertex);
 593 }
 594
 595 void AlidNdEtaSystematicsSelector::SlaveTerminate()
 596 {
 597   // The SlaveTerminate() function is called after all entries or objects
 598   // have been processed. When running with PROOF SlaveTerminate() is called
 599   // on each slave server.
 600
 601   AliSelector::SlaveTerminate();
 602
 603   // Add the histograms to the output on each slave server
 604   if (!fOutput)
 605   {
 606     AliDebug(AliLog::kError, Form("ERROR: Output list not initialized."));
 607     return;
 608   }
 609
 610   if (fSecondaries)
 611     fOutput->Add(fSecondaries);
 612
 613   for (Int_t i=0; i<4; ++i)
 614     if (fdNdEtaCorrectionSpecies[i])
 615       fOutput->Add(fdNdEtaCorrectionSpecies[i]);
 616
 617   if (fSigmaVertex)
 618     fOutput->Add(fSigmaVertex);
 619
 620   for (Int_t i=0; i<3; ++i) {
 621     if (fdNdEtaCorrectionVertexReco[i])
 622       fOutput->Add(fdNdEtaCorrectionVertexReco[i]);
 623
 624     if (fdNdEtaCorrectionTriggerBias[i])
 625       fOutput->Add(fdNdEtaCorrectionTriggerBias[i]);
 626   }
 627
 628 }
 629
 630 void AlidNdEtaSystematicsSelector::Terminate()
 631 {
 632   // The Terminate() function is the last function to be called during
 633   // a query. It always runs on the client, it can be used to present
 634   // the results graphically or save the results to file.
 635
 636   AliSelector::Terminate();
 637
 638   fSecondaries = dynamic_cast<TH2F*> (fOutput->FindObject("fSecondaries"));
 639   for (Int_t i=0; i<4; ++i)
 640     fdNdEtaCorrectionSpecies[i] = dynamic_cast<AlidNdEtaCorrection*> (fOutput->FindObject(Form("correction_%d", i)));
 641   fSigmaVertex = dynamic_cast<TH1F*> (fOutput->FindObject("fSigmaVertex"));
 642
 643   fdNdEtaCorrectionVertexReco[0] = dynamic_cast<AlidNdEtaCorrection*> (fOutput->FindObject("vertexRecoND"));
 644   fdNdEtaCorrectionVertexReco[1] = dynamic_cast<AlidNdEtaCorrection*> (fOutput->FindObject("vertexRecoSD"));
 645   fdNdEtaCorrectionVertexReco[2] = dynamic_cast<AlidNdEtaCorrection*> (fOutput->FindObject("vertexRecoDD"));
 646
 647   fdNdEtaCorrectionTriggerBias[0] = dynamic_cast<AlidNdEtaCorrection*> (fOutput->FindObject("triggerBiasND"));
 648   fdNdEtaCorrectionTriggerBias[1] = dynamic_cast<AlidNdEtaCorrection*> (fOutput->FindObject("triggerBiasSD"));
 649   fdNdEtaCorrectionTriggerBias[2] = dynamic_cast<AlidNdEtaCorrection*> (fOutput->FindObject("triggerBiasDD"));
 650
 651   if (fPIDParticles)
 652   {
 653     TDatabasePDG* pdgDB = new TDatabasePDG;
 654
 655     for (Int_t i=0; i <= fPIDParticles->GetNbinsX()+1; ++i)
 656       if (fPIDParticles->GetBinContent(i) > 0)
 657         printf("PDG = %d (%s): generated: %d, reconstructed: %d, ratio: %f\n",
 658                (Int_t) fPIDParticles->GetBinCenter(i),
 659                pdgDB->GetParticle((Int_t) fPIDParticles->GetBinCenter(i))->GetName(),
 660                (Int_t) fPIDParticles->GetBinContent(i),
 661                (Int_t) fPIDTracks->GetBinContent(i),
 662                ((fPIDTracks->GetBinContent(i) > 0) ? fPIDParticles->GetBinContent(i) / fPIDTracks->GetBinContent(i) : -1));
 663
 664     delete pdgDB;
 665     pdgDB = 0;
 666   }
 667
 668   for (Int_t i=0; i<3; ++i) {
 669     if (fdNdEtaCorrectionVertexReco[i])
 670       fdNdEtaCorrectionVertexReco[i]->Finish();
 671
 672     if (fdNdEtaCorrectionTriggerBias[i])
 673       fdNdEtaCorrectionTriggerBias[i]->Finish();
 674   }
 675
 676   TFile* fout = TFile::Open("systematics.root", "RECREATE");
 677
 678   if (fEsdTrackCuts)
 679     fEsdTrackCuts->SaveHistograms("esd_track_cuts");
 680
 681   if (fSecondaries)
 682     fSecondaries->Write();
 683
 684   if (fSigmaVertex)
 685     fSigmaVertex->Write();
 686
 687   for (Int_t i=0; i<4; ++i)
 688     if (fdNdEtaCorrectionSpecies[i])
 689       fdNdEtaCorrectionSpecies[i]->SaveHistograms();
 690
 691   for (Int_t i=0; i<3; ++i) {
 692     if (fdNdEtaCorrectionVertexReco[i])
 693       fdNdEtaCorrectionVertexReco[i]->SaveHistograms();
 694
 695     if (fdNdEtaCorrectionTriggerBias[i])
 696       fdNdEtaCorrectionTriggerBias[i]->SaveHistograms();
 697   }
 698
 699   fout->Write();
 700   fout->Close();
 701
 702   if (fSecondaries)
 703   {
 704     new TCanvas;
 705     fSecondaries->Draw("COLZ");
 706   }
 707
 708   if (fSigmaVertex)
 709   {
 710     new TCanvas;
 711     fSigmaVertex->Draw();
 712   }
 713 }
 714
 715 Bool_t AlidNdEtaSystematicsSelector::SignOK(TParticlePDG* particle)
 716 {
 717   // returns if a particle with this sign should be counted
 718   // this is determined by the value of fSignMode, which should have the same sign
 719   // as the charge
 720   // if fSignMode is 0 all particles are counted
 721
 722   if (fSignMode == 0)
 723     return kTRUE;
 724
 725   if (!particle)
 726   {
 727     printf("WARNING: not counting a particle that does not have a pdg particle\n");
 728     return kFALSE;
 729   }
 730
 731   Double_t charge = particle->Charge();
 732
 733   if (fSignMode > 0)
 734     if (charge < 0)
 735       return kFALSE;
 736
 737   if (fSignMode < 0)
 738     if (charge > 0)
 739       return kFALSE;
 740
 741   return kTRUE;
 742 }