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