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
  20 #include <AliLog.h>
  21 #include <AliESDVertex.h>
  22 #include <AliESDEvent.h>
  23 #include <AliStack.h>
  24 #include <AliHeader.h>
  25 #include <AliGenEventHeader.h>
  26 #include <AliMultiplicity.h>
  27 #include <AliAnalysisManager.h>
  28 #include <AliMCEventHandler.h>
  29 #include <AliMCEvent.h>
  30 #include <AliESDInputHandler.h>
  31
  32 #include "AliESDtrackCuts.h"
  33 #include "AliPWG0Helper.h"
  34 #include "AliCorrection.h"
  35 #include "AliCorrectionMatrix3D.h"
  36 #include "dNdEta/dNdEtaAnalysis.h"
  37
  38 ClassImp(AlidNdEtaTask)
  39
  40 AlidNdEtaTask::AlidNdEtaTask(const char* opt) :
  41   AliAnalysisTask("AlidNdEtaTask", ""),
  42   fESD(0),
  43   fOutput(0),
  44   fOption(opt),
  45   fAnalysisMode(AliPWG0Helper::kTPC),
  46   fTrigger(AliPWG0Helper::kMB1),
  47   fReadMC(kFALSE),
  48   fUseMCVertex(kFALSE),
  49   fUseMCKine(kFALSE),
  50   fEsdTrackCuts(0),
  51   fdNdEtaAnalysisESD(0),
  52   fMult(0),
  53   fMultVtx(0),
  54   fEvents(0),
  55   fVertexResolution(0),
  56   fdNdEtaAnalysis(0),
  57   fdNdEtaAnalysisNSD(0),
  58   fdNdEtaAnalysisTr(0),
  59   fdNdEtaAnalysisTrVtx(0),
  60   fdNdEtaAnalysisTracks(0),
  61   fVertex(0),
  62   fPartPt(0),
  63   fDeltaPhi(0)
  64 {
  65   //
  66   // Constructor. Initialization of pointers
  67   //
  68
  69   // Define input and output slots here
  70   DefineInput(0, TChain::Class());
  71   DefineOutput(0, TList::Class());
  72 }
  73
  74 AlidNdEtaTask::~AlidNdEtaTask()
  75 {
  76   //
  77   // Destructor
  78   //
  79
  80   // histograms are in the output list and deleted when the output
  81   // list is deleted by the TSelector dtor
  82
  83   if (fOutput) {
  84     delete fOutput;
  85     fOutput = 0;
  86   }
  87 }
  88
  89 //________________________________________________________________________
  90 void AlidNdEtaTask::ConnectInputData(Option_t *)
  91 {
  92   // Connect ESD
  93   // Called once
  94
  95   Printf("AlidNdEtaTask::ConnectInputData called");
  96
  97   TTree* tree = dynamic_cast<TTree*> (GetInputData(0));
  98   if (!tree) {
  99     Printf("ERROR: Could not read tree from input slot 0");
 100   } else {
 101     // Disable all branches and enable only the needed ones
 102     //tree->SetBranchStatus("*", 0);
 103
 104     tree->SetBranchStatus("TriggerMask", 1);
 105     tree->SetBranchStatus("SPDVertex*", 1);
 106     tree->SetBranchStatus("PrimaryVertex*", 1);
 107     tree->SetBranchStatus("TPCVertex*", 1);
 108
 109     if (fAnalysisMode == AliPWG0Helper::kSPD) {
 110       tree->SetBranchStatus("AliMultiplicity*", 1);
 111     }
 112
 113     if (fAnalysisMode == AliPWG0Helper::kTPC || fAnalysisMode == AliPWG0Helper::kTPCITS) {
 114       //AliESDtrackCuts::EnableNeededBranches(tree);
 115       tree->SetBranchStatus("Tracks*", 1);
 116     }
 117
 118     AliESDInputHandler *esdH = dynamic_cast<AliESDInputHandler*> (AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler());
 119
 120     if (!esdH) {
 121       Printf("ERROR: Could not get ESDInputHandler");
 122     } else
 123       fESD = esdH->GetEvent();
 124   }
 125
 126   // disable info messages of AliMCEvent (per event)
 127   AliLog::SetClassDebugLevel("AliMCEvent", AliLog::kWarning - AliLog::kDebug + 1);
 128 }
 129
 130 void AlidNdEtaTask::CreateOutputObjects()
 131 {
 132   // create result objects and add to output list
 133
 134   Printf("AlidNdEtaTask::CreateOutputObjects");
 135
 136   fOutput = new TList;
 137   fOutput->SetOwner();
 138
 139   fdNdEtaAnalysisESD = new dNdEtaAnalysis("fdNdEtaAnalysisESD", "fdNdEtaAnalysisESD", fAnalysisMode);
 140   fOutput->Add(fdNdEtaAnalysisESD);
 141
 142   fMult = new TH1F("fMult", "fMult;Ntracks;Count", 201, -0.5, 200.5);
 143   fOutput->Add(fMult);
 144
 145   fMultVtx = new TH1F("fMultVtx", "fMultVtx;Ntracks;Count", 201, -0.5, 200.5);
 146   fOutput->Add(fMultVtx);
 147
 148   for (Int_t i=0; i<3; ++i)
 149   {
 150     fPartEta[i] = new TH1F(Form("dndeta_check_%d", i), Form("dndeta_check_%d", i), 60, -6, 6);
 151     fPartEta[i]->Sumw2();
 152     fOutput->Add(fPartEta[i]);
 153   }
 154
 155   fEvents = new TH1F("dndeta_check_vertex", "dndeta_check_vertex", 160, -40, 40);
 156   fOutput->Add(fEvents);
 157
 158   fVertexResolution = new TH1F("dndeta_vertex_resolution_z", "dndeta_vertex_resolution_z", 1000, 0, 10);
 159   fOutput->Add(fVertexResolution);
 160
 161   if (fReadMC)
 162   {
 163     fdNdEtaAnalysis = new dNdEtaAnalysis("dndeta", "dndeta", fAnalysisMode);
 164     fOutput->Add(fdNdEtaAnalysis);
 165
 166     fdNdEtaAnalysisNSD = new dNdEtaAnalysis("dndetaNSD", "dndetaNSD", fAnalysisMode);
 167     fOutput->Add(fdNdEtaAnalysisNSD);
 168
 169     fdNdEtaAnalysisTr = new dNdEtaAnalysis("dndetaTr", "dndetaTr", fAnalysisMode);
 170     fOutput->Add(fdNdEtaAnalysisTr);
 171
 172     fdNdEtaAnalysisTrVtx = new dNdEtaAnalysis("dndetaTrVtx", "dndetaTrVtx", fAnalysisMode);
 173     fOutput->Add(fdNdEtaAnalysisTrVtx);
 174
 175     fdNdEtaAnalysisTracks = new dNdEtaAnalysis("dndetaTracks", "dndetaTracks", fAnalysisMode);
 176     fOutput->Add(fdNdEtaAnalysisTracks);
 177
 178     fVertex = new TH3F("vertex_check", "vertex_check", 50, -50, 50, 50, -50, 50, 50, -50, 50);
 179     fOutput->Add(fVertex);
 180
 181     fPartPt =  new TH1F("dndeta_check_pt", "dndeta_check_pt", 1000, 0, 10);
 182     fPartPt->Sumw2();
 183     fOutput->Add(fPartPt);
 184   }
 185
 186   if (fAnalysisMode == AliPWG0Helper::kSPD)
 187     fDeltaPhi = new TH1F("fDeltaPhi", "fDeltaPhi;#Delta #phi;Entries", 1000, -3.14, 3.14);
 188 }
 189
 190 void AlidNdEtaTask::Exec(Option_t*)
 191 {
 192   // process the event
 193
 194   // Check prerequisites
 195   if (!fESD)
 196   {
 197     AliDebug(AliLog::kError, "ESD branch not available");
 198     return;
 199   }
 200
 201   // trigger definition
 202   Bool_t eventTriggered = AliPWG0Helper::IsEventTriggered(fESD->GetTriggerMask(), fTrigger);
 203
 204   // get the ESD vertex
 205   const AliESDVertex* vtxESD = AliPWG0Helper::GetVertex(fESD, fAnalysisMode);
 206
 207   Double_t vtx[3];
 208   if (vtxESD) {
 209     vtxESD->GetXYZ(vtx);
 210   }
 211   else
 212     Printf("No vertex");
 213
 214   // fill z vertex resolution
 215   if (vtxESD)
 216     fVertexResolution->Fill(vtxESD->GetZRes());
 217
 218   // post the data already here
 219   PostData(0, fOutput);
 220
 221   // needed for syst. studies
 222   AliStack* stack = 0;
 223
 224   if (fUseMCVertex || fUseMCKine) {
 225     if (!fReadMC) {
 226       Printf("ERROR: fUseMCVertex or fUseMCKine set without fReadMC set!");
 227       return;
 228     }
 229
 230     AliMCEventHandler* eventHandler = dynamic_cast<AliMCEventHandler*> (AliAnalysisManager::GetAnalysisManager()->GetMCtruthEventHandler());
 231     if (!eventHandler) {
 232       Printf("ERROR: Could not retrieve MC event handler");
 233       return;
 234     }
 235
 236     AliMCEvent* mcEvent = eventHandler->MCEvent();
 237     if (!mcEvent) {
 238       Printf("ERROR: Could not retrieve MC event");
 239       return;
 240     }
 241
 242     AliHeader* header = mcEvent->Header();
 243     if (!header)
 244     {
 245       AliDebug(AliLog::kError, "Header not available");
 246       return;
 247     }
 248
 249     if (fUseMCVertex)
 250     {
 251       Printf("WARNING: Replacing vertex by MC vertex. This is for systematical checks only.");
 252       // get the MC vertex
 253       AliGenEventHeader* genHeader = header->GenEventHeader();
 254       TArrayF vtxMC(3);
 255       genHeader->PrimaryVertex(vtxMC);
 256
 257       vtx[2] = vtxMC[2];
 258     }
 259
 260     if (fUseMCKine)
 261     {
 262       stack = mcEvent->Stack();
 263       if (!stack)
 264       {
 265         AliDebug(AliLog::kError, "Stack not available");
 266         return;
 267       }
 268     }
 269   }
 270
 271   // create list of (label, eta, pt) tuples
 272   Int_t inputCount = 0;
 273   Int_t* labelArr = 0;
 274   Float_t* etaArr = 0;
 275   Float_t* ptArr = 0;
 276   if (fAnalysisMode == AliPWG0Helper::kSPD)
 277   {
 278     // get tracklets
 279     const AliMultiplicity* mult = fESD->GetMultiplicity();
 280     if (!mult)
 281     {
 282       AliDebug(AliLog::kError, "AliMultiplicity not available");
 283       return;
 284     }
 285
 286     labelArr = new Int_t[mult->GetNumberOfTracklets()];
 287     etaArr = new Float_t[mult->GetNumberOfTracklets()];
 288     ptArr = new Float_t[mult->GetNumberOfTracklets()];
 289
 290     if (fUseMCKine && stack)
 291       Printf("Processing only primaries (MC information used). This is for systematical checks only.");
 292
 293     // get multiplicity from ITS tracklets
 294     for (Int_t i=0; i<mult->GetNumberOfTracklets(); ++i)
 295     {
 296       //printf("%d %f %f %f\n", i, mult->GetTheta(i), mult->GetPhi(i), mult->GetDeltaPhi(i));
 297
 298       if (fUseMCKine && stack)
 299         if (mult->GetLabel(i, 0) < 0 || mult->GetLabel(i, 0) != mult->GetLabel(i, 1) || !stack->IsPhysicalPrimary(mult->GetLabel(i, 0)))
 300           continue;
 301
 302       Float_t deltaPhi = mult->GetDeltaPhi(i);
 303       // prevent values to be shifted by 2 Pi()
 304       if (deltaPhi < -TMath::Pi())
 305         deltaPhi += TMath::Pi() * 2;
 306       if (deltaPhi > TMath::Pi())
 307         deltaPhi -= TMath::Pi() * 2;
 308
 309       if (TMath::Abs(deltaPhi) > 1)
 310         printf("WARNING: Very high Delta Phi: %d %f %f %f\n", i, mult->GetTheta(i), mult->GetPhi(i), deltaPhi);
 311
 312       fDeltaPhi->Fill(deltaPhi);
 313
 314       etaArr[inputCount] = mult->GetEta(i);
 315       labelArr[inputCount] = mult->GetLabel(i, 0);
 316       ptArr[inputCount] = 0; // no pt for tracklets
 317       ++inputCount;
 318     }
 319
 320     Printf("Accepted %d tracklets", inputCount);
 321   }
 322   else if (fAnalysisMode == AliPWG0Helper::kTPC || fAnalysisMode == AliPWG0Helper::kTPCITS)
 323   {
 324     if (!fEsdTrackCuts)
 325     {
 326       AliDebug(AliLog::kError, "fESDTrackCuts not available");
 327       return;
 328     }
 329
 330     // get multiplicity from ESD tracks
 331     TObjArray* list = fEsdTrackCuts->GetAcceptedTracks(fESD);
 332     Int_t nGoodTracks = list->GetEntries();
 333
 334     labelArr = new Int_t[nGoodTracks];
 335     etaArr = new Float_t[nGoodTracks];
 336     ptArr = new Float_t[nGoodTracks];
 337
 338     // loop over esd tracks
 339     for (Int_t i=0; i<nGoodTracks; i++)
 340     {
 341       AliESDtrack* esdTrack = dynamic_cast<AliESDtrack*> (list->At(i));
 342       if (!esdTrack)
 343       {
 344         AliDebug(AliLog::kError, Form("ERROR: Could not retrieve track %d.", i));
 345         continue;
 346       }
 347
 348       etaArr[inputCount] = esdTrack->Eta();
 349       labelArr[inputCount] = TMath::Abs(esdTrack->GetLabel());
 350       ptArr[inputCount] = esdTrack->Pt();
 351       ++inputCount;
 352     }
 353
 354     Printf("Accepted %d tracks", nGoodTracks);
 355
 356     delete list;
 357   }
 358   else
 359     return;
 360
 361   // Processing of ESD information (always)
 362   if (eventTriggered)
 363   {
 364     // control hist
 365     fMult->Fill(inputCount);
 366     fdNdEtaAnalysisESD->FillTriggeredEvent(inputCount);
 367
 368     if (vtxESD)
 369     {
 370       // control hist
 371       fMultVtx->Fill(inputCount);
 372
 373       for (Int_t i=0; i<inputCount; ++i)
 374       {
 375         Float_t eta = etaArr[i];
 376         Float_t pt  = ptArr[i];
 377
 378         fdNdEtaAnalysisESD->FillTrack(vtx[2], eta, pt);
 379
 380         if (TMath::Abs(vtx[2]) < 20)
 381         {
 382           fPartEta[0]->Fill(eta);
 383
 384           if (vtx[2] < 0)
 385             fPartEta[1]->Fill(eta);
 386           else
 387             fPartEta[2]->Fill(eta);
 388         }
 389       }
 390
 391       // for event count per vertex
 392       fdNdEtaAnalysisESD->FillEvent(vtx[2], inputCount);
 393
 394       // control hist
 395       fEvents->Fill(vtx[2]);
 396     }
 397   }
 398
 399   if (fReadMC)   // Processing of MC information (optional)
 400   {
 401     AliMCEventHandler* eventHandler = dynamic_cast<AliMCEventHandler*> (AliAnalysisManager::GetAnalysisManager()->GetMCtruthEventHandler());
 402     if (!eventHandler) {
 403       Printf("ERROR: Could not retrieve MC event handler");
 404       return;
 405     }
 406
 407     AliMCEvent* mcEvent = eventHandler->MCEvent();
 408     if (!mcEvent) {
 409       Printf("ERROR: Could not retrieve MC event");
 410       return;
 411     }
 412
 413     stack = mcEvent->Stack();
 414     if (!stack)
 415     {
 416       AliDebug(AliLog::kError, "Stack not available");
 417       return;
 418     }
 419
 420     AliHeader* header = mcEvent->Header();
 421     if (!header)
 422     {
 423       AliDebug(AliLog::kError, "Header not available");
 424       return;
 425     }
 426
 427     // get the MC vertex
 428     AliGenEventHeader* genHeader = header->GenEventHeader();
 429     TArrayF vtxMC(3);
 430     genHeader->PrimaryVertex(vtxMC);
 431
 432     // get process type; NB: this only works for Pythia
 433     Int_t processType = AliPWG0Helper::GetPythiaEventProcessType(header);
 434     AliDebug(AliLog::kDebug+1, Form("Found pythia process type %d", processType));
 435
 436     if (processType<0)
 437       AliDebug(AliLog::kError, Form("Unknown Pythia process type %d.", processType));
 438
 439     // loop over mc particles
 440     Int_t nPrim  = stack->GetNprimary();
 441
 442     Int_t nAcceptedParticles = 0;
 443
 444     for (Int_t iMc = 0; iMc < nPrim; ++iMc)
 445     {
 446       //Printf("Getting particle %d", iMc);
 447       TParticle* particle = stack->Particle(iMc);
 448
 449       if (!particle)
 450         continue;
 451
 452       if (AliPWG0Helper::IsPrimaryCharged(particle, nPrim) == kFALSE)
 453         continue;
 454
 455       AliDebug(AliLog::kDebug+1, Form("Accepted primary %d, unique ID: %d", iMc, particle->GetUniqueID()));
 456       Float_t eta = particle->Eta();
 457       Float_t pt = particle->Pt();
 458
 459        // 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))
 460       if (TMath::Abs(eta) < 1.5 && pt > 0.3)
 461         nAcceptedParticles++;
 462
 463       fdNdEtaAnalysis->FillTrack(vtxMC[2], eta, pt);
 464       fVertex->Fill(particle->Vx(), particle->Vy(), particle->Vz());
 465
 466       if (processType != 92 && processType != 93)
 467         fdNdEtaAnalysisNSD->FillTrack(vtxMC[2], eta, pt);
 468
 469       if (eventTriggered)
 470       {
 471         fdNdEtaAnalysisTr->FillTrack(vtxMC[2], eta, pt);
 472         if (vtxESD)
 473           fdNdEtaAnalysisTrVtx->FillTrack(vtxMC[2], eta, pt);
 474       }
 475
 476       if (TMath::Abs(eta) < 0.8)
 477         fPartPt->Fill(pt);
 478     }
 479
 480     fdNdEtaAnalysis->FillEvent(vtxMC[2], nAcceptedParticles);
 481     if (processType != 92 && processType != 93)
 482       fdNdEtaAnalysisNSD->FillEvent(vtxMC[2], nAcceptedParticles);
 483
 484     if (eventTriggered)
 485     {
 486       fdNdEtaAnalysisTr->FillEvent(vtxMC[2], nAcceptedParticles);
 487       if (vtxESD)
 488         fdNdEtaAnalysisTrVtx->FillEvent(vtxMC[2], nAcceptedParticles);
 489     }
 490
 491     if (eventTriggered && vtxESD)
 492     {
 493       // from tracks is only done for triggered and vertex reconstructed events
 494
 495       for (Int_t i=0; i<inputCount; ++i)
 496       {
 497         Int_t label = labelArr[i];
 498
 499         if (label < 0)
 500           continue;
 501
 502         //Printf("Getting particle of track %d", label);
 503         TParticle* particle = stack->Particle(label);
 504
 505         if (!particle)
 506         {
 507           AliDebug(AliLog::kError, Form("ERROR: Could not retrieve particle %d.", label));
 508           continue;
 509         }
 510
 511         fdNdEtaAnalysisTracks->FillTrack(vtxMC[2], particle->Eta(), particle->Pt());
 512       } // end of track loop
 513
 514       // for event count per vertex
 515       fdNdEtaAnalysisTracks->FillEvent(vtxMC[2], inputCount);
 516     }
 517   }
 518
 519   delete[] etaArr;
 520   delete[] labelArr;
 521   delete[] ptArr;
 522 }
 523
 524 void AlidNdEtaTask::Terminate(Option_t *)
 525 {
 526   // The Terminate() function is the last function to be called during
 527   // a query. It always runs on the client, it can be used to present
 528   // the results graphically or save the results to file.
 529
 530   fOutput = dynamic_cast<TList*> (GetOutputData(0));
 531   if (!fOutput) {
 532     Printf("ERROR: fOutput not available");
 533     return;
 534   }
 535
 536   fdNdEtaAnalysisESD = dynamic_cast<dNdEtaAnalysis*> (fOutput->FindObject("fdNdEtaAnalysisESD"));
 537   fMult = dynamic_cast<TH1F*> (fOutput->FindObject("fMult"));
 538   fMultVtx = dynamic_cast<TH1F*> (fOutput->FindObject("fMultVtx"));
 539   fEvents = dynamic_cast<TH1F*> (fOutput->FindObject("dndeta_check_vertex"));
 540   fVertexResolution = dynamic_cast<TH1F*> (fOutput->FindObject("dndeta_vertex_resolution_z"));
 541   for (Int_t i=0; i<3; ++i)
 542     fPartEta[i] = dynamic_cast<TH1F*> (fOutput->FindObject(Form("dndeta_check_%d", i)));
 543
 544   if (!fdNdEtaAnalysisESD)
 545   {
 546     AliDebug(AliLog::kError, "ERROR: fdNdEtaAnalysisESD not available");
 547     return;
 548   }
 549
 550   if (fMult && fMultVtx)
 551   {
 552     new TCanvas;
 553     fMult->Draw();
 554     fMultVtx->SetLineColor(2);
 555     fMultVtx->Draw("SAME");
 556   }
 557
 558   if (fPartEta[0])
 559   {
 560     Int_t events1 = (Int_t) fEvents->Integral(fEvents->GetXaxis()->FindBin(-19.9), fEvents->GetXaxis()->FindBin(-0.001));
 561     Int_t events2 = (Int_t) fEvents->Integral(fEvents->GetXaxis()->FindBin(0.001), fEvents->GetXaxis()->FindBin(19.9));
 562
 563     Printf("%d events with vertex used", events1 + events2);
 564
 565     if (events1 > 0 && events2 > 0)
 566     {
 567     fPartEta[0]->Scale(1.0 / (events1 + events2));
 568     fPartEta[1]->Scale(1.0 / events1);
 569     fPartEta[2]->Scale(1.0 / events2);
 570
 571     for (Int_t i=0; i<3; ++i)
 572       fPartEta[i]->Scale(1.0 / fPartEta[i]->GetBinWidth(1));
 573
 574     new TCanvas("control", "control", 500, 500);
 575     for (Int_t i=0; i<3; ++i)
 576     {
 577       fPartEta[i]->SetLineColor(i+1);
 578       fPartEta[i]->Draw((i==0) ? "" : "SAME");
 579     }
 580     }
 581   }
 582
 583   if (fEvents)
 584   {
 585     new TCanvas("control3", "control3", 500, 500);
 586     fEvents->Draw();
 587   }
 588
 589   TFile* fout = new TFile("analysis_esd_raw.root", "RECREATE");
 590
 591   if (fdNdEtaAnalysisESD)
 592     fdNdEtaAnalysisESD->SaveHistograms();
 593
 594   if (fEsdTrackCuts)
 595     fEsdTrackCuts->SaveHistograms("esd_tracks_cuts");
 596
 597   if (fMult)
 598     fMult->Write();
 599
 600   if (fMultVtx)
 601     fMultVtx->Write();
 602
 603   for (Int_t i=0; i<3; ++i)
 604     if (fPartEta[i])
 605       fPartEta[i]->Write();
 606
 607   if (fEvents)
 608     fEvents->Write();
 609
 610   if (fVertexResolution)
 611     fVertexResolution->Write();
 612
 613   if (fDeltaPhi)
 614     fDeltaPhi->Write();
 615
 616   fout->Write();
 617   fout->Close();
 618
 619   Printf("Writting result to analysis_esd_raw.root");
 620
 621   if (fReadMC)
 622   {
 623     fdNdEtaAnalysis = dynamic_cast<dNdEtaAnalysis*> (fOutput->FindObject("dndeta"));
 624     fdNdEtaAnalysisNSD = dynamic_cast<dNdEtaAnalysis*> (fOutput->FindObject("dndetaNSD"));
 625     fdNdEtaAnalysisTr = dynamic_cast<dNdEtaAnalysis*> (fOutput->FindObject("dndetaTr"));
 626     fdNdEtaAnalysisTrVtx = dynamic_cast<dNdEtaAnalysis*> (fOutput->FindObject("dndetaTrVtx"));
 627     fdNdEtaAnalysisTracks = dynamic_cast<dNdEtaAnalysis*> (fOutput->FindObject("dndetaTracks"));
 628     fPartPt = dynamic_cast<TH1F*> (fOutput->FindObject("dndeta_check_pt"));
 629
 630     if (!fdNdEtaAnalysis || !fdNdEtaAnalysisTr || !fdNdEtaAnalysisTrVtx || !fPartPt)
 631     {
 632       AliDebug(AliLog::kError, Form("ERROR: Histograms not available %p %p", (void*) fdNdEtaAnalysis, (void*) fPartPt));
 633       return;
 634     }
 635
 636     fdNdEtaAnalysis->Finish(0, -1, AlidNdEtaCorrection::kNone);
 637     fdNdEtaAnalysisNSD->Finish(0, -1, AlidNdEtaCorrection::kNone);
 638     fdNdEtaAnalysisTr->Finish(0, -1, AlidNdEtaCorrection::kNone);
 639     fdNdEtaAnalysisTrVtx->Finish(0, -1, AlidNdEtaCorrection::kNone);
 640     fdNdEtaAnalysisTracks->Finish(0, -1, AlidNdEtaCorrection::kNone);
 641
 642     Int_t events = (Int_t) fdNdEtaAnalysis->GetData()->GetEventCorrection()->GetMeasuredHistogram()->Integral();
 643     fPartPt->Scale(1.0/events);
 644     fPartPt->Scale(1.0/fPartPt->GetBinWidth(1));
 645
 646     fout = new TFile("analysis_mc.root","RECREATE");
 647
 648     fdNdEtaAnalysis->SaveHistograms();
 649     fdNdEtaAnalysisNSD->SaveHistograms();
 650     fdNdEtaAnalysisTr->SaveHistograms();
 651     fdNdEtaAnalysisTrVtx->SaveHistograms();
 652     fdNdEtaAnalysisTracks->SaveHistograms();
 653     fPartPt->Write();
 654
 655     fout->Write();
 656     fout->Close();
 657
 658     Printf("Writting result to analysis_mc.root");
 659
 660     if (fPartPt)
 661     {
 662       new TCanvas("control2", "control2", 500, 500);
 663       fPartPt->Draw();
 664     }
 665   }
 666 }