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