PWG0/dNdEta/AlidNdEtaTask.cxx

   1 /* $Id$ */
   2
   3 #include "AlidNdEtaTask.h"
   4
   5 #include <TStyle.h>
   6 #include <TSystem.h>
   7 #include <TCanvas.h>
   8 #include <TVector3.h>
   9 #include <TChain.h>
  10 #include <TFile.h>
  11 #include <TH1F.h>
  12 #include <TH2F.h>
  13 #include <TH3F.h>
  14 #include <TParticle.h>
  15 #include <TRandom.h>
  16 #include <TNtuple.h>
  17 #include <TObjString.h>
  18 #include <TF1.h>
  19 #include <TGraph.h>
  20
  21 #include <AliLog.h>
  22 #include <AliESDVertex.h>
  23 #include <AliESDEvent.h>
  24 #include <AliStack.h>
  25 #include <AliHeader.h>
  26 #include <AliGenEventHeader.h>
  27 #include <AliMultiplicity.h>
  28 #include <AliAnalysisManager.h>
  29 #include <AliMCEventHandler.h>
  30 #include <AliMCEvent.h>
  31 #include <AliESDInputHandler.h>
  32 #include <AliESDInputHandlerRP.h>
  33 #include <AliESDHeader.h>
  34
  35 #include "AliESDtrackCuts.h"
  36 #include "AliPWG0Helper.h"
  37 #include "AliCorrection.h"
  38 #include "AliCorrectionMatrix3D.h"
  39 #include "dNdEta/dNdEtaAnalysis.h"
  40 #include "AliTriggerAnalysis.h"
  41 #include "AliPhysicsSelection.h"
  42
  43 //#define FULLALIROOT
  44
  45 #ifdef FULLALIROOT
  46   #include "../ITS/AliITSRecPoint.h"
  47   #include "AliCDBManager.h"
  48   #include "AliCDBEntry.h"
  49   #include "AliGeomManager.h"
  50   #include "TGeoManager.h"
  51 #endif
  52
  53 ClassImp(AlidNdEtaTask)
  54
  55 AlidNdEtaTask::AlidNdEtaTask(const char* opt) :
  56   AliAnalysisTaskSE("AlidNdEtaTask"),
  57   fESD(0),
  58   fOutput(0),
  59   fOption(opt),
  60   fAnalysisMode((AliPWG0Helper::AnalysisMode) (AliPWG0Helper::kTPC | AliPWG0Helper::kFieldOn)),
  61   fTrigger(AliTriggerAnalysis::kMB1),
  62   fRequireTriggerClass(),
  63   fRejectTriggerClass(),
  64   fFillPhi(kFALSE),
  65   fDeltaPhiCut(-1),
  66   fReadMC(kFALSE),
  67   fUseMCVertex(kFALSE),
  68   fOnlyPrimaries(kFALSE),
  69   fUseMCKine(kFALSE),
  70   fCheckEventType(kFALSE),
  71   fSymmetrize(kFALSE),
  72   fMultAxisEta1(kFALSE),
  73   fDiffTreatment(AliPWG0Helper::kMCFlags),
  74   fEsdTrackCuts(0),
  75   fdNdEtaAnalysisESD(0),
  76   fMult(0),
  77   fMultVtx(0),
  78   fEvents(0),
  79   fVertexResolution(0),
  80   fdNdEtaAnalysis(0),
  81   fdNdEtaAnalysisND(0),
  82   fdNdEtaAnalysisNSD(0),
  83   fdNdEtaAnalysisOnePart(0),
  84   fdNdEtaAnalysisTr(0),
  85   fdNdEtaAnalysisTrVtx(0),
  86   fdNdEtaAnalysisTracks(0),
  87   fPartPt(0),
  88   fVertex(0),
  89   fVertexVsMult(0),
  90   fPhi(0),
  91   fRawPt(0),
  92   fEtaPhi(0),
  93   fModuleMap(0),
  94   fDeltaPhi(0),
  95   fDeltaTheta(0),
  96   fFiredChips(0),
  97   fTrackletsVsClusters(0),
  98   fTrackletsVsUnassigned(0),
  99   fStats(0),
 100   fStats2(0)
 101 {
 102   //
 103   // Constructor. Initialization of pointers
 104   //
 105
 106   // Define input and output slots here
 107   //DefineInput(0, TChain::Class());
 108   DefineOutput(1, TList::Class());
 109
 110   fZPhi[0] = 0;
 111   fZPhi[1] = 0;
 112
 113   AliLog::SetClassDebugLevel("AlidNdEtaTask", AliLog::kWarning);
 114 }
 115
 116 AlidNdEtaTask::~AlidNdEtaTask()
 117 {
 118   //
 119   // Destructor
 120   //
 121
 122   // histograms are in the output list and deleted when the output
 123   // list is deleted by the TSelector dtor
 124
 125   if (fOutput) {
 126     delete fOutput;
 127     fOutput = 0;
 128   }
 129 }
 130
 131 Bool_t AlidNdEtaTask::UserNotify()
 132 {
 133   static Int_t count = 0;
 134   count++;
 135   Printf("Processing %d. file: %s", count, ((TTree*) GetInputData(0))->GetCurrentFile()->GetName());
 136   return kTRUE;
 137 }
 138
 139 //________________________________________________________________________
 140 void AlidNdEtaTask::ConnectInputData(Option_t *opt)
 141 {
 142   // Connect ESD
 143   // Called once
 144
 145   AliAnalysisTaskSE::ConnectInputData(opt);
 146
 147   Printf("AlidNdEtaTask::ConnectInputData called");
 148
 149   /*
 150   AliESDInputHandler *esdH = dynamic_cast<AliESDInputHandler*> (AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler());
 151
 152   if (!esdH) {
 153     Printf("ERROR: Could not get ESDInputHandler");
 154   } else {
 155     fESD = esdH->GetEvent();
 156
 157     TString branches("AliESDHeader Vertex ");
 158
 159     if (fAnalysisMode & AliPWG0Helper::kSPD || fTrigger & AliTriggerAnalysis::kOfflineFlag)
 160       branches += "AliMultiplicity ";
 161
 162     if (fAnalysisMode & AliPWG0Helper::kTPC || fAnalysisMode & AliPWG0Helper::kTPCITS)
 163       branches += "Tracks ";
 164
 165     // Enable only the needed branches
 166     esdH->SetActiveBranches(branches);
 167   }
 168   */
 169
 170   // disable info messages of AliMCEvent (per event)
 171   AliLog::SetClassDebugLevel("AliMCEvent", AliLog::kWarning - AliLog::kDebug + 1);
 172
 173   #ifdef FULLALIROOT
 174     if (fCheckEventType)
 175       AliCDBManager::Instance()->SetDefaultStorage("raw://");
 176     else
 177       AliCDBManager::Instance()->SetDefaultStorage("MC", "Residual");
 178     AliCDBManager::Instance()->SetRun(0);
 179
 180     AliCDBManager* mgr = AliCDBManager::Instance();
 181     AliCDBEntry* obj = mgr->Get(AliCDBPath("GRP", "Geometry", "Data"));
 182     AliGeomManager::SetGeometry((TGeoManager*) obj->GetObject());
 183
 184     AliGeomManager::GetNalignable("ITS");
 185     AliGeomManager::ApplyAlignObjsFromCDB("ITS");
 186   #endif
 187 }
 188
 189 void AlidNdEtaTask::UserCreateOutputObjects()
 190 {
 191   // create result objects and add to output list
 192
 193   Printf("AlidNdEtaTask::CreateOutputObjects");
 194
 195   if (fOnlyPrimaries)
 196     Printf("WARNING: Processing only primaries (MC information used). This is for systematical checks only.");
 197
 198   if (fUseMCKine)
 199     Printf("WARNING: Using MC kine information. This is for systematical checks only.");
 200
 201   if (fUseMCVertex)
 202     Printf("WARNING: Replacing vertex by MC vertex. This is for systematical checks only.");
 203
 204   fOutput = new TList;
 205   fOutput->SetOwner();
 206
 207   fdNdEtaAnalysisESD = new dNdEtaAnalysis("fdNdEtaAnalysisESD", "fdNdEtaAnalysisESD", fAnalysisMode);
 208   fOutput->Add(fdNdEtaAnalysisESD);
 209
 210   fMult = new TH1F("fMult", "fMult;Ntracks;Count", 201, -0.5, 200.5);
 211   fOutput->Add(fMult);
 212
 213   fMultVtx = new TH1F("fMultVtx", "fMultVtx;Ntracks;Count", 201, -0.5, 200.5);
 214   fOutput->Add(fMultVtx);
 215
 216   for (Int_t i=0; i<3; ++i)
 217   {
 218     fPartEta[i] = new TH1F(Form("dndeta_check_%d", i), Form("dndeta_check_%d", i), 60, -3, 3);
 219     fPartEta[i]->Sumw2();
 220     fOutput->Add(fPartEta[i]);
 221   }
 222
 223   fEvents = new TH1F("dndeta_check_vertex", "dndeta_check_vertex", 800, -40, 40);
 224   fOutput->Add(fEvents);
 225
 226   fVertexResolution = new TH2F("dndeta_vertex_resolution_z", ";z resolution;#Delta #phi", 1000, 0, 2, 100, 0, 0.2);
 227   fOutput->Add(fVertexResolution);
 228
 229   fPhi = new TH1F("fPhi", "fPhi;#phi in rad.;count", 720, 0, 2 * TMath::Pi());
 230   fOutput->Add(fPhi);
 231
 232   fEtaPhi = new TH2F("fEtaPhi", "fEtaPhi;#eta;#phi in rad.;count", 80, -4, 4, 18*20, 0, 2 * TMath::Pi());
 233   fOutput->Add(fEtaPhi);
 234
 235   fStats = new TH1F("fStats", "fStats", 5, 0.5, 5.5);
 236   fStats->GetXaxis()->SetBinLabel(1, "vertexer 3d");
 237   fStats->GetXaxis()->SetBinLabel(2, "vertexer z");
 238   fStats->GetXaxis()->SetBinLabel(3, "trigger");
 239   fStats->GetXaxis()->SetBinLabel(4, "physics events");
 240   fStats->GetXaxis()->SetBinLabel(5, "physics events after veto");
 241   fOutput->Add(fStats);
 242
 243   fStats2 = new TH2F("fStats2", "fStats2", 7, -0.5, 6.5, 12, -0.5, 11.5);
 244
 245   fStats2->GetXaxis()->SetBinLabel(1, "No trigger");
 246   fStats2->GetXaxis()->SetBinLabel(2, "No Vertex");
 247   fStats2->GetXaxis()->SetBinLabel(3, "Vtx res. too bad");
 248   fStats2->GetXaxis()->SetBinLabel(4, "|z-vtx| > 15");
 249   fStats2->GetXaxis()->SetBinLabel(5, "|z-vtx| > 10");
 250   fStats2->GetXaxis()->SetBinLabel(6, "0 tracklets");
 251   fStats2->GetXaxis()->SetBinLabel(7, "Selected");
 252
 253   fStats2->GetYaxis()->SetBinLabel(1, "n/a");
 254   fStats2->GetYaxis()->SetBinLabel(2, "empty");
 255   fStats2->GetYaxis()->SetBinLabel(3, "BBA");
 256   fStats2->GetYaxis()->SetBinLabel(4, "BBC");
 257   fStats2->GetYaxis()->SetBinLabel(5, "BBA BBC");
 258   fStats2->GetYaxis()->SetBinLabel(6, "BGA");
 259   fStats2->GetYaxis()->SetBinLabel(7, "BGC");
 260   fStats2->GetYaxis()->SetBinLabel(8, "BGA BGC");
 261   fStats2->GetYaxis()->SetBinLabel(9, "BBA BGC");
 262   fStats2->GetYaxis()->SetBinLabel(10, "BGA BBC");
 263   fStats2->GetYaxis()->SetBinLabel(11, "GFO");
 264   fStats2->GetYaxis()->SetBinLabel(12, "NO GFO");
 265   fOutput->Add(fStats2);
 266
 267   fTrackletsVsClusters = new TH2F("fTrackletsVsClusters", ";tracklets;clusters in ITS", 50, -0.5, 49.5, 1000, -0.5, 999.5);
 268   fOutput->Add(fTrackletsVsClusters);
 269
 270   if (fAnalysisMode & AliPWG0Helper::kSPD)
 271   {
 272     fDeltaPhi = new TH1F("fDeltaPhi", "fDeltaPhi;#Delta #phi;Entries", 500, -0.16, 0.16);
 273     fOutput->Add(fDeltaPhi);
 274     fDeltaTheta = new TH1F("fDeltaTheta", "fDeltaTheta;#Delta #theta;Entries", 500, -0.05, 0.05);
 275     fOutput->Add(fDeltaTheta);
 276     fFiredChips = new TH2F("fFiredChips", "fFiredChips;Chips L1 + L2;tracklets", 1201, -0.5, 1201, 50, -0.5, 49.5);
 277     fOutput->Add(fFiredChips);
 278     fTrackletsVsUnassigned = new TH2F("fTrackletsVsUnassigned", ";tracklets;unassigned clusters in L0", 50, -0.5, 49.5, 200, -0.5, 199.5);
 279     fOutput->Add(fTrackletsVsUnassigned);
 280     for (Int_t i=0; i<2; i++)
 281     {
 282       fZPhi[i] = new TH2F(Form("fZPhi_%d", i), Form("fZPhi Layer %d;z (cm);#phi (rad.)", i), 200, -15, 15, 200, 0, TMath::Pi() * 2);
 283       fOutput->Add(fZPhi[i]);
 284     }
 285
 286     fModuleMap = new TH1F("fModuleMap", "fModuleMap;module number;cluster count", 240, -0.5, 239.5);
 287     fOutput->Add(fModuleMap);
 288   }
 289
 290   if (fAnalysisMode & AliPWG0Helper::kTPC || fAnalysisMode & AliPWG0Helper::kTPCITS)
 291   {
 292     fRawPt =  new TH1F("fRawPt", "raw pt;p_{T};Count", 2000, 0, 100);
 293     fOutput->Add(fRawPt);
 294   }
 295
 296   fVertex = new TH3F("vertex_check", "vertex_check;x;y;z", 100, -1, 1, 100, -1, 1, 100, -30, 30);
 297   fOutput->Add(fVertex);
 298
 299   fVertexVsMult = new TH3F("fVertexVsMult", "fVertexVsMult;x;y;multiplicity", 100, -1, 1, 100, -1, 1, 100, -0.5, 99.5);
 300   fOutput->Add(fVertexVsMult);
 301
 302   if (fReadMC)
 303   {
 304     fdNdEtaAnalysis = new dNdEtaAnalysis("dndeta", "dndeta", fAnalysisMode);
 305     fOutput->Add(fdNdEtaAnalysis);
 306
 307     fdNdEtaAnalysisND = new dNdEtaAnalysis("dndetaND", "dndetaND", fAnalysisMode);
 308     fOutput->Add(fdNdEtaAnalysisND);
 309
 310     fdNdEtaAnalysisNSD = new dNdEtaAnalysis("dndetaNSD", "dndetaNSD", fAnalysisMode);
 311     fOutput->Add(fdNdEtaAnalysisNSD);
 312
 313     fdNdEtaAnalysisOnePart = new dNdEtaAnalysis("dndetaOnePart", "dndetaOnePart", fAnalysisMode);
 314     fOutput->Add(fdNdEtaAnalysisOnePart);
 315
 316     fdNdEtaAnalysisTr = new dNdEtaAnalysis("dndetaTr", "dndetaTr", fAnalysisMode);
 317     fOutput->Add(fdNdEtaAnalysisTr);
 318
 319     fdNdEtaAnalysisTrVtx = new dNdEtaAnalysis("dndetaTrVtx", "dndetaTrVtx", fAnalysisMode);
 320     fOutput->Add(fdNdEtaAnalysisTrVtx);
 321
 322     fdNdEtaAnalysisTracks = new dNdEtaAnalysis("dndetaTracks", "dndetaTracks", fAnalysisMode);
 323     fOutput->Add(fdNdEtaAnalysisTracks);
 324
 325     fPartPt =  new TH1F("dndeta_check_pt", "dndeta_check_pt", 1000, 0, 10);
 326     fPartPt->Sumw2();
 327     fOutput->Add(fPartPt);
 328   }
 329
 330   if (fEsdTrackCuts)
 331   {
 332     fEsdTrackCuts->SetName("fEsdTrackCuts");
 333     fOutput->Add(fEsdTrackCuts);
 334   }
 335
 336   //AliLog::SetClassDebugLevel("AliPhysicsSelection", AliLog::kDebug);
 337 }
 338
 339 Bool_t AlidNdEtaTask::IsEventInBinZero()
 340 {
 341   // checks if the event goes to the 0 bin
 342
 343   fESD = (AliESDEvent*) fInputEvent;
 344
 345   AliInputEventHandler* inputHandler = static_cast<AliInputEventHandler*> (AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler());
 346   if (!inputHandler)
 347   {
 348     Printf("ERROR: Could not receive input handler");
 349     return kFALSE;
 350   }
 351
 352   static AliTriggerAnalysis* triggerAnalysis = 0;
 353   if (!triggerAnalysis)
 354   {
 355     AliPhysicsSelection* physicsSelection = dynamic_cast<AliPhysicsSelection*> (inputHandler->GetEventSelection());
 356     if (physicsSelection)
 357       triggerAnalysis = physicsSelection->GetTriggerAnalysis();
 358   }
 359
 360   if (!triggerAnalysis)
 361   {
 362     Printf("ERROR: Could not receive trigger analysis object");
 363     return kFALSE;
 364   }
 365
 366   if (!triggerAnalysis->IsTriggerFired(fESD, fTrigger))
 367     return kFALSE;
 368
 369   // TODO 0 bin check
 370
 371   return kTRUE;
 372 }
 373
 374 void AlidNdEtaTask::UserExec(Option_t*)
 375 {
 376   // process the event
 377
 378   fESD = (AliESDEvent*) fInputEvent;
 379
 380   // these variables are also used in the MC section below; however, if ESD is off they stay with their default values
 381   Bool_t eventTriggered = kFALSE;
 382   const AliESDVertex* vtxESD = 0;
 383
 384   // post the data already here
 385   PostData(1, fOutput);
 386
 387   // ESD analysis
 388   if (fESD)
 389   {
 390     AliESDHeader* esdHeader = fESD->GetHeader();
 391     if (!esdHeader)
 392     {
 393       Printf("ERROR: esdHeader could not be retrieved");
 394       return;
 395     }
 396
 397     AliInputEventHandler* inputHandler = (AliInputEventHandler*) AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler();
 398     if (!inputHandler)
 399     {
 400       Printf("ERROR: Could not receive input handler");
 401       return;
 402     }
 403
 404     eventTriggered = inputHandler->IsEventSelected();
 405
 406     static AliTriggerAnalysis* triggerAnalysis = 0;
 407     if (!triggerAnalysis)
 408     {
 409       AliPhysicsSelection* physicsSelection = dynamic_cast<AliPhysicsSelection*> (inputHandler->GetEventSelection());
 410       if (physicsSelection)
 411         triggerAnalysis = physicsSelection->GetTriggerAnalysis();
 412     }
 413
 414     if (eventTriggered)
 415       eventTriggered = triggerAnalysis->IsTriggerFired(fESD, fTrigger);
 416
 417     AliTriggerAnalysis::V0Decision v0A = triggerAnalysis->V0Trigger(fESD, AliTriggerAnalysis::kASide, kFALSE);
 418     AliTriggerAnalysis::V0Decision v0C = triggerAnalysis->V0Trigger(fESD, AliTriggerAnalysis::kCSide, kFALSE);
 419     Bool_t fastORHW = (triggerAnalysis->SPDFiredChips(fESD, 1) > 0);
 420
 421     Int_t vZero = 0;
 422     if (v0A != AliTriggerAnalysis::kV0Invalid && v0C != AliTriggerAnalysis::kV0Invalid)
 423     {
 424       if (v0A == AliTriggerAnalysis::kV0Empty && v0C == AliTriggerAnalysis::kV0Empty) vZero = 1;
 425       if (v0A == AliTriggerAnalysis::kV0BB    && v0C == AliTriggerAnalysis::kV0Empty) vZero = 2;
 426       if (v0A == AliTriggerAnalysis::kV0Empty && v0C == AliTriggerAnalysis::kV0BB)    vZero = 3;
 427       if (v0A == AliTriggerAnalysis::kV0BB    && v0C == AliTriggerAnalysis::kV0BB)    vZero = 4;
 428       if (v0A == AliTriggerAnalysis::kV0BG    && v0C == AliTriggerAnalysis::kV0Empty) vZero = 5;
 429       if (v0A == AliTriggerAnalysis::kV0Empty && v0C == AliTriggerAnalysis::kV0BG)    vZero = 6;
 430       if (v0A == AliTriggerAnalysis::kV0BG    && v0C == AliTriggerAnalysis::kV0BG)    vZero = 7;
 431       if (v0A == AliTriggerAnalysis::kV0BB    && v0C == AliTriggerAnalysis::kV0BG)    vZero = 8;
 432       if (v0A == AliTriggerAnalysis::kV0BG    && v0C == AliTriggerAnalysis::kV0BB)    vZero = 9;
 433     }
 434
 435     if (vZero == 0)
 436       Printf("VZERO: %d %d %d %d", vZero, eventTriggered, v0A, v0C);
 437
 438     Bool_t filled = kFALSE;
 439
 440     if (!eventTriggered)
 441     {
 442       fStats2->Fill(0.0, vZero);
 443       fStats2->Fill(0.0, (fastORHW) ? 10 : 11);
 444       filled = kTRUE;
 445     }
 446
 447     if (eventTriggered)
 448       fStats->Fill(3);
 449
 450     /*
 451     // ITS cluster tree
 452     AliESDInputHandlerRP* handlerRP = dynamic_cast<AliESDInputHandlerRP*> (AliAnalysisManager::GetAnalys
 453 isManager()->GetInputEventHandler());
 454     if (handlerRP)
 455     {
 456       TTree* itsClusterTree = handlerRP->GetTreeR("ITS");
 457       if (!itsClusterTree)
 458         return;
 459
 460       TClonesArray* itsClusters = new TClonesArray("AliITSRecPoint");
 461       TBranch* itsClusterBranch=itsClusterTree->GetBranch("ITSRecPoints");
 462
 463       itsClusterBranch->SetAddress(&itsClusters);
 464
 465       Int_t nItsSubs = (Int_t)itsClusterTree->GetEntries();
 466
 467       Int_t totalClusters = 0;
 468
 469       // loop over the its subdetectors
 470       for (Int_t iIts=0; iIts < nItsSubs; iIts++) {
 471
 472         if (!itsClusterTree->GetEvent(iIts))
 473           continue;
 474
 475         Int_t nClusters = itsClusters->GetEntriesFast();
 476         totalClusters += nClusters;
 477
 478         #ifdef FULLALIROOT
 479           if (fAnalysisMode & AliPWG0Helper::kSPD)
 480           {
 481             // loop over clusters
 482             while (nClusters--) {
 483               AliITSRecPoint* cluster = (AliITSRecPoint*) itsClusters->UncheckedAt(nClusters);
 484
 485               Int_t layer = cluster->GetLayer();
 486
 487               if (layer > 1)
 488                 continue;
 489
 490               Float_t xyz[3] = {0., 0., 0.};
 491               cluster->GetGlobalXYZ(xyz);
 492
 493               Float_t phi = TMath::Pi() + TMath::ATan2(-xyz[1], -xyz[0]);
 494               Float_t z = xyz[2];
 495
 496               fZPhi[layer]->Fill(z, phi);
 497               fModuleMap->Fill(layer * 80 + cluster->GetDetectorIndex());
 498             }
 499           }
 500         #endif
 501       }
 502     }
 503     */
 504
 505     vtxESD = AliPWG0Helper::GetVertex(fESD, fAnalysisMode);
 506     if (!vtxESD)
 507     {
 508       if (!filled)
 509       {
 510         fStats2->Fill(1, vZero);
 511         fStats2->Fill(1, (fastORHW) ? 10 : 11);
 512         filled = kTRUE;
 513       }
 514     }
 515     else
 516     {
 517       if (!AliPWG0Helper::TestVertex(vtxESD, fAnalysisMode))
 518       {
 519         if (!filled)
 520         {
 521           fStats2->Fill(2, vZero);
 522           fStats2->Fill(2, (fastORHW) ? 10 : 11);
 523           filled = kTRUE;
 524         }
 525       }
 526
 527       Double_t vtx[3];
 528       vtxESD->GetXYZ(vtx);
 529
 530       // try to compare spd vertex and vertexer from tracks
 531       // remove vertices outside +- 15 cm
 532       if (TMath::Abs(vtx[2]) > 15)
 533       {
 534         if (!filled)
 535         {
 536           fStats2->Fill(3, vZero);
 537           fStats2->Fill(3, (fastORHW) ? 10 : 11);
 538           filled = kTRUE;
 539         }
 540       }
 541
 542       if (TMath::Abs(vtx[2]) > 10)
 543       {
 544         if (!filled)
 545         {
 546           fStats2->Fill(4, vZero);
 547           fStats2->Fill(4, (fastORHW) ? 10 : 11);
 548           filled = kTRUE;
 549         }
 550       }
 551
 552       const AliMultiplicity* mult = fESD->GetMultiplicity();
 553       if (!mult)
 554         return;
 555
 556       if (mult->GetNumberOfTracklets() == 0)
 557       {
 558         if (!filled)
 559         {
 560           fStats2->Fill(5, vZero);
 561           fStats2->Fill(5, (fastORHW) ? 10 : 11);
 562           filled = kTRUE;
 563         }
 564       }
 565     }
 566
 567     if (!filled)
 568     {
 569       fStats2->Fill(6, vZero);
 570       fStats2->Fill(6, (fastORHW) ? 10 : 11);
 571       //Printf("File: %s, IEV: %d, TRG: ---, Orbit: 0x%x, Period: %d, BC: %d", ((TTree*) GetInputData(0))->GetCurrentFile()->GetName(), fESD->GetEventNumberInFile(), fESD->GetOrbitNumber(),fESD->GetPeriodNumber(),fESD->GetBunchCrossNumber());
 572     }
 573
 574     if (eventTriggered)
 575       fStats->Fill(3);
 576
 577     fStats->Fill(5);
 578
 579     // get the ESD vertex
 580     vtxESD = AliPWG0Helper::GetVertex(fESD, fAnalysisMode);
 581
 582     Double_t vtx[3];
 583
 584     // fill z vertex resolution
 585     if (vtxESD)
 586     {
 587       if (strcmp(vtxESD->GetTitle(), "vertexer: Z") == 0)
 588         fVertexResolution->Fill(vtxESD->GetZRes(), vtxESD->GetDispersion());
 589       else
 590         fVertexResolution->Fill(vtxESD->GetZRes(), 0);
 591
 592       //if (strcmp(vtxESD->GetTitle(), "vertexer: 3D") == 0)
 593       {
 594         fVertex->Fill(vtxESD->GetXv(), vtxESD->GetYv(), vtxESD->GetZv());
 595       }
 596
 597       if (AliPWG0Helper::TestVertex(vtxESD, fAnalysisMode))
 598       {
 599           vtxESD->GetXYZ(vtx);
 600
 601           // vertex stats
 602           if (strcmp(vtxESD->GetTitle(), "vertexer: 3D") == 0)
 603           {
 604             fStats->Fill(1);
 605             if (fCheckEventType && TMath::Abs(vtx[0]) > 0.3)
 606             {
 607               Printf("Suspicious x-vertex x=%f y=%f z=%f (period: %d orbit %x)", vtx[0], vtx[1], vtx[2], fESD->GetPeriodNumber(), fESD->GetOrbitNumber());
 608             }
 609             if (fCheckEventType && (vtx[1] < 0.05 || vtx[1] > 0.5))
 610             {
 611               Printf("Suspicious y-vertex x=%f y=%f z=%f (period: %d orbit %x)", vtx[0], vtx[1], vtx[2], fESD->GetPeriodNumber(), fESD->GetOrbitNumber());
 612             }
 613           }
 614
 615           if (strcmp(vtxESD->GetTitle(), "vertexer: Z") == 0)
 616             fStats->Fill(2);
 617
 618           // remove vertices outside +- 15 cm
 619           if (TMath::Abs(vtx[2]) > 15)
 620             vtxESD = 0;
 621       }
 622       else
 623         vtxESD = 0;
 624     }
 625
 626     // needed for syst. studies
 627     AliStack* stack = 0;
 628     TArrayF vtxMC(3);
 629
 630     if (fUseMCVertex || fUseMCKine || fOnlyPrimaries || fReadMC) {
 631       if (!fReadMC) {
 632         Printf("ERROR: fUseMCVertex or fUseMCKine or fOnlyPrimaries set without fReadMC set!");
 633         return;
 634       }
 635
 636       AliMCEventHandler* eventHandler = dynamic_cast<AliMCEventHandler*> (AliAnalysisManager::GetAnalysisManager()->GetMCtruthEventHandler());
 637       if (!eventHandler) {
 638         Printf("ERROR: Could not retrieve MC event handler");
 639         return;
 640       }
 641
 642       AliMCEvent* mcEvent = eventHandler->MCEvent();
 643       if (!mcEvent) {
 644         Printf("ERROR: Could not retrieve MC event");
 645         return;
 646       }
 647
 648       AliHeader* header = mcEvent->Header();
 649       if (!header)
 650       {
 651         AliDebug(AliLog::kError, "Header not available");
 652         return;
 653       }
 654
 655       // get the MC vertex
 656       AliGenEventHeader* genHeader = header->GenEventHeader();
 657       if (!genHeader)
 658       {
 659         AliDebug(AliLog::kError, "Could not retrieve genHeader from Header");
 660         return;
 661       }
 662       genHeader->PrimaryVertex(vtxMC);
 663
 664       if (fUseMCVertex)
 665         vtx[2] = vtxMC[2];
 666
 667       stack = mcEvent->Stack();
 668       if (!stack)
 669       {
 670         AliDebug(AliLog::kError, "Stack not available");
 671         return;
 672       }
 673     }
 674
 675     // create list of (label, eta, pt) tuples
 676     Int_t inputCount = 0;
 677     Int_t* labelArr = 0;
 678     Float_t* etaArr = 0;
 679     Float_t* thirdDimArr = 0;
 680     if (fAnalysisMode & AliPWG0Helper::kSPD)
 681     {
 682       if (vtxESD)
 683       {
 684         // get tracklets
 685         const AliMultiplicity* mult = fESD->GetMultiplicity();
 686         if (!mult)
 687         {
 688           AliDebug(AliLog::kError, "AliMultiplicity not available");
 689           return;
 690         }
 691
 692         Int_t arrayLength = mult->GetNumberOfTracklets();
 693         if (fAnalysisMode & AliPWG0Helper::kSPDOnlyL0)
 694           arrayLength += mult->GetNumberOfSingleClusters();
 695
 696         labelArr = new Int_t[arrayLength];
 697         etaArr = new Float_t[arrayLength];
 698         thirdDimArr = new Float_t[arrayLength];
 699
 700         // get multiplicity from SPD tracklets
 701         for (Int_t i=0; i<mult->GetNumberOfTracklets(); ++i)
 702         {
 703           //printf("%d %f %f %f\n", i, mult->GetTheta(i), mult->GetPhi(i), mult->GetDeltaPhi(i));
 704
 705           if (fOnlyPrimaries)
 706             if (mult->GetLabel(i, 0) < 0 || mult->GetLabel(i, 0) != mult->GetLabel(i, 1) || !stack->IsPhysicalPrimary(mult->GetLabel(i, 0)))
 707               continue;
 708
 709           Float_t deltaPhi = mult->GetDeltaPhi(i);
 710           // prevent values to be shifted by 2 Pi()
 711           if (deltaPhi < -TMath::Pi())
 712             deltaPhi += TMath::Pi() * 2;
 713           if (deltaPhi > TMath::Pi())
 714             deltaPhi -= TMath::Pi() * 2;
 715
 716           if (TMath::Abs(deltaPhi) > 1)
 717             printf("WARNING: Very high Delta Phi: %d %f %f %f\n", i, mult->GetTheta(i), mult->GetPhi(i), deltaPhi);
 718
 719           Int_t label = mult->GetLabel(i, 0);
 720           Float_t eta = mult->GetEta(i);
 721
 722           // control histograms
 723           Float_t phi = mult->GetPhi(i);
 724           if (phi < 0)
 725             phi += TMath::Pi() * 2;
 726
 727           // TEST exclude probably inefficient phi region
 728           //if (phi > 5.70 || phi < 0.06)
 729           //  continue;
 730
 731           fPhi->Fill(phi);
 732
 733           if (vtxESD && TMath::Abs(vtx[2]) < 10)
 734           {
 735             fEtaPhi->Fill(eta, phi);
 736             fDeltaPhi->Fill(deltaPhi);
 737             fDeltaTheta->Fill(mult->GetDeltaTheta(i));
 738           }
 739
 740           if (fDeltaPhiCut > 0 && (TMath::Abs(deltaPhi) > fDeltaPhiCut || TMath::Abs(mult->GetDeltaTheta(i)) > fDeltaPhiCut / 0.08 * 0.025))
 741             continue;
 742
 743           if (fUseMCKine)
 744           {
 745             if (label > 0)
 746             {
 747               TParticle* particle = stack->Particle(label);
 748               eta = particle->Eta();
 749               phi = particle->Phi();
 750             }
 751             else
 752               Printf("WARNING: fUseMCKine set without fOnlyPrimaries and no label found");
 753           }
 754
 755           if (fSymmetrize)
 756             eta = TMath::Abs(eta);
 757
 758           etaArr[inputCount] = eta;
 759           labelArr[inputCount] = label;
 760           thirdDimArr[inputCount] = phi;
 761           ++inputCount;
 762         }
 763
 764         if (fAnalysisMode & AliPWG0Helper::kSPDOnlyL0)
 765         {
 766           // get additional clusters from L0
 767           for (Int_t i=0; i<mult->GetNumberOfSingleClusters(); ++i)
 768           {
 769             etaArr[inputCount] = -TMath::Log(TMath::Tan(mult->GetThetaSingle(i)/2.));
 770             labelArr[inputCount] = -1;
 771             thirdDimArr[inputCount] = mult->GetPhiSingle(i);
 772
 773             ++inputCount;
 774           }
 775         }
 776
 777         if (!fFillPhi)
 778         {
 779           // fill multiplicity in third axis
 780           for (Int_t i=0; i<inputCount; ++i)
 781             thirdDimArr[i] = inputCount;
 782         }
 783
 784         Int_t firedChips = mult->GetNumberOfFiredChips(0) + mult->GetNumberOfFiredChips(1);
 785         fFiredChips->Fill(firedChips, inputCount);
 786         Printf("Accepted %d tracklets (%d fired chips)", inputCount, firedChips);
 787
 788         fTrackletsVsUnassigned->Fill(inputCount, mult->GetNumberOfSingleClusters());
 789       }
 790     }
 791     else if (fAnalysisMode & AliPWG0Helper::kTPC || fAnalysisMode & AliPWG0Helper::kTPCITS)
 792     {
 793       if (!fEsdTrackCuts)
 794       {
 795         AliDebug(AliLog::kError, "fESDTrackCuts not available");
 796         return;
 797       }
 798
 799       Bool_t foundInEta10 = kFALSE;
 800
 801       if (vtxESD)
 802       {
 803         // get multiplicity from ESD tracks
 804         TObjArray* list = fEsdTrackCuts->GetAcceptedTracks(fESD, fAnalysisMode & AliPWG0Helper::kTPC);
 805         Int_t nGoodTracks = list->GetEntries();
 806         Printf("Accepted %d tracks out of %d total ESD tracks", nGoodTracks, fESD->GetNumberOfTracks());
 807
 808         labelArr = new Int_t[nGoodTracks];
 809         etaArr = new Float_t[nGoodTracks];
 810         thirdDimArr = new Float_t[nGoodTracks];
 811
 812         // loop over esd tracks
 813         for (Int_t i=0; i<nGoodTracks; i++)
 814         {
 815           AliESDtrack* esdTrack = dynamic_cast<AliESDtrack*> (list->At(i));
 816           if (!esdTrack)
 817           {
 818             AliDebug(AliLog::kError, Form("ERROR: Could not retrieve track %d.", i));
 819             continue;
 820           }
 821
 822           Float_t phi = esdTrack->Phi();
 823           if (phi < 0)
 824             phi += TMath::Pi() * 2;
 825
 826           Float_t eta = esdTrack->Eta();
 827           Int_t label = TMath::Abs(esdTrack->GetLabel());
 828           Float_t pT  = esdTrack->Pt();
 829
 830           // force pT to fixed value without B field
 831           if (!(fAnalysisMode & AliPWG0Helper::kFieldOn))
 832             pT = 1;
 833
 834           fPhi->Fill(phi);
 835           fEtaPhi->Fill(eta, phi);
 836           if (eventTriggered && vtxESD)
 837             fRawPt->Fill(pT);
 838
 839           if (esdTrack->Pt() < 0.15)
 840             continue;
 841
 842           // INEL>0 trigger
 843           if (TMath::Abs(esdTrack->Eta()) < 1 && esdTrack->Pt() > 0.15)
 844             foundInEta10 = kTRUE;
 845
 846           if (fOnlyPrimaries)
 847           {
 848             if (label == 0)
 849               continue;
 850
 851             if (stack->IsPhysicalPrimary(label) == kFALSE)
 852               continue;
 853           }
 854
 855           if (fUseMCKine)
 856           {
 857             if (label > 0)
 858             {
 859               TParticle* particle = stack->Particle(label);
 860               eta = particle->Eta();
 861               pT = particle->Pt();
 862             }
 863             else
 864               Printf("WARNING: fUseMCKine set without fOnlyPrimaries and no label found");
 865           }
 866
 867           if (fSymmetrize)
 868             eta = TMath::Abs(eta);
 869           etaArr[inputCount] = eta;
 870           labelArr[inputCount] = TMath::Abs(esdTrack->GetLabel());
 871           thirdDimArr[inputCount] = pT;
 872           ++inputCount;
 873         }
 874
 875         //if (inputCount > 30)
 876         //  Printf("Event with %d accepted TPC tracks. Period number: %d Orbit number: %x Bunch crossing number: %d", inputCount, fESD->GetPeriodNumber(), fESD->GetOrbitNumber(), fESD->GetBunchCrossNumber());
 877
 878         // TODO restrict inputCount used as measure for the multiplicity to |eta| < 1
 879
 880         delete list;
 881       }
 882
 883       if (!foundInEta10)
 884         eventTriggered = kFALSE;
 885     }
 886     else
 887       return;
 888
 889     // Processing of ESD information (always)
 890     if (eventTriggered)
 891     {
 892       // control hist
 893       fMult->Fill(inputCount);
 894       fdNdEtaAnalysisESD->FillTriggeredEvent(inputCount);
 895
 896       if (vtxESD)
 897       {
 898         // control hist
 899
 900         if (strcmp(vtxESD->GetTitle(), "vertexer: 3D") == 0)
 901           fVertexVsMult->Fill(vtxESD->GetXv(), vtxESD->GetYv(), inputCount);
 902
 903         Int_t part05 = 0;
 904         Int_t part10 = 0;
 905         for (Int_t i=0; i<inputCount; ++i)
 906         {
 907           Float_t eta = etaArr[i];
 908           Float_t thirdDim = thirdDimArr[i];
 909
 910           fdNdEtaAnalysisESD->FillTrack(vtx[2], eta, thirdDim);
 911
 912           if (TMath::Abs(vtx[2]) < 10)
 913           {
 914             fPartEta[0]->Fill(eta);
 915
 916             if (vtx[2] < 0)
 917               fPartEta[1]->Fill(eta);
 918             else
 919               fPartEta[2]->Fill(eta);
 920           }
 921
 922           if (TMath::Abs(eta) < 0.5)
 923             part05++;
 924           if (TMath::Abs(eta) < 1.0)
 925             part10++;
 926         }
 927
 928         Int_t multAxis = inputCount;
 929         if (fMultAxisEta1)
 930           multAxis = part10;
 931
 932         fMultVtx->Fill(multAxis);
 933         //if (TMath::Abs(vtx[2]) < 10)
 934         //  fMultVtx->Fill(part05);
 935
 936         // for event count per vertex
 937         fdNdEtaAnalysisESD->FillEvent(vtx[2], multAxis);
 938
 939         // control hist
 940         if (multAxis > 0)
 941           fEvents->Fill(vtx[2]);
 942
 943         if (fReadMC)
 944         {
 945           // from tracks is only done for triggered and vertex reconstructed events
 946           for (Int_t i=0; i<inputCount; ++i)
 947           {
 948             Int_t label = labelArr[i];
 949
 950             if (label < 0)
 951               continue;
 952
 953             //Printf("Getting particle of track %d", label);
 954             TParticle* particle = stack->Particle(label);
 955
 956             if (!particle)
 957             {
 958               AliDebug(AliLog::kError, Form("ERROR: Could not retrieve particle %d.", label));
 959               continue;
 960             }
 961
 962             Float_t thirdDim = -1;
 963             if (fAnalysisMode & AliPWG0Helper::kSPD)
 964             {
 965               if (fFillPhi)
 966               {
 967                 thirdDim = particle->Phi();
 968               }
 969               else
 970                 thirdDim = multAxis;
 971             }
 972             else
 973               thirdDim = particle->Pt();
 974
 975             Float_t eta = particle->Eta();
 976             if (fSymmetrize)
 977               eta = TMath::Abs(eta);
 978             fdNdEtaAnalysisTracks->FillTrack(vtxMC[2], eta, thirdDim);
 979           } // end of track loop
 980
 981           // for event count per vertex
 982           fdNdEtaAnalysisTracks->FillEvent(vtxMC[2], multAxis);
 983         }
 984       }
 985     }
 986     if (etaArr)
 987       delete[] etaArr;
 988     if (labelArr)
 989       delete[] labelArr;
 990     if (thirdDimArr)
 991       delete[] thirdDimArr;
 992   }
 993
 994   if (fReadMC)   // Processing of MC information (optional)
 995   {
 996     AliMCEventHandler* eventHandler = dynamic_cast<AliMCEventHandler*> (AliAnalysisManager::GetAnalysisManager()->GetMCtruthEventHandler());
 997     if (!eventHandler) {
 998       Printf("ERROR: Could not retrieve MC event handler");
 999       return;
1000     }
1001
1002     AliMCEvent* mcEvent = eventHandler->MCEvent();
1003     if (!mcEvent) {
1004       Printf("ERROR: Could not retrieve MC event");
1005       return;
1006     }
1007
1008     AliStack* stack = mcEvent->Stack();
1009     if (!stack)
1010     {
1011       AliDebug(AliLog::kError, "Stack not available");
1012       return;
1013     }
1014
1015     AliHeader* header = mcEvent->Header();
1016     if (!header)
1017     {
1018       AliDebug(AliLog::kError, "Header not available");
1019       return;
1020     }
1021
1022     // get the MC vertex
1023     AliGenEventHeader* genHeader = header->GenEventHeader();
1024     if (!genHeader)
1025     {
1026       AliDebug(AliLog::kError, "Could not retrieve genHeader from Header");
1027       return;
1028     }
1029
1030     TArrayF vtxMC(3);
1031     genHeader->PrimaryVertex(vtxMC);
1032
1033     // get process type
1034     Int_t processType = AliPWG0Helper::GetEventProcessType(fESD, header, stack, fDiffTreatment);
1035     AliDebug(AliLog::kDebug+1, Form("Found process type %d", processType));
1036
1037     if (processType == AliPWG0Helper::kInvalidProcess)
1038       AliDebug(AliLog::kError, Form("Unknown process type %d.", processType));
1039
1040     // loop over mc particles
1041     Int_t nPrim  = stack->GetNprimary();
1042
1043     Int_t nAcceptedParticles = 0;
1044     Bool_t oneParticleEvent = kFALSE;
1045
1046     // count particles first, then fill
1047     for (Int_t iMc = 0; iMc < nPrim; ++iMc)
1048     {
1049       if (!stack->IsPhysicalPrimary(iMc))
1050         continue;
1051
1052       //Printf("Getting particle %d", iMc);
1053       TParticle* particle = stack->Particle(iMc);
1054
1055       if (!particle)
1056         continue;
1057
1058       if (particle->GetPDG()->Charge() == 0)
1059         continue;
1060
1061       //AliDebug(AliLog::kDebug+1, Form("Accepted primary %d, unique ID: %d", iMc, particle->GetUniqueID()));
1062       Float_t eta = particle->Eta();
1063
1064       if (TMath::Abs(eta) < 1.0)
1065         oneParticleEvent = kTRUE;
1066
1067       // make a rough eta cut (so that nAcceptedParticles is not too far off the true measured value (NB: these histograms are only gathered for comparison))
1068       if (TMath::Abs(eta) < 1.5) // && pt > 0.3)
1069         nAcceptedParticles++;
1070     }
1071
1072     for (Int_t iMc = 0; iMc < nPrim; ++iMc)
1073     {
1074       if (!stack->IsPhysicalPrimary(iMc))
1075         continue;
1076
1077       //Printf("Getting particle %d", iMc);
1078       TParticle* particle = stack->Particle(iMc);
1079
1080       if (!particle)
1081         continue;
1082
1083       if (particle->GetPDG()->Charge() == 0)
1084         continue;
1085
1086       Float_t eta = particle->Eta();
1087       if (fSymmetrize)
1088         eta = TMath::Abs(eta);
1089
1090       Float_t thirdDim = -1;
1091
1092       if (fAnalysisMode & AliPWG0Helper::kSPD)
1093       {
1094         if (fFillPhi)
1095         {
1096           thirdDim = particle->Phi();
1097         }
1098         else
1099           thirdDim = nAcceptedParticles;
1100       }
1101       else
1102         thirdDim = particle->Pt();
1103
1104       fdNdEtaAnalysis->FillTrack(vtxMC[2], eta, thirdDim);
1105
1106       if (processType != AliPWG0Helper::kSD)
1107         fdNdEtaAnalysisNSD->FillTrack(vtxMC[2], eta, thirdDim);
1108
1109       if (processType == AliPWG0Helper::kND)
1110         fdNdEtaAnalysisND->FillTrack(vtxMC[2], eta, thirdDim);
1111
1112       if (oneParticleEvent)
1113         fdNdEtaAnalysisOnePart->FillTrack(vtxMC[2], eta, thirdDim);
1114
1115       if (eventTriggered)
1116       {
1117         fdNdEtaAnalysisTr->FillTrack(vtxMC[2], eta, thirdDim);
1118         if (vtxESD)
1119           fdNdEtaAnalysisTrVtx->FillTrack(vtxMC[2], eta, thirdDim);
1120       }
1121
1122       if (TMath::Abs(eta) < 1.0 && particle->Pt() > 0 && particle->P() > 0)
1123       {
1124         //Float_t value = 1. / TMath::TwoPi() / particle->Pt() * particle->Energy() / particle->P();
1125         Float_t value = 1;
1126         fPartPt->Fill(particle->Pt(), value);
1127       }
1128     }
1129
1130     fdNdEtaAnalysis->FillEvent(vtxMC[2], nAcceptedParticles);
1131     if (processType != AliPWG0Helper::kSD)
1132       fdNdEtaAnalysisNSD->FillEvent(vtxMC[2], nAcceptedParticles);
1133     if (processType == AliPWG0Helper::kND)
1134       fdNdEtaAnalysisND->FillEvent(vtxMC[2], nAcceptedParticles);
1135     if (oneParticleEvent)
1136       fdNdEtaAnalysisOnePart->FillEvent(vtxMC[2], nAcceptedParticles);
1137
1138     if (eventTriggered)
1139     {
1140       fdNdEtaAnalysisTr->FillEvent(vtxMC[2], nAcceptedParticles);
1141       if (vtxESD)
1142         fdNdEtaAnalysisTrVtx->FillEvent(vtxMC[2], nAcceptedParticles);
1143     }
1144   }
1145 }
1146
1147 void AlidNdEtaTask::Terminate(Option_t *)
1148 {
1149   // The Terminate() function is the last function to be called during
1150   // a query. It always runs on the client, it can be used to present
1151   // the results graphically or save the results to file.
1152
1153   fOutput = dynamic_cast<TList*> (GetOutputData(1));
1154   if (!fOutput)
1155     Printf("ERROR: fOutput not available");
1156
1157   if (fOutput)
1158   {
1159     fdNdEtaAnalysisESD = dynamic_cast<dNdEtaAnalysis*> (fOutput->FindObject("fdNdEtaAnalysisESD"));
1160     fMult = dynamic_cast<TH1F*> (fOutput->FindObject("fMult"));
1161     fMultVtx = dynamic_cast<TH1F*> (fOutput->FindObject("fMultVtx"));
1162     for (Int_t i=0; i<3; ++i)
1163       fPartEta[i] = dynamic_cast<TH1F*> (fOutput->FindObject(Form("dndeta_check_%d", i)));
1164     fEvents = dynamic_cast<TH1F*> (fOutput->FindObject("dndeta_check_vertex"));
1165     fVertexResolution = dynamic_cast<TH2F*> (fOutput->FindObject("dndeta_vertex_resolution_z"));
1166
1167     fVertex = dynamic_cast<TH3F*> (fOutput->FindObject("vertex_check"));
1168     fVertexVsMult = dynamic_cast<TH3F*> (fOutput->FindObject("fVertexVsMult"));
1169     fPhi = dynamic_cast<TH1F*> (fOutput->FindObject("fPhi"));
1170     fRawPt = dynamic_cast<TH1F*> (fOutput->FindObject("fRawPt"));
1171     fEtaPhi = dynamic_cast<TH2F*> (fOutput->FindObject("fEtaPhi"));
1172     for (Int_t i=0; i<2; ++i)
1173       fZPhi[i] = dynamic_cast<TH2F*> (fOutput->FindObject(Form("fZPhi_%d", i)));
1174     fModuleMap = dynamic_cast<TH1F*> (fOutput->FindObject("fModuleMap"));
1175     fDeltaPhi = dynamic_cast<TH1F*> (fOutput->FindObject("fDeltaPhi"));
1176     fDeltaTheta = dynamic_cast<TH1F*> (fOutput->FindObject("fDeltaTheta"));
1177     fFiredChips = dynamic_cast<TH2F*> (fOutput->FindObject("fFiredChips"));
1178     fTrackletsVsClusters = dynamic_cast<TH2F*> (fOutput->FindObject("fTrackletsVsClusters"));
1179     fTrackletsVsUnassigned = dynamic_cast<TH2F*> (fOutput->FindObject("fTrackletsVsUnassigned"));
1180     fStats = dynamic_cast<TH1F*> (fOutput->FindObject("fStats"));
1181     fStats2 = dynamic_cast<TH2F*> (fOutput->FindObject("fStats2"));
1182
1183     fEsdTrackCuts = dynamic_cast<AliESDtrackCuts*> (fOutput->FindObject("fEsdTrackCuts"));
1184   }
1185
1186   if (!fdNdEtaAnalysisESD)
1187   {
1188     AliDebug(AliLog::kError, "ERROR: fdNdEtaAnalysisESD not available");
1189   }
1190   else
1191   {
1192     if (fMult && fMultVtx)
1193     {
1194       new TCanvas;
1195       fMult->Draw();
1196       fMultVtx->SetLineColor(2);
1197       fMultVtx->Draw("SAME");
1198     }
1199
1200     if (fFiredChips)
1201     {
1202       new TCanvas;
1203       fFiredChips->Draw("COLZ");
1204     }
1205
1206     if (fPartEta[0])
1207     {
1208       Int_t events1 = (Int_t) fEvents->Integral(fEvents->GetXaxis()->FindBin(-9.9), fEvents->GetXaxis()->FindBin(-0.001));
1209       Int_t events2 = (Int_t) fEvents->Integral(fEvents->GetXaxis()->FindBin(0.001), fEvents->GetXaxis()->FindBin(9.9));
1210
1211       Printf("%d events with vertex used", events1 + events2);
1212       Printf("%d total events with vertex, %d total triggered events, %d triggered events with 0 multiplicity", (Int_t) fMultVtx->Integral(), (Int_t) fMult->Integral(), (Int_t) fMult->GetBinContent(1));
1213
1214       if (events1 > 0 && events2 > 0)
1215       {
1216         fPartEta[0]->Scale(1.0 / (events1 + events2));
1217         fPartEta[1]->Scale(1.0 / events1);
1218         fPartEta[2]->Scale(1.0 / events2);
1219
1220         for (Int_t i=0; i<3; ++i)
1221           fPartEta[i]->Scale(1.0 / fPartEta[i]->GetBinWidth(1));
1222
1223         new TCanvas("control", "control", 500, 500);
1224         for (Int_t i=0; i<3; ++i)
1225         {
1226           fPartEta[i]->SetLineColor(i+1);
1227           fPartEta[i]->Draw((i==0) ? "" : "SAME");
1228         }
1229       }
1230     }
1231
1232     if (fEvents)
1233     {
1234         new TCanvas("control3", "control3", 500, 500);
1235         fEvents->Draw();
1236     }
1237
1238     if (fStats2)
1239     {
1240       new TCanvas;
1241       fStats2->Draw("TEXT");
1242     }
1243
1244     TFile* fout = new TFile("analysis_esd_raw.root", "RECREATE");
1245
1246     if (fdNdEtaAnalysisESD)
1247       fdNdEtaAnalysisESD->SaveHistograms();
1248
1249     if (fEsdTrackCuts)
1250       fEsdTrackCuts->SaveHistograms("esd_track_cuts");
1251
1252     if (fMult)
1253       fMult->Write();
1254
1255     if (fMultVtx)
1256       fMultVtx->Write();
1257
1258     for (Int_t i=0; i<3; ++i)
1259       if (fPartEta[i])
1260         fPartEta[i]->Write();
1261
1262     if (fEvents)
1263       fEvents->Write();
1264
1265     if (fVertexResolution)
1266     {
1267       fVertexResolution->Write();
1268       fVertexResolution->ProjectionX();
1269       fVertexResolution->ProjectionY();
1270     }
1271
1272     if (fDeltaPhi)
1273       fDeltaPhi->Write();
1274
1275     if (fDeltaTheta)
1276       fDeltaTheta->Write();
1277
1278     if (fPhi)
1279       fPhi->Write();
1280
1281     if (fRawPt)
1282       fRawPt->Write();
1283
1284     if (fEtaPhi)
1285       fEtaPhi->Write();
1286
1287     for (Int_t i=0; i<2; ++i)
1288       if (fZPhi[i])
1289         fZPhi[i]->Write();
1290
1291     if (fModuleMap)
1292       fModuleMap->Write();
1293
1294     if (fFiredChips)
1295       fFiredChips->Write();
1296
1297     if (fTrackletsVsClusters)
1298       fTrackletsVsClusters->Write();
1299
1300     if (fTrackletsVsUnassigned)
1301       fTrackletsVsUnassigned->Write();
1302
1303     if (fStats)
1304       fStats->Write();
1305
1306     if (fStats2)
1307       fStats2->Write();
1308
1309     if (fVertex)
1310       fVertex->Write();
1311
1312     if (fVertexVsMult)
1313       fVertexVsMult->Write();
1314
1315     fout->Write();
1316     fout->Close();
1317
1318     Printf("Writting result to analysis_esd_raw.root");
1319   }
1320
1321   if (fReadMC)
1322   {
1323     if (fOutput)
1324     {
1325       fdNdEtaAnalysis = dynamic_cast<dNdEtaAnalysis*> (fOutput->FindObject("dndeta"));
1326       fdNdEtaAnalysisND = dynamic_cast<dNdEtaAnalysis*> (fOutput->FindObject("dndetaND"));
1327       fdNdEtaAnalysisNSD = dynamic_cast<dNdEtaAnalysis*> (fOutput->FindObject("dndetaNSD"));
1328       fdNdEtaAnalysisOnePart = dynamic_cast<dNdEtaAnalysis*> (fOutput->FindObject("dndetaOnePart"));
1329       fdNdEtaAnalysisTr = dynamic_cast<dNdEtaAnalysis*> (fOutput->FindObject("dndetaTr"));
1330       fdNdEtaAnalysisTrVtx = dynamic_cast<dNdEtaAnalysis*> (fOutput->FindObject("dndetaTrVtx"));
1331       fdNdEtaAnalysisTracks = dynamic_cast<dNdEtaAnalysis*> (fOutput->FindObject("dndetaTracks"));
1332       fPartPt = dynamic_cast<TH1F*> (fOutput->FindObject("dndeta_check_pt"));
1333     }
1334
1335     if (!fdNdEtaAnalysis || !fdNdEtaAnalysisTr || !fdNdEtaAnalysisTrVtx || !fPartPt)
1336     {
1337       AliDebug(AliLog::kError, Form("ERROR: Histograms not available %p %p", (void*) fdNdEtaAnalysis, (void*) fPartPt));
1338       return;
1339     }
1340
1341     fdNdEtaAnalysis->Finish(0, -1, AlidNdEtaCorrection::kNone);
1342     fdNdEtaAnalysisND->Finish(0, -1, AlidNdEtaCorrection::kNone);
1343     fdNdEtaAnalysisNSD->Finish(0, -1, AlidNdEtaCorrection::kNone);
1344     fdNdEtaAnalysisOnePart->Finish(0, -1, AlidNdEtaCorrection::kNone);
1345     fdNdEtaAnalysisTr->Finish(0, -1, AlidNdEtaCorrection::kNone);
1346     fdNdEtaAnalysisTrVtx->Finish(0, -1, AlidNdEtaCorrection::kNone);
1347     fdNdEtaAnalysisTracks->Finish(0, -1, AlidNdEtaCorrection::kNone);
1348
1349     Int_t events = (Int_t) fdNdEtaAnalysis->GetData()->GetEventCorrection()->GetMeasuredHistogram()->Integral();
1350     fPartPt->Scale(1.0/events);
1351     fPartPt->Scale(1.0/fPartPt->GetBinWidth(1));
1352
1353     TFile* fout = new TFile("analysis_mc.root","RECREATE");
1354
1355     fdNdEtaAnalysis->SaveHistograms();
1356     fdNdEtaAnalysisND->SaveHistograms();
1357     fdNdEtaAnalysisNSD->SaveHistograms();
1358     fdNdEtaAnalysisOnePart->SaveHistograms();
1359     fdNdEtaAnalysisTr->SaveHistograms();
1360     fdNdEtaAnalysisTrVtx->SaveHistograms();
1361     fdNdEtaAnalysisTracks->SaveHistograms();
1362
1363     if (fPartPt)
1364       fPartPt->Write();
1365
1366     fout->Write();
1367     fout->Close();
1368
1369     Printf("Writting result to analysis_mc.root");
1370
1371     if (fPartPt)
1372     {
1373       new TCanvas("control2", "control2", 500, 500);
1374       fPartPt->Draw();
1375     }
1376   }
1377 }