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