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