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