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