PWG0/dNdEta/AlidNdEtaCorrectionTask.cxx

   1 /* $Id$ */
   2
   3 #include "AlidNdEtaCorrectionTask.h"
   4
   5 #include <TCanvas.h>
   6 #include <TChain.h>
   7 #include <TFile.h>
   8 #include <TH1F.h>
   9 #include <TH2F.h>
  10 #include <TH3F.h>
  11 #include <TProfile.h>
  12 #include <TParticle.h>
  13 #include <TParticlePDG.h>
  14 #include <TDatabasePDG.h>
  15 #include <TRandom.h>
  16
  17 #include <AliLog.h>
  18 #include <AliESDVertex.h>
  19 #include <AliESDEvent.h>
  20 #include <AliESDRun.h>
  21 #include <AliStack.h>
  22 #include <AliHeader.h>
  23 #include <AliGenEventHeader.h>
  24 #include <AliMultiplicity.h>
  25 #include <AliAnalysisManager.h>
  26 #include <AliMCEventHandler.h>
  27 #include <AliMCEvent.h>
  28 #include <AliESDInputHandler.h>
  29
  30 #include "AliESDtrackCuts.h"
  31 #include "AliPWG0Helper.h"
  32 #include "dNdEta/dNdEtaAnalysis.h"
  33 #include "dNdEta/AlidNdEtaCorrection.h"
  34 #include "AliTriggerAnalysis.h"
  35 #include "AliPhysicsSelection.h"
  36
  37 ClassImp(AlidNdEtaCorrectionTask)
  38
  39 AlidNdEtaCorrectionTask::AlidNdEtaCorrectionTask() :
  40   AliAnalysisTask(),
  41   fESD(0),
  42   fOutput(0),
  43   fOption(),
  44   fAnalysisMode((AliPWG0Helper::AnalysisMode) (AliPWG0Helper::kTPC | AliPWG0Helper::kFieldOn)),
  45   fTrigger(AliTriggerAnalysis::kMB1),
  46   fFillPhi(kFALSE),
  47   fDeltaPhiCut(-1),
  48   fSymmetrize(kFALSE),
  49   fMultAxisEta1(kFALSE),
  50   fDiffTreatment(AliPWG0Helper::kMCFlags),
  51   fSignMode(0),
  52   fOnlyPrimaries(kFALSE),
  53   fStatError(0),
  54   fSystSkipParticles(kFALSE),
  55   fEsdTrackCuts(0),
  56   fdNdEtaCorrection(0),
  57   fdNdEtaAnalysisMC(0),
  58   fdNdEtaAnalysisESD(0),
  59   fPIDParticles(0),
  60   fPIDTracks(0),
  61   fVertexCorrelation(0),
  62   fVertexCorrelationShift(0),
  63   fVertexProfile(0),
  64   fVertexShift(0),
  65   fVertexShiftNorm(0),
  66   fEtaCorrelation(0),
  67   fEtaCorrelationShift(0),
  68   fEtaProfile(0),
  69   fEtaResolution(0),
  70   fDeltaPhiCorrelation(0),
  71   fpTResolution(0),
  72   fEsdTrackCutsPrim(0),
  73   fEsdTrackCutsSec(0),
  74   fTemp1(0),
  75   fTemp2(0),
  76   fMultAll(0),
  77   fMultTr(0),
  78   fMultVtx(0),
  79   fEventStats(0),
  80   fEsdTrackCutsCheck(0),
  81   fEtaCorrelationAllESD(0),
  82   fpTCorrelation(0),
  83   fpTCorrelationShift(0),
  84   fpTCorrelationAllESD(0),
  85   fpTCorrelationShiftAllESD(0),
  86   fPtMin(0.15),
  87   fPtMC(0),
  88   fEtaMC(0),
  89   fPtESD(0),
  90   fEtaESD(0),
  91   fVtxMC(0),
  92   fNumberEventMC(0),
  93   fNumberEvent(0),
  94   fEventNumber(-1),
  95   fWeightSecondaries(kFALSE)
  96 {
  97   //
  98   // Constructor. Initialization of pointers
  99   //
 100
 101   for (Int_t i=0; i<4; i++)
 102     fdNdEtaCorrectionSpecial[i] = 0;
 103
 104   for (Int_t i=0; i<8; i++)
 105     fDeltaPhi[i] = 0;
 106
 107   AliLog::SetClassDebugLevel("AlidNdEtaCorrectionTask", AliLog::kWarning);
 108 }
 109
 110 AlidNdEtaCorrectionTask::AlidNdEtaCorrectionTask(const char* opt) :
 111   AliAnalysisTask("AlidNdEtaCorrectionTask", ""),
 112   fESD(0),
 113   fOutput(0),
 114   fOption(opt),
 115   fAnalysisMode((AliPWG0Helper::AnalysisMode) (AliPWG0Helper::kTPC | AliPWG0Helper::kFieldOn)),
 116   fTrigger(AliTriggerAnalysis::kMB1),
 117   fFillPhi(kFALSE),
 118   fDeltaPhiCut(0),
 119   fSymmetrize(kFALSE),
 120   fMultAxisEta1(kFALSE),
 121   fDiffTreatment(AliPWG0Helper::kMCFlags),
 122   fSignMode(0),
 123   fOnlyPrimaries(kFALSE),
 124   fStatError(0),
 125   fSystSkipParticles(kFALSE),
 126   fEsdTrackCuts(0),
 127   fdNdEtaCorrection(0),
 128   fdNdEtaAnalysisMC(0),
 129   fdNdEtaAnalysisESD(0),
 130   fPIDParticles(0),
 131   fPIDTracks(0),
 132   fVertexCorrelation(0),
 133   fVertexCorrelationShift(0),
 134   fVertexProfile(0),
 135   fVertexShift(0),
 136   fVertexShiftNorm(0),
 137   fEtaCorrelation(0),
 138   fEtaCorrelationShift(0),
 139   fEtaProfile(0),
 140   fEtaResolution(0),
 141   fDeltaPhiCorrelation(0),
 142   fpTResolution(0),
 143   fEsdTrackCutsPrim(0),
 144   fEsdTrackCutsSec(0),
 145   fTemp1(0),
 146   fTemp2(0),
 147   fMultAll(0),
 148   fMultTr(0),
 149   fMultVtx(0),
 150   fEventStats(0),
 151   fEsdTrackCutsCheck(0),
 152   fEtaCorrelationAllESD(0),
 153   fpTCorrelation(0),
 154   fpTCorrelationShift(0),
 155   fpTCorrelationAllESD(0),
 156   fpTCorrelationShiftAllESD(0),
 157   fPtMin(0.15),
 158   fPtMC(0),
 159   fEtaMC(0),
 160   fPtESD(0),
 161   fEtaESD(0),
 162   fVtxMC(0),
 163   fNumberEventMC(0),
 164   fNumberEvent(0),
 165   fEventNumber(-1),
 166   fWeightSecondaries(kFALSE)
 167 {
 168   //
 169   // Constructor. Initialization of pointers
 170   //
 171
 172   // Define input and output slots here
 173   DefineInput(0, TChain::Class());
 174   DefineOutput(0, TList::Class());
 175
 176   for (Int_t i=0; i<4; i++)
 177     fdNdEtaCorrectionSpecial[i] = 0;
 178
 179   for (Int_t i=0; i<8; i++)
 180     fDeltaPhi[i] = 0;
 181
 182   AliLog::SetClassDebugLevel("AlidNdEtaCorrectionTask", AliLog::kWarning);
 183 }
 184
 185 AlidNdEtaCorrectionTask::~AlidNdEtaCorrectionTask()
 186 {
 187   //
 188   // Destructor
 189   //
 190
 191   // histograms are in the output list and deleted when the output
 192   // list is deleted by the TSelector dtor
 193
 194   if (fOutput && !AliAnalysisManager::GetAnalysisManager()->IsProofMode()) {
 195     delete fOutput;
 196     fOutput = 0;
 197   }
 198
 199 }
 200
 201 //________________________________________________________________________
 202 void AlidNdEtaCorrectionTask::ConnectInputData(Option_t *)
 203 {
 204   // Connect ESD
 205   // Called once
 206
 207   Printf("AlidNdEtaCorrectionTask::ConnectInputData called");
 208
 209   AliESDInputHandler *esdH = dynamic_cast<AliESDInputHandler*> (AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler());
 210
 211   if (!esdH) {
 212     Printf("ERROR: Could not get ESDInputHandler");
 213   } else {
 214     fESD = esdH->GetEvent();
 215
 216     // Enable only the needed branches
 217     esdH->SetActiveBranches("AliESDHeader Vertex");
 218
 219     if (fAnalysisMode & AliPWG0Helper::kSPD)
 220       esdH->SetActiveBranches("AliESDHeader Vertex AliMultiplicity");
 221
 222     if (fAnalysisMode & AliPWG0Helper::kTPC || fAnalysisMode & AliPWG0Helper::kTPCITS) {
 223       esdH->SetActiveBranches("AliESDHeader Vertex Tracks");
 224     }
 225   }
 226
 227   // disable info messages of AliMCEvent (per event)
 228   AliLog::SetClassDebugLevel("AliMCEvent", AliLog::kWarning - AliLog::kDebug + 1);
 229 }
 230
 231 void AlidNdEtaCorrectionTask::CreateOutputObjects()
 232 {
 233   // create result objects and add to output list
 234
 235   AliDebug(2,Form("*********************************** fOption = %s", fOption.Data()));
 236   if (fOption.Contains("only-positive"))
 237   {
 238     Printf("INFO: Processing only positive particles.");
 239     fSignMode = 1;
 240   }
 241   else if (fOption.Contains("only-negative"))
 242   {
 243     Printf("INFO: Processing only negative particles.");
 244     fSignMode = -1;
 245   }
 246
 247   if (fOption.Contains("stat-error-1"))
 248   {
 249     Printf("INFO: Evaluation statistical errors. Mode: 1.");
 250     fStatError = 1;
 251   }
 252   else if (fOption.Contains("stat-error-2"))
 253   {
 254     Printf("INFO: Evaluation statistical errors. Mode: 2.");
 255     fStatError = 2;
 256   }
 257
 258   fOutput = new TList;
 259   fOutput->SetOwner();
 260
 261   fdNdEtaCorrection = new AlidNdEtaCorrection("dndeta_correction", "dndeta_correction", fAnalysisMode);
 262   fOutput->Add(fdNdEtaCorrection);
 263
 264   fPIDParticles = new TH1F("fPIDParticles", "PID of generated primary particles", 10001, -5000.5, 5000.5);
 265   fOutput->Add(fPIDParticles);
 266
 267   fPIDTracks = new TH1F("fPIDTracks", "MC PID of reconstructed tracks", 10001, -5000.5, 5000.5);
 268   fOutput->Add(fPIDTracks);
 269
 270   fdNdEtaAnalysisMC = new dNdEtaAnalysis("dndetaMC", "dndetaMC", fAnalysisMode);
 271   fOutput->Add(fdNdEtaAnalysisMC);
 272
 273   fdNdEtaAnalysisESD = new dNdEtaAnalysis("dndetaESD", "dndetaESD", fAnalysisMode);
 274   fOutput->Add(fdNdEtaAnalysisESD);
 275
 276   if (fEsdTrackCuts)
 277   {
 278     fEsdTrackCutsPrim = static_cast<AliESDtrackCuts*> (fEsdTrackCuts->Clone("fEsdTrackCutsPrim"));
 279     fEsdTrackCutsSec  = static_cast<AliESDtrackCuts*> (fEsdTrackCuts->Clone("fEsdTrackCutsSec"));
 280     fEsdTrackCutsCheck  = static_cast<AliESDtrackCuts*> (fEsdTrackCuts->Clone("fEsdTrackCutsCheck"));
 281     fEsdTrackCutsCheck->SetPtRange(0.15);
 282     fEsdTrackCutsCheck->SetEtaRange(-1.,1.);
 283     fOutput->Add(fEsdTrackCutsPrim);
 284     fOutput->Add(fEsdTrackCutsSec);
 285   }
 286
 287   if (fOption.Contains("process-types")) {
 288     fdNdEtaCorrectionSpecial[0] = new AlidNdEtaCorrection("dndeta_correction_ND", "dndeta_correction_ND", fAnalysisMode);
 289     fdNdEtaCorrectionSpecial[1] = new AlidNdEtaCorrection("dndeta_correction_SD", "dndeta_correction_SD", fAnalysisMode);
 290     fdNdEtaCorrectionSpecial[2] = new AlidNdEtaCorrection("dndeta_correction_DD", "dndeta_correction_DD", fAnalysisMode);
 291
 292     fOutput->Add(fdNdEtaCorrectionSpecial[0]);
 293     fOutput->Add(fdNdEtaCorrectionSpecial[1]);
 294     fOutput->Add(fdNdEtaCorrectionSpecial[2]);
 295   }
 296
 297   if (fOption.Contains("particle-species")) {
 298     fdNdEtaCorrectionSpecial[0] = new AlidNdEtaCorrection("dndeta_correction_pi", "dndeta_correction_pi", fAnalysisMode);
 299     fdNdEtaCorrectionSpecial[1] = new AlidNdEtaCorrection("dndeta_correction_K", "dndeta_correction_K", fAnalysisMode);
 300     fdNdEtaCorrectionSpecial[2] = new AlidNdEtaCorrection("dndeta_correction_p", "dndeta_correction_p", fAnalysisMode);
 301     fdNdEtaCorrectionSpecial[3] = new AlidNdEtaCorrection("dndeta_correction_other", "dndeta_correction_other", fAnalysisMode);
 302
 303     for (Int_t i=0; i<4; i++)
 304       fOutput->Add(fdNdEtaCorrectionSpecial[i]);
 305   }
 306
 307
 308   //fTemp1 = new TH2F("fTemp1", "fTemp1", 4, 0.5, 4.5, 101, -1.5, 99.5); // nsd study
 309   fTemp1 = new TH2F("fTemp1", "fTemp1", 300, -15, 15, 80, -2.0, 2.0);
 310   fOutput->Add(fTemp1);
 311
 312   fTemp2 = new TH2F("fTemp2", "fTemp2", 300, -15, 15, 80, -2.0, 2.0);
 313   fOutput->Add(fTemp2);
 314
 315   fVertexCorrelation = new TH2F("fVertexCorrelation", "fVertexCorrelation;MC z-vtx;ESD z-vtx", 120, -30, 30, 120, -30, 30);
 316   fOutput->Add(fVertexCorrelation);
 317   fVertexCorrelationShift = new TH3F("fVertexCorrelationShift", "fVertexCorrelationShift;MC z-vtx;MC z-vtx - ESD z-vtx;rec. tracks", 120, -30, 30, 100, -1, 1, 100, -0.5, 99.5);
 318   fOutput->Add(fVertexCorrelationShift);
 319   fVertexProfile = new TProfile("fVertexProfile", "fVertexProfile;MC z-vtx;MC z-vtx - ESD z-vtx", 40, -20, 20);
 320   fOutput->Add(fVertexProfile);
 321   fVertexShift = new TH1F("fVertexShift", "fVertexShift;(MC z-vtx - ESD z-vtx);Entries", 201, -2, 2);
 322   fOutput->Add(fVertexShift);
 323   fVertexShiftNorm = new TH2F("fVertexShiftNorm", "fVertexShiftNorm;(MC z-vtx - ESD z-vtx);rec. tracks;Entries", 200, -100, 100, 100, -0.5, 99.5);
 324   fOutput->Add(fVertexShiftNorm);
 325
 326   fEtaCorrelation = new TH2F("fEtaCorrelation", "fEtaCorrelation;MC #eta;ESD #eta", 120, -3, 3, 120, -3, 3);
 327   fOutput->Add(fEtaCorrelation);
 328   fEtaCorrelationAllESD = new TH2F("fEtaCorrelationAllESD", "fEtaCorrelationAllESD;MC #eta;ESD #eta", 120, -3, 3, 120, -3, 3);
 329   fOutput->Add(fEtaCorrelationAllESD);
 330   fEtaCorrelationShift = new TH2F("fEtaCorrelationShift", "fEtaCorrelationShift;MC #eta;MC #eta - ESD #eta", 120, -3, 3, 100, -0.1, 0.1);
 331   fOutput->Add(fEtaCorrelationShift);
 332   fEtaProfile = new TProfile("fEtaProfile", "fEtaProfile;MC #eta;MC #eta - ESD #eta", 120, -3, 3);
 333   fOutput->Add(fEtaProfile);
 334   fEtaResolution = new TH1F("fEtaResolution", "fEtaResolution;MC #eta - ESD #eta", 201, -0.2, 0.2);
 335   fOutput->Add(fEtaResolution);
 336
 337   fpTResolution = new TH2F("fpTResolution", ";MC p_{T} (GeV/c);(MC p_{T} - ESD p_{T}) / MC p_{T}", 160, 0, 20, 201, -0.2, 0.2);
 338   fOutput->Add(fpTResolution);
 339
 340   fpTCorrelation = new TH2F("fpTCorrelation", "fpTCorrelation;MC p_{T} (GeV/c);ESD p_{T}", 160, 0, 20, 160, 0, 20);
 341   fOutput->Add(fpTCorrelation);
 342   fpTCorrelationShift = new TH2F("fpTCorrelationShift", "fpTCorrelationShift;MC p_{T} (GeV/c);MC p_{T} - ESD p_{T}", 160, 0, 20, 100, -1, 1);
 343   fOutput->Add(fpTCorrelationShift);
 344   fpTCorrelationAllESD = new TH2F("fpTCorrelationAllESD", "fpTCorrelationAllESD;MC p_{T} (GeV/c);ESD p_{T}", 160, 0, 20, 160, 0, 20);
 345   fOutput->Add(fpTCorrelationAllESD);
 346   fpTCorrelationShiftAllESD = new TH2F("fpTCorrelationShiftAllESD", "fpTCorrelationShiftAllESD;MC p_{T} (GeV/c);MC p_{T} - ESD p_{T}", 160, 0, 20, 100, -1, 1);
 347   fOutput->Add(fpTCorrelationShiftAllESD);
 348
 349   fMultAll = new TH1F("fMultAll", "fMultAll", 500, -0.5, 499.5);
 350   fOutput->Add(fMultAll);
 351   fMultVtx = new TH1F("fMultVtx", "fMultVtx", 500, -0.5, 499.5);
 352   fOutput->Add(fMultVtx);
 353   fMultTr = new TH1F("fMultTr", "fMultTr", 500, -0.5, 499.5);
 354   fOutput->Add(fMultTr);
 355
 356   for (Int_t i=0; i<8; i++)
 357   {
 358     fDeltaPhi[i] = new TH2F(Form("fDeltaPhi_%d", i), ";#Delta phi;#phi;Entries", 2000, -0.1, 0.1, 18 * 5, 0, TMath::TwoPi());
 359     fOutput->Add(fDeltaPhi[i]);
 360   }
 361
 362   fEventStats = new TH2F("fEventStats", "fEventStats;event type;status;count", 109, -6.5, 102.5, 4, -0.5, 3.5);
 363   fOutput->Add(fEventStats);
 364   fEventStats->GetXaxis()->SetBinLabel(1, "INEL"); // x = -5
 365   fEventStats->GetXaxis()->SetBinLabel(2, "NSD");  // x = -4
 366   fEventStats->GetXaxis()->SetBinLabel(3, "ND");   // x = -3
 367   fEventStats->GetXaxis()->SetBinLabel(4, "SD");   // x = -2
 368   fEventStats->GetXaxis()->SetBinLabel(5, "DD");   // x = -1
 369
 370   fEventStats->GetXaxis()->SetBinLabel(108, "INEL=0");   // x = -101
 371   fEventStats->GetXaxis()->SetBinLabel(109, "INEL>0");   // x = -102
 372
 373   for (Int_t i=-1; i<100; i++)
 374     fEventStats->GetXaxis()->SetBinLabel(7+i, Form("%d", i));
 375
 376   fEventStats->GetYaxis()->SetBinLabel(1, "nothing");
 377   fEventStats->GetYaxis()->SetBinLabel(2, "trg");
 378   fEventStats->GetYaxis()->SetBinLabel(3, "vtx");
 379   fEventStats->GetYaxis()->SetBinLabel(4, "trgvtx");
 380
 381   if (fEsdTrackCuts)
 382   {
 383     fEsdTrackCuts->SetName("fEsdTrackCuts");
 384     // TODO like this we send an empty object back...
 385     fOutput->Add(fEsdTrackCuts->Clone());
 386   }
 387   fPtMC = new TH1F("fPtMC", "Pt from MC for selected tracks;MC p_{T} (GeV/c)", 160, 0, 20);
 388   fOutput->Add(fPtMC);
 389   fEtaMC = new TH1F("fEtaMC", "Eta from MC for selected tracks;MC #eta", 120, -3, 3);
 390   fOutput->Add(fEtaMC);
 391   fPtESD = new TH1F("fPtESD", "Pt from ESD for selected tracks;ESD p_{T} (GeV/c)", 160, 0, 20);
 392   fOutput->Add(fPtESD);
 393   fEtaESD = new TH1F("fEtaESD", "Eta from ESD for selected tracks;ESD #eta", 120, -3, 3);
 394   fOutput->Add(fEtaESD);
 395
 396   fVtxMC = new TH1F("fVtxMC", "Vtx,z from MC for all events;MC vtx_z (cm)", 100, -30, 30);
 397   fOutput->Add(fVtxMC);
 398
 399   fNumberEventMC = new TH1F("fNumberEventMC","Number of event accepted at MC level",600000,-0.5,600000-0.5);
 400   fOutput->Add(fNumberEventMC);
 401
 402   fNumberEvent = new TH1F("fNumberEvent","Number of event accepted at Reco level",600000,-0.5,600000-0.5);
 403   fOutput->Add(fNumberEvent);
 404 }
 405
 406 void AlidNdEtaCorrectionTask::Exec(Option_t*)
 407 {
 408   // process the event
 409
 410   fEventNumber++;
 411   // Check prerequisites
 412   if (!fESD)
 413   {
 414     AliDebug(AliLog::kError, "ESD branch not available");
 415     return;
 416   }
 417
 418   if (fOnlyPrimaries)
 419     Printf("WARNING: Processing only primaries. For systematical studies only!");
 420
 421   if (fStatError > 0)
 422     Printf("WARNING: Statistical error evaluation active!");
 423
 424   AliInputEventHandler* inputHandler = dynamic_cast<AliInputEventHandler*> (AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler());
 425   if (!inputHandler)
 426   {
 427     Printf("ERROR: Could not receive input handler");
 428     return;
 429   }
 430
 431   Bool_t eventTriggered = inputHandler->IsEventSelected();
 432
 433   static AliTriggerAnalysis* triggerAnalysis = 0;
 434   if (!triggerAnalysis)
 435   {
 436     AliPhysicsSelection* physicsSelection = dynamic_cast<AliPhysicsSelection*> (inputHandler->GetEventSelection());
 437     if (physicsSelection)
 438       triggerAnalysis = physicsSelection->GetTriggerAnalysis();
 439   }
 440
 441   if (eventTriggered)
 442     eventTriggered = triggerAnalysis->IsTriggerFired(fESD, fTrigger);
 443
 444   if (!eventTriggered)
 445     Printf("No trigger");
 446
 447   // post the data already here
 448   PostData(0, fOutput);
 449
 450   // MC info
 451   AliMCEventHandler* eventHandler = dynamic_cast<AliMCEventHandler*> (AliAnalysisManager::GetAnalysisManager()->GetMCtruthEventHandler());
 452   if (!eventHandler) {
 453     Printf("ERROR: Could not retrieve MC event handler");
 454     return;
 455   }
 456
 457   AliMCEvent* mcEvent = eventHandler->MCEvent();
 458   if (!mcEvent) {
 459     Printf("ERROR: Could not retrieve MC event");
 460     return;
 461   }
 462
 463   AliStack* stack = mcEvent->Stack();
 464   if (!stack)
 465   {
 466     AliDebug(AliLog::kError, "Stack not available");
 467     return;
 468   }
 469
 470   AliHeader* header = mcEvent->Header();
 471   if (!header)
 472   {
 473     AliDebug(AliLog::kError, "Header not available");
 474     return;
 475   }
 476
 477   // get process type
 478   Int_t processType = AliPWG0Helper::GetEventProcessType(fESD, header, stack, fDiffTreatment);
 479
 480   AliDebug(AliLog::kDebug+1, Form("Found process type %d", processType));
 481
 482   if (processType == AliPWG0Helper::kInvalidProcess)
 483   {
 484     AliDebug(AliLog::kWarning, "Unknown process type. Setting to ND");
 485     processType = AliPWG0Helper::kND;
 486   }
 487
 488   // get the MC vertex
 489   AliGenEventHeader* genHeader = header->GenEventHeader();
 490   TArrayF vtxMC(3);
 491   genHeader->PrimaryVertex(vtxMC);
 492   fVtxMC->Fill(vtxMC[2]);
 493   AliDebug(2,Form("MC vtx: x = %.6f, y = %.6f, z = %.6f",vtxMC[0],vtxMC[1],vtxMC[2]));
 494
 495   // get the ESD vertex
 496   Bool_t eventVertex = kFALSE;
 497   Double_t vtx[3];
 498   const AliESDVertex* vtxESD = AliPWG0Helper::GetVertex(fESD, fAnalysisMode);
 499   if (vtxESD && AliPWG0Helper::TestVertex(vtxESD, fAnalysisMode))
 500   {
 501     vtxESD->GetXYZ(vtx);
 502     eventVertex = kTRUE;
 503
 504     // remove vertices outside +- 15 cm
 505     if (TMath::Abs(vtx[2]) > 15)
 506     {
 507       eventVertex = kFALSE;
 508       vtxESD = 0;
 509     }
 510   }
 511   else
 512     vtxESD = 0;
 513
 514   // create list of (label, eta, pt) tuples
 515   Int_t inputCount = 0;
 516   Int_t* labelArr = 0;
 517   Int_t* labelArr2 = 0; // only for case of SPD
 518   Float_t* etaArr = 0;
 519   Float_t* thirdDimArr = 0;
 520   Float_t* deltaPhiArr = 0;
 521   if (fAnalysisMode & AliPWG0Helper::kSPD)
 522   {
 523     if (vtxESD)
 524     {
 525       // get tracklets
 526       const AliMultiplicity* mult = fESD->GetMultiplicity();
 527       if (!mult)
 528       {
 529         AliDebug(AliLog::kError, "AliMultiplicity not available");
 530         return;
 531       }
 532
 533       labelArr = new Int_t[mult->GetNumberOfTracklets()];
 534       labelArr2 = new Int_t[mult->GetNumberOfTracklets()];
 535       etaArr = new Float_t[mult->GetNumberOfTracklets()];
 536       thirdDimArr = new Float_t[mult->GetNumberOfTracklets()];
 537       deltaPhiArr = new Float_t[mult->GetNumberOfTracklets()];
 538
 539       Bool_t foundInEta10 = kFALSE;
 540
 541       // get multiplicity from SPD tracklets
 542       for (Int_t i=0; i<mult->GetNumberOfTracklets(); ++i)
 543       {
 544         //printf("%d %f %f %f\n", i, mult->GetTheta(i), mult->GetPhi(i), mult->GetDeltaPhi(i));
 545
 546         Float_t phi = mult->GetPhi(i);
 547         if (phi < 0)
 548           phi += TMath::Pi() * 2;
 549         Float_t deltaPhi = mult->GetDeltaPhi(i);
 550
 551         if (TMath::Abs(deltaPhi) > 1)
 552           printf("WARNING: Very high Delta Phi: %d %f %f %f\n", i, mult->GetTheta(i), mult->GetPhi(i), deltaPhi);
 553
 554         if (fOnlyPrimaries)
 555           if (mult->GetLabel(i, 0) < 0 || mult->GetLabel(i, 0) != mult->GetLabel(i, 1) || !stack->IsPhysicalPrimary(mult->GetLabel(i, 0)))
 556             continue;
 557
 558         if (fDeltaPhiCut > 0 && (TMath::Abs(deltaPhi) > fDeltaPhiCut || TMath::Abs(mult->GetDeltaTheta(i)) > fDeltaPhiCut / 0.08 * 0.025))
 559           continue;
 560
 561         if (fSystSkipParticles && gRandom->Uniform() < 0.0153)
 562         {
 563           Printf("Skipped tracklet!");
 564           continue;
 565         }
 566
 567         // TEST exclude potentially inefficient phi region
 568         //if (phi > 5.70 || phi < 0.06)
 569         //  continue;
 570
 571         // we have to repeat the trigger here, because the tracklet might have been kicked out fSystSkipParticles
 572         if (TMath::Abs(mult->GetEta(i)) < 1)
 573           foundInEta10 = kTRUE;
 574
 575         etaArr[inputCount] = mult->GetEta(i);
 576         if (fSymmetrize)
 577           etaArr[inputCount] = TMath::Abs(etaArr[inputCount]);
 578         labelArr[inputCount] = mult->GetLabel(i, 0);
 579         labelArr2[inputCount] = mult->GetLabel(i, 1);
 580         thirdDimArr[inputCount] = phi;
 581         deltaPhiArr[inputCount] = deltaPhi;
 582         ++inputCount;
 583       }
 584
 585       /*
 586       for (Int_t i=0; i<mult->GetNumberOfSingleClusters(); ++i)
 587       {
 588         if (TMath::Abs(TMath::Log(TMath::Tan(mult->GetThetaSingle(i)/2.))) < 1);
 589         {
 590           foundInEta10 = kTRUE;
 591           break;
 592         }
 593       }
 594       */
 595
 596       if (fSystSkipParticles && (fTrigger & AliTriggerAnalysis::kOneParticle) && !foundInEta10)
 597         eventTriggered = kFALSE;
 598     }
 599   }
 600   else if (fAnalysisMode & AliPWG0Helper::kTPC || fAnalysisMode & AliPWG0Helper::kTPCITS)
 601   {
 602     if (!fEsdTrackCuts)
 603     {
 604       AliDebug(AliLog::kError, "fESDTrackCuts not available");
 605       return;
 606     }
 607
 608     Bool_t foundInEta10 = kFALSE;
 609
 610     if (vtxESD)
 611     {
 612       // control histograms on pT
 613       Int_t nfromstack = stack->GetNtrack();
 614       AliDebug(3,Form(" from stack we have %d tracks\n",nfromstack));
 615       for (Int_t itrack = 0; itrack < fESD->GetNumberOfTracks(); itrack++){
 616         AliESDtrack* esdTrackcheck = dynamic_cast<AliESDtrack*> (fESD->GetTrack(itrack));
 617         if (!esdTrackcheck){
 618           AliDebug(AliLog::kError, Form("ERROR: Could not retrieve track %d.", itrack));
 619           continue;
 620         }
 621         if (fOnlyPrimaries){
 622           Int_t label = TMath::Abs(esdTrackcheck->GetLabel());
 623           AliDebug(4,Form("label = %d\n",label));
 624           if (label == 0 || label > nfromstack) continue;
 625           if (stack->IsPhysicalPrimary(label) == kFALSE) continue;
 626         }
 627
 628         Int_t label = TMath::Abs(esdTrackcheck->GetLabel());
 629         if (label == 0 || label > nfromstack) continue;
 630         if (!stack->Particle(label)){
 631           AliDebug(4,Form("WARNING: No particle for %d", label));
 632         }
 633         else{
 634           if (!fEsdTrackCuts->AcceptTrack(esdTrackcheck)){
 635             TParticle* particle = stack->Particle(label);
 636             if (fEsdTrackCutsCheck->AcceptTrack(esdTrackcheck)){
 637               //if (TMath::Abs(particle->Eta() < 0.8) && particle->Pt() > fPtMin){
 638               Float_t ptMC = particle->Pt();
 639               Float_t etaMC = particle->Eta();
 640               Float_t ptESD = esdTrackcheck->Pt();
 641               Float_t etaESD = esdTrackcheck->Eta();
 642               fEtaCorrelationAllESD->Fill(etaMC,etaESD);
 643               fpTCorrelationAllESD->Fill(ptMC,ptESD);
 644               fpTCorrelationShiftAllESD->Fill(ptMC,ptMC-ptESD);
 645             }
 646           }
 647         }
 648       } // end loop over all ESDs
 649
 650       // get multiplicity from ESD tracks
 651       TObjArray* list = fEsdTrackCuts->GetAcceptedTracks(fESD, (fAnalysisMode & AliPWG0Helper::kTPC));
 652       Int_t nGoodTracks = list->GetEntries();
 653
 654       Printf("Accepted %d tracks", nGoodTracks);
 655
 656       labelArr = new Int_t[nGoodTracks];
 657       labelArr2 = new Int_t[nGoodTracks];
 658       etaArr = new Float_t[nGoodTracks];
 659       thirdDimArr = new Float_t[nGoodTracks];
 660       deltaPhiArr = new Float_t[nGoodTracks];
 661
 662       // loop over esd tracks
 663       for (Int_t i=0; i<nGoodTracks; i++)
 664       {
 665         AliESDtrack* esdTrack = dynamic_cast<AliESDtrack*> (list->At(i));
 666         if (!esdTrack)
 667         {
 668           AliDebug(AliLog::kError, Form("ERROR: Could not retrieve track %d.", i));
 669           continue;
 670         }
 671
 672         AliDebug(3,Form("particle %d: pt = %.6f, eta = %.6f",i,esdTrack->Pt(), esdTrack->Eta()));
 673         // 2 Options for INEL>0 trigger - choose one
 674         // 1. HL
 675         //if (esdTrack->Pt() < 0.15)
 676         // foundInEta10 = kTRUE;
 677         // 2. MB Working Group definition
 678         if (esdTrack->Pt() < fPtMin)
 679           continue;
 680
 681         if (fOnlyPrimaries)
 682         {
 683           Int_t label = TMath::Abs(esdTrack->GetLabel());
 684           if (label == 0)
 685             continue;
 686
 687           if (stack->IsPhysicalPrimary(label) == kFALSE)
 688             continue;
 689         }
 690
 691         Int_t label = TMath::Abs(esdTrack->GetLabel());
 692         if (!stack->Particle(label)){
 693           AliDebug(3,Form("WARNING: No particle for %d", label));
 694         }
 695         else{
 696           TParticle* particle = stack->Particle(label);
 697           Float_t ptMC = particle->Pt();
 698           Float_t etaMC = particle->Eta();
 699           fPtMC->Fill(ptMC);
 700           fEtaMC->Fill(etaMC);
 701           fPtESD->Fill(esdTrack->Pt());
 702           fEtaESD->Fill(esdTrack->Eta());
 703         }
 704
 705         // 2 Options for INEL>0 trigger - choose one
 706         // 1. HL
 707         //if (TMath::Abs(esdTrack->Eta()) < 1 && esdTrack->Pt() > 0.15)
 708         // foundInEta10 = kTRUE;
 709         // 2. MB Working Group definition
 710         if (TMath::Abs(esdTrack->Eta()) < 0.8 && esdTrack->Pt() > fPtMin)
 711           foundInEta10 = kTRUE;
 712
 713         etaArr[inputCount] = esdTrack->Eta();
 714         if (fSymmetrize)
 715           etaArr[inputCount] = TMath::Abs(etaArr[inputCount]);
 716         labelArr[inputCount] = TMath::Abs(esdTrack->GetLabel());
 717         labelArr2[inputCount] = labelArr[inputCount]; // no second label for tracks
 718         thirdDimArr[inputCount] = esdTrack->Pt();
 719         deltaPhiArr[inputCount] = 0; // no delta phi for tracks
 720         ++inputCount;
 721       }
 722
 723       delete list;
 724
 725       // TODO this code crashes for TPCITS because particles are requested from the stack for some labels that are out of bound
 726       if (0 && eventTriggered)
 727       {
 728         // collect values for primaries and secondaries
 729         for (Int_t iTrack = 0; iTrack < fESD->GetNumberOfTracks(); iTrack++)
 730         {
 731           AliESDtrack* track = 0;
 732
 733           if (fAnalysisMode & AliPWG0Helper::kTPC)
 734             track = AliESDtrackCuts::GetTPCOnlyTrack(fESD, iTrack);
 735           else if (fAnalysisMode & AliPWG0Helper::kTPCITS)
 736             track = fESD->GetTrack(iTrack);
 737
 738           if (!track)
 739             continue;
 740
 741           Int_t label = TMath::Abs(track->GetLabel());
 742           if (!stack->Particle(label) || !stack->Particle(label)->GetPDG())
 743           {
 744             Printf("WARNING: No particle for %d", label);
 745             if (stack->Particle(label))
 746               stack->Particle(label)->Print();
 747             continue;
 748           }
 749
 750           if (stack->Particle(label)->GetPDG()->Charge() == 0)
 751             continue;
 752
 753           if (TMath::Abs(track->Eta()) < 0.8 && track->Pt() > 0.15)
 754           {
 755             if (stack->IsPhysicalPrimary(label))
 756             {
 757               // primary
 758               if (fEsdTrackCutsPrim->AcceptTrack(track))
 759               {
 760 //                 if (AliESDtrackCuts::GetSigmaToVertex(track) > 900)
 761 //                 {
 762 //                   Printf("Track %d has nsigma of %f. Printing track and vertex...", iTrack, AliESDtrackCuts::GetSigmaToVertex(track));
 763 //                   Float_t b[2];
 764 //                   Float_t r[3];
 765 //                   track->GetImpactParameters(b, r);
 766 //                   Printf("Impact parameter %f %f and resolution: %f %f %f", b[0], b[1], r[0], r[1], r[2]);
 767 //                   track->Print("");
 768 //                   if (vtxESD)
 769 //                     vtxESD->Print();
 770 //                 }
 771               }
 772             }
 773             else
 774             {
 775               // secondary
 776               fEsdTrackCutsSec->AcceptTrack(track);
 777             }
 778           }
 779
 780           // TODO mem leak in the continue statements above
 781           if (fAnalysisMode & AliPWG0Helper::kTPC)
 782             delete track;
 783         }
 784       }
 785     }
 786
 787     if (!foundInEta10)
 788       eventTriggered = kFALSE;
 789     else{
 790      //Printf("The event triggered. Its number in file is %d",fESD->GetEventNumberInFile());
 791       fNumberEvent->Fill(fESD->GetEventNumberInFile());
 792     }
 793   }
 794   else
 795     return;
 796
 797   // fill process type
 798   Int_t biny = (Int_t) eventTriggered + 2 * (Int_t) eventVertex;
 799   // INEL
 800   fEventStats->Fill(-6, biny);
 801   // NSD
 802   if (processType != AliPWG0Helper::kSD)
 803     fEventStats->Fill(-5, biny);
 804   // SD, ND, DD
 805   if (processType == AliPWG0Helper::kND)
 806     fEventStats->Fill(-4, biny);
 807   if (processType == AliPWG0Helper::kSD)
 808     fEventStats->Fill(-3, biny);
 809   if (processType == AliPWG0Helper::kDD)
 810     fEventStats->Fill(-2, biny);
 811
 812   // loop over mc particles
 813   Int_t nPrim  = stack->GetNprimary();
 814   Int_t nAccepted = 0;
 815
 816   Bool_t oneParticleEvent = kFALSE;
 817   for (Int_t iMc = 0; iMc < nPrim; ++iMc)
 818   {
 819     //Printf("Getting particle %d", iMc);
 820     TParticle* particle = stack->Particle(iMc);
 821
 822     if (!particle)
 823       continue;
 824
 825     if (AliPWG0Helper::IsPrimaryCharged(particle, nPrim) == kFALSE)
 826       continue;
 827
 828     // for INEL > 0, MB Working Group definition use the second option
 829     // 1. standard
 830     //if (TMath::Abs(particle->Eta()) < 1.0)
 831     // 2. MB Working Group definition
 832     if (TMath::Abs(particle->Eta()) < 0.8 && particle->Pt() > fPtMin)
 833     {
 834       oneParticleEvent = kTRUE;
 835       fNumberEventMC->Fill(fESD->GetEventNumberInFile());
 836       break;
 837     }
 838   }
 839
 840   if (TMath::Abs(vtxMC[2]) < 5.5)
 841   {
 842     if (oneParticleEvent)
 843       fEventStats->Fill(102, biny);
 844     else
 845       fEventStats->Fill(101, biny);
 846   }
 847
 848   for (Int_t iMc = 0; iMc < nPrim; ++iMc)
 849   {
 850     //Printf("Getting particle %d", iMc);
 851     TParticle* particle = stack->Particle(iMc);
 852
 853     if (!particle)
 854       continue;
 855
 856     if (AliPWG0Helper::IsPrimaryCharged(particle, nPrim) == kFALSE)
 857     {
 858       //if (TMath::Abs(particle->GetPdgCode()) > 3000 && TMath::Abs(particle->Eta()) < 1.0)
 859       //  fPIDParticles->Fill(particle->GetPdgCode());
 860       continue;
 861     }
 862
 863     if (SignOK(particle->GetPDG()) == kFALSE)
 864       continue;
 865
 866     // for INEL > 0, MB Working Group definition use the second option
 867     // 1. standard
 868     //if (fPIDParticles && TMath::Abs(particle->Eta()) < 1.0)
 869     // 2. MB Working Group definition
 870     if (fPIDParticles && TMath::Abs(particle->Eta()) < 0.8 && particle->Pt() > fPtMin)
 871       fPIDParticles->Fill(particle->GetPdgCode());
 872
 873     Float_t eta = particle->Eta();
 874     if (fSymmetrize)
 875       eta = TMath::Abs(eta);
 876
 877     Float_t thirdDim = -1;
 878     if (fAnalysisMode & AliPWG0Helper::kSPD)
 879     {
 880       if (fFillPhi)
 881       {
 882         thirdDim = particle->Phi();
 883       }
 884       else
 885         thirdDim = inputCount;
 886     }
 887     else
 888       thirdDim = particle->Pt();
 889
 890     // calculate y
 891     //Float_t y = 0.5 * TMath::Log((particle->Energy() + particle->Pz()) / (particle->Energy() - particle->Pz()));
 892     //fTemp1->Fill(eta);
 893     //fTemp2->Fill(y);
 894
 895     Int_t processType2 = processType;
 896     if (oneParticleEvent)
 897       processType2 |= AliPWG0Helper::kOnePart;
 898
 899     fdNdEtaCorrection->FillMCParticle(vtxMC[2], eta, thirdDim, eventTriggered, eventVertex, processType2);
 900
 901     if (fOption.Contains("process-types"))
 902     {
 903       // non diffractive
 904       if (processType==AliPWG0Helper::kND)
 905         fdNdEtaCorrectionSpecial[0]->FillMCParticle(vtxMC[2], eta, thirdDim, eventTriggered, eventVertex, processType2);
 906
 907       // single diffractive
 908       if (processType==AliPWG0Helper::kSD)
 909         fdNdEtaCorrectionSpecial[1]->FillMCParticle(vtxMC[2], eta, thirdDim, eventTriggered, eventVertex, processType2);
 910
 911       // double diffractive
 912       if (processType==AliPWG0Helper::kDD)
 913         fdNdEtaCorrectionSpecial[2]->FillMCParticle(vtxMC[2], eta, thirdDim, eventTriggered, eventVertex, processType2);
 914     }
 915
 916     if (fOption.Contains("particle-species"))
 917     {
 918       Int_t id = -1;
 919       switch (TMath::Abs(particle->GetPdgCode()))
 920       {
 921         case 211:   id = 0; break;
 922         case 321:   id = 1; break;
 923         case 2212:  id = 2; break;
 924         default:    id = 3; break;
 925       }
 926       fdNdEtaCorrectionSpecial[id]->FillMCParticle(vtxMC[2], eta, thirdDim, eventTriggered, eventVertex, processType2);
 927     }
 928
 929     if (eventTriggered)
 930       if (eventVertex)
 931         fdNdEtaAnalysisMC->FillTrack(vtxMC[2], eta, thirdDim);
 932
 933     // TODO this value might be needed lower for the SPD study (only used for control histograms anyway)
 934     if (TMath::Abs(eta) < 1.5) // && pt > 0.2)
 935       nAccepted++;
 936   }
 937
 938   if (AliPWG0Helper::GetLastProcessType() >= -1)
 939     fEventStats->Fill(AliPWG0Helper::GetLastProcessType(), biny);
 940
 941   fMultAll->Fill(nAccepted);
 942   if (eventTriggered) {
 943     fMultTr->Fill(nAccepted);
 944     if (eventVertex)
 945       fMultVtx->Fill(nAccepted);
 946   }
 947
 948   Int_t processed = 0;
 949
 950   Bool_t* primCount = 0;
 951   if (fStatError > 0)
 952   {
 953     primCount = new Bool_t[nPrim];
 954     for (Int_t i=0; i<nPrim; ++i)
 955       primCount[i] = kFALSE;
 956   }
 957
 958   Int_t nEta05 = 0;
 959   Int_t nEta10 = 0;
 960   for (Int_t i=0; i<inputCount; ++i)
 961   {
 962     if (TMath::Abs(etaArr[i]) < 0.5)
 963       nEta05++;
 964     if (TMath::Abs(etaArr[i]) < 1.0)
 965       nEta10++;
 966   }
 967
 968   for (Int_t i=0; i<inputCount; ++i)
 969   {
 970     Int_t label = labelArr[i];
 971     Int_t label2 = labelArr2[i];
 972
 973     if (label < 0)
 974     {
 975       Printf("WARNING: cannot find corresponding mc particle for track(let) %d with label %d.", i, label);
 976       continue;
 977     }
 978
 979     TParticle* particle = stack->Particle(label);
 980     if (!particle)
 981     {
 982       AliDebug(AliLog::kError, Form("UNEXPECTED: particle with label %d not found in stack (track loop).", label));
 983       continue;
 984     }
 985
 986     // for INEL > 0, MB Working Group definition use the second option
 987     // 1. standard
 988     //if (TMath::Abs(particle->Eta()) < 1.0)
 989     // 2. INEL > 0 MB Working Group definition
 990     if (TMath::Abs(particle->Eta()) < 0.8 && particle->Pt()>fPtMin)
 991       fPIDTracks->Fill(particle->GetPdgCode());
 992
 993     // find particle that is filled in the correction map
 994     // this should be the particle that has been reconstructed
 995     // for TPC this is clear and always identified by the track's label
 996     // for SPD the labels might not be identical. In this case the particle closest to the beam line that is a primary is taken:
 997     //  L1 L2 (P = primary, S = secondary)
 998     //   P P' : --> P
 999     //   P S  : --> P
1000     //   S P  : --> P
1001     //   S S' : --> S
1002
1003     if (label != label2 && label2 >= 0)
1004     {
1005       if (stack->IsPhysicalPrimary(label) == kFALSE && stack->IsPhysicalPrimary(label2))
1006       {
1007         particle = stack->Particle(label2);
1008         if (!particle)
1009         {
1010           AliDebug(AliLog::kError, Form("UNEXPECTED: particle with label %d not found in stack (track loop).", label2));
1011           continue;
1012         }
1013       }
1014     }
1015
1016     if (eventTriggered && eventVertex)
1017     {
1018       if (SignOK(particle->GetPDG()) == kFALSE)
1019         continue;
1020
1021       processed++;
1022
1023       // resolutions
1024       if (fSymmetrize)
1025         fEtaResolution->Fill(TMath::Abs(particle->Eta()) - etaArr[i]);
1026       else
1027         fEtaResolution->Fill(particle->Eta() - etaArr[i]);
1028
1029       if (fAnalysisMode & AliPWG0Helper::kTPC || fAnalysisMode & AliPWG0Helper::kTPCITS){
1030         // for INEL > 0, MB Working Group definition use the second option
1031         // 1. standard
1032         //if (TMath::Abs(particle->Eta() < 0.9) && particle->Pt() > 0)
1033         // 2. INEL > 0 MB WOrking Group definition
1034         if (TMath::Abs(particle->Eta() < 0.8) && particle->Pt() > 0){
1035           fpTResolution->Fill(particle->Pt(), (particle->Pt() - thirdDimArr[i]) / particle->Pt());
1036           fpTCorrelation->Fill(particle->Pt(),thirdDimArr[i]);
1037           fpTCorrelationShift->Fill(particle->Pt(),particle->Pt()-thirdDimArr[i]);
1038         }
1039       }
1040
1041       Float_t eta = -999;
1042       Float_t thirdDim = -1;
1043
1044       Bool_t firstIsPrim = stack->IsPhysicalPrimary(label);
1045       // in case of same label the MC values are filled, otherwise (background) the reconstructed values
1046       if (label == label2)
1047       {
1048         eta = particle->Eta();
1049         if (fSymmetrize)
1050           eta = TMath::Abs(eta);
1051
1052         if (fAnalysisMode & AliPWG0Helper::kSPD)
1053         {
1054           if (fFillPhi)
1055           {
1056             thirdDim = particle->Phi();
1057           }
1058           else
1059             thirdDim = inputCount;
1060         }
1061         else
1062           thirdDim = particle->Pt();
1063       }
1064       else
1065       {
1066         if (fAnalysisMode & AliPWG0Helper::kSPD && !fFillPhi)
1067         {
1068           thirdDim = (fMultAxisEta1) ? nEta10 : inputCount;
1069         }
1070         else
1071           thirdDim = thirdDimArr[i];
1072
1073         eta = etaArr[i];
1074       }
1075
1076       Bool_t fillTrack = kTRUE;
1077
1078       // statistical error evaluation active?
1079       if (fStatError > 0)
1080       {
1081         Bool_t statErrorDecision = kFALSE;
1082
1083         // primary and not yet counted
1084         if (label == label2 && firstIsPrim && !primCount[label])
1085         {
1086           statErrorDecision = kTRUE;
1087           primCount[label] = kTRUE;
1088         }
1089
1090         // in case of 1 we count only unique primaries, in case of 2 all the rest
1091         if (fStatError == 1)
1092         {
1093           fillTrack = statErrorDecision;
1094         }
1095         else if (fStatError == 2)
1096           fillTrack = !statErrorDecision;
1097       }
1098
1099       if (fillTrack)
1100       {
1101         Double_t weight = 1.;
1102         if (fWeightSecondaries){
1103           if (!firstIsPrim){
1104             weight = GetSecondaryCorrection(thirdDim);
1105           }
1106         }
1107         fdNdEtaCorrection->FillTrackedParticle(vtxMC[2], eta, thirdDim, weight);
1108         fTemp2->Fill(vtxMC[2], eta);
1109       }
1110
1111       // eta comparison for tracklets with the same label (others are background)
1112       if (label == label2)
1113       {
1114         Float_t eta2 = particle->Eta();
1115         if (fSymmetrize)
1116           eta2 = TMath::Abs(eta2);
1117
1118         fEtaProfile->Fill(eta2, eta2 - etaArr[i]);
1119         fEtaCorrelation->Fill(eta2, etaArr[i]);
1120         fEtaCorrelationShift->Fill(eta2, eta2 - etaArr[i]);
1121       }
1122
1123       if (fSymmetrize)
1124         fdNdEtaAnalysisESD->FillTrack(vtxMC[2], TMath::Abs(particle->Eta()), thirdDim);
1125       else
1126         fdNdEtaAnalysisESD->FillTrack(vtxMC[2], particle->Eta(), thirdDim);
1127
1128       if (fOption.Contains("process-types"))
1129       {
1130         // non diffractive
1131         if (processType == AliPWG0Helper::kND)
1132           fdNdEtaCorrectionSpecial[0]->FillTrackedParticle(vtxMC[2], eta, thirdDim);
1133
1134         // single diffractive
1135         if (processType == AliPWG0Helper::kSD)
1136           fdNdEtaCorrectionSpecial[1]->FillTrackedParticle(vtxMC[2], eta, thirdDim);
1137
1138         // double diffractive
1139         if (processType == AliPWG0Helper::kDD)
1140           fdNdEtaCorrectionSpecial[2]->FillTrackedParticle(vtxMC[2], eta, thirdDim);
1141       }
1142
1143       if (fOption.Contains("particle-species"))
1144       {
1145         // find mother first
1146         TParticle* mother = AliPWG0Helper::FindPrimaryMother(stack, label);
1147
1148         Int_t id = -1;
1149         switch (TMath::Abs(mother->GetPdgCode()))
1150         {
1151           case 211:   id = 0; break;
1152           case 321:   id = 1; break;
1153           case 2212:  id = 2; break;
1154           default:    id = 3; break;
1155         }
1156         fdNdEtaCorrectionSpecial[id]->FillTrackedParticle(vtxMC[2], eta, thirdDim);
1157       }
1158
1159       // control histograms
1160       Int_t hist = -1;
1161       if (label == label2)
1162       {
1163         if (firstIsPrim)
1164         {
1165           hist = 0;
1166         }
1167         else
1168           hist = 1;
1169       }
1170       else if (label2 >= 0)
1171       {
1172         Bool_t secondIsPrim = stack->IsPhysicalPrimary(label2);
1173         if (firstIsPrim && secondIsPrim)
1174         {
1175           hist = 2;
1176         }
1177         else if (firstIsPrim && !secondIsPrim)
1178         {
1179           hist = 3;
1180
1181           // check if secondary is caused by the primary or it is a fake combination
1182           //Printf("PS case --> Label 1 is %d, label 2 is %d", label, label2);
1183           Int_t mother = label2;
1184           while (stack->Particle(mother) && stack->Particle(mother)->GetMother(0) >= 0)
1185           {
1186             //Printf("  %d created by %d, %d", mother, stack->Particle(mother)->GetMother(0), stack->Particle(mother)->GetMother(1));
1187             if (stack->Particle(mother)->GetMother(0) == label)
1188             {
1189               /*Printf("The untraceable decay was:");
1190               Printf("   from:");
1191               particle->Print();
1192               Printf("   to (status code %d):", stack->Particle(mother)->GetStatusCode());
1193               stack->Particle(mother)->Print();*/
1194               hist = 4;
1195             }
1196             mother = stack->Particle(mother)->GetMother(0);
1197           }
1198         }
1199         else if (!firstIsPrim && secondIsPrim)
1200         {
1201           hist = 5;
1202         }
1203         else if (!firstIsPrim && !secondIsPrim)
1204         {
1205           hist = 6;
1206         }
1207
1208       }
1209       else
1210         hist = 7;
1211
1212       fDeltaPhi[hist]->Fill(deltaPhiArr[i], thirdDimArr[i]);
1213     }
1214   }
1215
1216   if (primCount)
1217   {
1218     delete[] primCount;
1219     primCount = 0;
1220   }
1221
1222   if (processed < inputCount)
1223     Printf("Only %d out of %d track(let)s were used", processed, inputCount);
1224
1225   if (eventTriggered && eventVertex)
1226   {
1227     Double_t diff = vtxMC[2] - vtx[2];
1228     fVertexShift->Fill(diff);
1229
1230     fVertexCorrelation->Fill(vtxMC[2], vtx[2]);
1231     fVertexCorrelationShift->Fill(vtxMC[2], vtxMC[2] - vtx[2], inputCount);
1232     fVertexProfile->Fill(vtxMC[2], vtxMC[2] - vtx[2]);
1233
1234     if (vtxESD->IsFromVertexerZ() && inputCount > 0)
1235       fVertexShiftNorm->Fill(diff, vtxESD->GetNContributors());
1236   }
1237
1238   Int_t multAxis = inputCount;
1239   if (fMultAxisEta1)
1240     multAxis = nEta10;
1241
1242   if (eventTriggered && eventVertex)
1243   {
1244     fdNdEtaAnalysisMC->FillEvent(vtxMC[2], multAxis);
1245     fdNdEtaAnalysisESD->FillEvent(vtxMC[2], multAxis);
1246   }
1247
1248   Int_t processType2 = processType;
1249   if (oneParticleEvent)
1250     processType2 |= AliPWG0Helper::kOnePart;
1251
1252   // stuff regarding the vertex reco correction and trigger bias correction
1253   fdNdEtaCorrection->FillEvent(vtxMC[2], multAxis, eventTriggered, eventVertex, processType2);
1254
1255   if (fOption.Contains("process-types"))
1256   {
1257     // non diffractive
1258     if (processType == AliPWG0Helper::kND)
1259       fdNdEtaCorrectionSpecial[0]->FillEvent(vtxMC[2], multAxis, eventTriggered, eventVertex, processType2);
1260
1261     // single diffractive
1262     if (processType == AliPWG0Helper::kSD)
1263       fdNdEtaCorrectionSpecial[1]->FillEvent(vtxMC[2], multAxis, eventTriggered, eventVertex, processType2);
1264
1265     // double diffractive
1266     if (processType == AliPWG0Helper::kDD)
1267       fdNdEtaCorrectionSpecial[2]->FillEvent(vtxMC[2], multAxis, eventTriggered, eventVertex, processType2);
1268   }
1269
1270   if (fOption.Contains("particle-species"))
1271     for (Int_t id=0; id<4; id++)
1272       fdNdEtaCorrectionSpecial[id]->FillEvent(vtxMC[2], multAxis, eventTriggered, eventVertex, processType2);
1273
1274   if (etaArr)
1275     delete[] etaArr;
1276   if (labelArr)
1277     delete[] labelArr;
1278   if (labelArr2)
1279     delete[] labelArr2;
1280   if (thirdDimArr)
1281     delete[] thirdDimArr;
1282   if (deltaPhiArr)
1283     delete[] deltaPhiArr;
1284 }
1285
1286 void AlidNdEtaCorrectionTask::Terminate(Option_t *)
1287 {
1288   // The Terminate() function is the last function to be called during
1289   // a query. It always runs on the client, it can be used to present
1290   // the results graphically or save the results to file.
1291
1292   fOutput = dynamic_cast<TList*> (GetOutputData(0));
1293   if (!fOutput) {
1294     Printf("ERROR: fOutput not available");
1295     return;
1296   }
1297
1298   fdNdEtaCorrection = dynamic_cast<AlidNdEtaCorrection*> (fOutput->FindObject("dndeta_correction"));
1299   fdNdEtaAnalysisMC = dynamic_cast<dNdEtaAnalysis*> (fOutput->FindObject("dndetaMC"));
1300   fdNdEtaAnalysisESD = dynamic_cast<dNdEtaAnalysis*> (fOutput->FindObject("dndetaESD"));
1301   if (!fdNdEtaCorrection || !fdNdEtaAnalysisMC || !fdNdEtaAnalysisESD)
1302   {
1303     AliDebug(AliLog::kError, "Could not read object from output list");
1304     return;
1305   }
1306
1307   fdNdEtaCorrection->Finish();
1308
1309   TString fileName;
1310   fileName.Form("correction_map%s.root", fOption.Data());
1311   TFile* fout = new TFile(fileName, "RECREATE");
1312
1313   fEsdTrackCuts = dynamic_cast<AliESDtrackCuts*> (fOutput->FindObject("fEsdTrackCuts"));
1314   if (fEsdTrackCuts)
1315     fEsdTrackCuts->SaveHistograms("esd_track_cuts");
1316
1317   fEsdTrackCutsPrim = dynamic_cast<AliESDtrackCuts*> (fOutput->FindObject("fEsdTrackCutsPrim"));
1318   if (fEsdTrackCutsPrim)
1319     fEsdTrackCutsPrim->SaveHistograms("esd_track_cuts_primaries");
1320
1321   fEsdTrackCutsSec = dynamic_cast<AliESDtrackCuts*> (fOutput->FindObject("fEsdTrackCutsSec"));
1322   if (fEsdTrackCutsSec)
1323     fEsdTrackCutsSec->SaveHistograms("esd_track_cuts_secondaries");
1324
1325   fdNdEtaCorrection->SaveHistograms();
1326   fdNdEtaAnalysisMC->SaveHistograms();
1327   fdNdEtaAnalysisESD->SaveHistograms();
1328
1329   if (fOutput->FindObject("dndeta_correction_ND"))
1330   {
1331     fdNdEtaCorrectionSpecial[0] = dynamic_cast<AlidNdEtaCorrection*> (fOutput->FindObject("dndeta_correction_ND"));
1332     fdNdEtaCorrectionSpecial[1] = dynamic_cast<AlidNdEtaCorrection*> (fOutput->FindObject("dndeta_correction_SD"));
1333     fdNdEtaCorrectionSpecial[2] = dynamic_cast<AlidNdEtaCorrection*> (fOutput->FindObject("dndeta_correction_DD"));
1334   }
1335   else
1336   {
1337     fdNdEtaCorrectionSpecial[0] = dynamic_cast<AlidNdEtaCorrection*> (fOutput->FindObject("dndeta_correction_pi"));
1338     fdNdEtaCorrectionSpecial[1] = dynamic_cast<AlidNdEtaCorrection*> (fOutput->FindObject("dndeta_correction_K"));
1339     fdNdEtaCorrectionSpecial[2] = dynamic_cast<AlidNdEtaCorrection*> (fOutput->FindObject("dndeta_correction_p"));
1340     fdNdEtaCorrectionSpecial[3] = dynamic_cast<AlidNdEtaCorrection*> (fOutput->FindObject("dndeta_correction_other"));
1341   }
1342   for (Int_t i=0; i<4; ++i)
1343     if (fdNdEtaCorrectionSpecial[i])
1344       fdNdEtaCorrectionSpecial[i]->SaveHistograms();
1345
1346   fTemp1 = dynamic_cast<TH1*> (fOutput->FindObject("fTemp1"));
1347   if (fTemp1)
1348     fTemp1->Write();
1349   fTemp2 = dynamic_cast<TH1*> (fOutput->FindObject("fTemp2"));
1350   if (fTemp2)
1351     fTemp2->Write();
1352
1353   fVertexCorrelation = dynamic_cast<TH2F*> (fOutput->FindObject("fVertexCorrelation"));
1354   if (fVertexCorrelation)
1355     fVertexCorrelation->Write();
1356   fVertexCorrelationShift = dynamic_cast<TH3F*> (fOutput->FindObject("fVertexCorrelationShift"));
1357   if (fVertexCorrelationShift)
1358   {
1359     ((TH2*) fVertexCorrelationShift->Project3D("yx"))->FitSlicesY();
1360     fVertexCorrelationShift->Write();
1361   }
1362   fVertexProfile = dynamic_cast<TProfile*> (fOutput->FindObject("fVertexProfile"));
1363   if (fVertexProfile)
1364     fVertexProfile->Write();
1365   fVertexShift = dynamic_cast<TH1F*> (fOutput->FindObject("fVertexShift"));
1366   if (fVertexShift)
1367     fVertexShift->Write();
1368   fVertexShiftNorm = dynamic_cast<TH2F*> (fOutput->FindObject("fVertexShiftNorm"));
1369   if (fVertexShiftNorm)
1370   {
1371     fVertexShiftNorm->ProjectionX();
1372     fVertexShiftNorm->Write();
1373   }
1374
1375   fEtaCorrelation = dynamic_cast<TH2F*> (fOutput->FindObject("fEtaCorrelation"));
1376   if (fEtaCorrelation)
1377     fEtaCorrelation->Write();
1378   fEtaCorrelationAllESD = dynamic_cast<TH2F*> (fOutput->FindObject("fEtaCorrelationAllESD"));
1379   if (fEtaCorrelationAllESD)
1380     fEtaCorrelationAllESD->Write();
1381   fEtaCorrelationShift = dynamic_cast<TH2F*> (fOutput->FindObject("fEtaCorrelationShift"));
1382   if (fEtaCorrelationShift)
1383   {
1384     fEtaCorrelationShift->FitSlicesY();
1385     fEtaCorrelationShift->Write();
1386   }
1387   fEtaProfile = dynamic_cast<TProfile*> (fOutput->FindObject("fEtaProfile"));
1388   if (fEtaProfile)
1389     fEtaProfile->Write();
1390   fEtaResolution = dynamic_cast<TH1F*> (fOutput->FindObject("fEtaResolution"));
1391   if (fEtaResolution)
1392     fEtaResolution->Write();
1393   fpTCorrelation = dynamic_cast<TH2F*> (fOutput->FindObject("fpTCorrelation"));
1394   if (fpTCorrelation)
1395     fpTCorrelation->Write();
1396   fpTCorrelationShift = dynamic_cast<TH2F*> (fOutput->FindObject("fpTCorrelationShift"));
1397   if (fpTCorrelationShift)
1398   {
1399     fpTCorrelationShift->FitSlicesY();
1400     fpTCorrelationShift->Write();
1401   }
1402   fpTResolution = dynamic_cast<TH2F*> (fOutput->FindObject("fpTResolution"));
1403   if (fpTResolution)
1404   {
1405     fpTResolution->FitSlicesY();
1406     fpTResolution->Write();
1407   }
1408
1409   fpTCorrelationAllESD = dynamic_cast<TH2F*> (fOutput->FindObject("fpTCorrelationAllESD"));
1410   if (fpTCorrelationAllESD)
1411     fpTCorrelationAllESD->Write();
1412   fpTCorrelationShiftAllESD = dynamic_cast<TH2F*> (fOutput->FindObject("fpTCorrelationShiftAllESD"));
1413   if (fpTCorrelationShiftAllESD)
1414   {
1415     fpTCorrelationShiftAllESD->FitSlicesY();
1416     fpTCorrelationShiftAllESD->Write();
1417   }
1418   fMultAll = dynamic_cast<TH1F*> (fOutput->FindObject("fMultAll"));
1419   if (fMultAll)
1420     fMultAll->Write();
1421
1422   fMultTr = dynamic_cast<TH1F*> (fOutput->FindObject("fMultTr"));
1423   if (fMultTr)
1424     fMultTr->Write();
1425
1426   fMultVtx = dynamic_cast<TH1F*> (fOutput->FindObject("fMultVtx"));
1427   if (fMultVtx)
1428     fMultVtx->Write();
1429
1430   for (Int_t i=0; i<8; ++i)
1431   {
1432     fDeltaPhi[i] = dynamic_cast<TH2*> (fOutput->FindObject(Form("fDeltaPhi_%d", i)));
1433     if (fDeltaPhi[i])
1434       fDeltaPhi[i]->Write();
1435   }
1436
1437   fEventStats = dynamic_cast<TH2F*> (fOutput->FindObject("fEventStats"));
1438   if (fEventStats)
1439     fEventStats->Write();
1440
1441   fPIDParticles = dynamic_cast<TH1F*> (fOutput->FindObject("fPIDParticles"));
1442   if (fPIDParticles)
1443     fPIDParticles->Write();
1444
1445   fPIDTracks = dynamic_cast<TH1F*> (fOutput->FindObject("fPIDTracks"));
1446   if (fPIDTracks)
1447     fPIDTracks->Write();
1448
1449   fPtMC = dynamic_cast<TH1F*> (fOutput->FindObject("fPtMC"));
1450   if (fPtMC)
1451     fPtMC->Write();
1452   fEtaMC = dynamic_cast<TH1F*> (fOutput->FindObject("fEtaMC"));
1453   if (fEtaMC)
1454     fEtaMC->Write();
1455   fPtESD = dynamic_cast<TH1F*> (fOutput->FindObject("fPtESD"));
1456   if (fPtESD)
1457     fPtESD->Write();
1458   fEtaESD = dynamic_cast<TH1F*> (fOutput->FindObject("fEtaESD"));
1459   if (fEtaESD)
1460     fEtaESD->Write();
1461   fVtxMC = dynamic_cast<TH1F*> (fOutput->FindObject("fVtxMC"));
1462   if (fVtxMC)
1463     fVtxMC->Write();
1464
1465   fNumberEventMC = dynamic_cast<TH1F*> (fOutput->FindObject("fNumberEventMC"));
1466   if (fNumberEventMC)
1467     fNumberEventMC->Write();
1468   fNumberEvent = dynamic_cast<TH1F*> (fOutput->FindObject("fNumberEvent"));
1469   if (fNumberEvent)
1470     fNumberEvent->Write();
1471
1472  //fdNdEtaCorrection->DrawHistograms();
1473
1474   Printf("Writing result to %s", fileName.Data());
1475
1476   if (fPIDParticles && fPIDTracks)
1477   {
1478     TDatabasePDG* pdgDB = new TDatabasePDG;
1479
1480     for (Int_t i=0; i <= fPIDParticles->GetNbinsX()+1; ++i)
1481       if (fPIDParticles->GetBinContent(i) > 0)
1482       {
1483         TObject* pdgParticle = pdgDB->GetParticle((Int_t) fPIDParticles->GetBinCenter(i));
1484         printf("PDG = %d (%s): generated: %d, reconstructed: %d, ratio: %f\n", (Int_t) fPIDParticles->GetBinCenter(i), (pdgParticle) ? pdgParticle->GetName() : "not found", (Int_t) fPIDParticles->GetBinContent(i), (Int_t) fPIDTracks->GetBinContent(i), ((fPIDTracks->GetBinContent(i) > 0) ? fPIDParticles->GetBinContent(i) / fPIDTracks->GetBinContent(i) : -1));
1485       }
1486
1487     delete pdgDB;
1488     pdgDB = 0;
1489   }
1490
1491   fout->Write();
1492   fout->Close();
1493 }
1494
1495 Bool_t AlidNdEtaCorrectionTask::SignOK(TParticlePDG* particle)
1496 {
1497   // returns if a particle with this sign should be counted
1498   // this is determined by the value of fSignMode, which should have the same sign
1499   // as the charge
1500   // if fSignMode is 0 all particles are counted
1501
1502   if (fSignMode == 0)
1503     return kTRUE;
1504
1505   if (!particle)
1506   {
1507     printf("WARNING: not counting a particle that does not have a pdg particle\n");
1508     return kFALSE;
1509   }
1510
1511   Double_t charge = particle->Charge();
1512
1513   if (fSignMode > 0)
1514     if (charge < 0)
1515       return kFALSE;
1516
1517   if (fSignMode < 0)
1518     if (charge > 0)
1519       return kFALSE;
1520
1521   return kTRUE;
1522 }
1523
1524 Double_t AlidNdEtaCorrectionTask::GetSecondaryCorrection(Double_t pt){
1525
1526         // getting the data driven correction factor to correct for
1527         // the underestimate of secondaries in Pythia
1528         // (A. Dainese + J. Otwinowski
1529
1530         if (pt <= 0.17) return 1.0;
1531         if (pt <= 0.4) return GetLinearInterpolationValue(0.17,1.0,0.4,1.07, pt);
1532         if (pt <= 0.6) return GetLinearInterpolationValue(0.4,1.07,0.6,1.25, pt);
1533         if (pt <= 1.2) return GetLinearInterpolationValue(0.6,1.25,1.2,1.5,  pt);
1534         return 1.5;
1535
1536 }
1537
1538 Double_t AlidNdEtaCorrectionTask::GetLinearInterpolationValue(Double_t const x1, Double_t const y1, Double_t const x2, Double_t const y2, const Double_t pt)
1539 {
1540
1541         //
1542         // linear interpolation to be used to get the secondary correction (see AlidNdEtaCorrectionTask::GetSecondaryCorrection)
1543         //
1544
1545         return ((y2-y1)/(x2-x1))*pt+(y2-(((y2-y1)/(x2-x1))*x2));
1546 }
1547