[u/mrichter/AliRoot.git] / PWG0 / dNdEta / AlidNdEtaTask.cxx

/* $Id$ */

#include "AlidNdEtaTask.h"

#include <TStyle.h>
#include <TSystem.h>
#include <TCanvas.h>
#include <TVector3.h>
#include <TChain.h>
#include <TFile.h>
#include <TH1F.h>
#include <TH2F.h>
#include <TH3F.h>
#include <TParticle.h>
#include <TRandom.h>
#include <TNtuple.h>
#include <TObjString.h>
#include <TF1.h>

#include <AliLog.h>
#include <AliESDVertex.h>
#include <AliESDEvent.h>
#include <AliStack.h>
#include <AliHeader.h>
#include <AliGenEventHeader.h>
#include <AliMultiplicity.h>
#include <AliAnalysisManager.h>
#include <AliMCEventHandler.h>
#include <AliMCEvent.h>
#include <AliESDInputHandler.h>

#include "AliESDtrackCuts.h"
#include "AliPWG0Helper.h"
#include "AliCorrection.h"
#include "AliCorrectionMatrix3D.h"
#include "dNdEta/dNdEtaAnalysis.h"

ClassImp(AlidNdEtaTask)

AlidNdEtaTask::AlidNdEtaTask(const char* opt) :
  AliAnalysisTask("AlidNdEtaTask", ""),
  fESD(0),
  fOutput(0),
  fOption(opt),
  fAnalysisMode(AliPWG0Helper::kTPC),
  fTrigger(AliPWG0Helper::kMB1),
  fReadMC(kFALSE),
  fUseMCVertex(kFALSE),
  fUseMCKine(kFALSE),
  fEsdTrackCuts(0),
  fdNdEtaAnalysisESD(0),
  fMult(0),
  fMultVtx(0),
  fEvents(0),
  fVertexResolution(0),
  fdNdEtaAnalysis(0),
  fdNdEtaAnalysisNSD(0),
  fdNdEtaAnalysisTr(0),
  fdNdEtaAnalysisTrVtx(0),
  fdNdEtaAnalysisTracks(0),
  fVertex(0),
  fPartPt(0),
  fDeltaPhi(0)
{
  //
  // Constructor. Initialization of pointers
  //

  // Define input and output slots here
  DefineInput(0, TChain::Class());
  DefineOutput(0, TList::Class());

  AliLog::SetClassDebugLevel("AlidNdEtaTask", AliLog::kWarning);
}

AlidNdEtaTask::~AlidNdEtaTask()
{
  //
  // Destructor
  //

  // histograms are in the output list and deleted when the output
  // list is deleted by the TSelector dtor

  if (fOutput) {
    delete fOutput;
    fOutput = 0;
  }
}

Bool_t AlidNdEtaTask::Notify()
{
  static Int_t count = 0;
  count++;
  Printf("Processing %d. file", count);
  return kTRUE;
}

//________________________________________________________________________
void AlidNdEtaTask::ConnectInputData(Option_t *)
{
  // Connect ESD
  // Called once

  Printf("AlidNdEtaTask::ConnectInputData called");

  AliESDInputHandler *esdH = dynamic_cast<AliESDInputHandler*> (AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler());

  if (!esdH) {
    Printf("ERROR: Could not get ESDInputHandler");
  } else {
    fESD = esdH->GetEvent();

    TTree* tree = esdH->GetTree();
    if (!tree) {
      Printf("ERROR: Could not read tree");
    } else {
      // Disable all branches and enable only the needed ones
      tree->SetBranchStatus("*", 0);

      tree->SetBranchStatus("AliESDHeader*", 1);
      tree->SetBranchStatus("*Vertex*", 1);

      if (fAnalysisMode == AliPWG0Helper::kSPD) {
        tree->SetBranchStatus("AliMultiplicity*", 1);
      }

      if (fAnalysisMode == AliPWG0Helper::kTPC || fAnalysisMode == AliPWG0Helper::kTPCITS) {
        //AliESDtrackCuts::EnableNeededBranches(tree);
        tree->SetBranchStatus("Tracks*", 1);
      }

    }
  }

  // disable info messages of AliMCEvent (per event)
  AliLog::SetClassDebugLevel("AliMCEvent", AliLog::kWarning - AliLog::kDebug + 1);
}

void AlidNdEtaTask::CreateOutputObjects()
{
  // create result objects and add to output list

  Printf("AlidNdEtaTask::CreateOutputObjects");

  fOutput = new TList;
  fOutput->SetOwner();

  fdNdEtaAnalysisESD = new dNdEtaAnalysis("fdNdEtaAnalysisESD", "fdNdEtaAnalysisESD", fAnalysisMode);
  fOutput->Add(fdNdEtaAnalysisESD);

  fMult = new TH1F("fMult", "fMult;Ntracks;Count", 201, -0.5, 200.5);
  fOutput->Add(fMult);

  fMultVtx = new TH1F("fMultVtx", "fMultVtx;Ntracks;Count", 201, -0.5, 200.5);
  fOutput->Add(fMultVtx);

  for (Int_t i=0; i<3; ++i)
  {
    fPartEta[i] = new TH1F(Form("dndeta_check_%d", i), Form("dndeta_check_%d", i), 60, -6, 6);
    fPartEta[i]->Sumw2();
    fOutput->Add(fPartEta[i]);
  }

  fEvents = new TH1F("dndeta_check_vertex", "dndeta_check_vertex", 160, -40, 40);
  fOutput->Add(fEvents);

  fVertexResolution = new TH1F("dndeta_vertex_resolution_z", "dndeta_vertex_resolution_z", 1000, 0, 10);
  fOutput->Add(fVertexResolution);

  if (fReadMC)
  {
    fdNdEtaAnalysis = new dNdEtaAnalysis("dndeta", "dndeta", fAnalysisMode);
    fOutput->Add(fdNdEtaAnalysis);

    fdNdEtaAnalysisNSD = new dNdEtaAnalysis("dndetaNSD", "dndetaNSD", fAnalysisMode);
    fOutput->Add(fdNdEtaAnalysisNSD);

    fdNdEtaAnalysisTr = new dNdEtaAnalysis("dndetaTr", "dndetaTr", fAnalysisMode);
    fOutput->Add(fdNdEtaAnalysisTr);

    fdNdEtaAnalysisTrVtx = new dNdEtaAnalysis("dndetaTrVtx", "dndetaTrVtx", fAnalysisMode);
    fOutput->Add(fdNdEtaAnalysisTrVtx);

    fdNdEtaAnalysisTracks = new dNdEtaAnalysis("dndetaTracks", "dndetaTracks", fAnalysisMode);
    fOutput->Add(fdNdEtaAnalysisTracks);

    fVertex = new TH3F("vertex_check", "vertex_check", 50, -50, 50, 50, -50, 50, 50, -50, 50);
    fOutput->Add(fVertex);

    fPartPt =  new TH1F("dndeta_check_pt", "dndeta_check_pt", 1000, 0, 10);
    fPartPt->Sumw2();
    fOutput->Add(fPartPt);
  }

  if (fAnalysisMode == AliPWG0Helper::kSPD)
    fDeltaPhi = new TH1F("fDeltaPhi", "fDeltaPhi;#Delta #phi;Entries", 1000, -3.14, 3.14);
}

void AlidNdEtaTask::Exec(Option_t*)
{
  // process the event

  // Check prerequisites
  if (!fESD)
  {
    AliDebug(AliLog::kError, "ESD branch not available");
    return;
  }

  // trigger definition
  Bool_t eventTriggered = AliPWG0Helper::IsEventTriggered(fESD->GetTriggerMask(), fTrigger);

  // get the ESD vertex
  const AliESDVertex* vtxESD = AliPWG0Helper::GetVertex(fESD, fAnalysisMode);
  
  Double_t vtx[3];
  if (vtxESD) {
    vtxESD->GetXYZ(vtx);
  }
  else
    Printf("No vertex");
  
  // fill z vertex resolution
  if (vtxESD)
    fVertexResolution->Fill(vtxESD->GetZRes());

  // post the data already here
  PostData(0, fOutput);

  // needed for syst. studies
  AliStack* stack = 0;
  
  if (fUseMCVertex || fUseMCKine) {
    if (!fReadMC) {
      Printf("ERROR: fUseMCVertex or fUseMCKine set without fReadMC set!");
      return;
    }

    AliMCEventHandler* eventHandler = dynamic_cast<AliMCEventHandler*> (AliAnalysisManager::GetAnalysisManager()->GetMCtruthEventHandler());
    if (!eventHandler) {
      Printf("ERROR: Could not retrieve MC event handler");
      return;
    }

    AliMCEvent* mcEvent = eventHandler->MCEvent();
    if (!mcEvent) {
      Printf("ERROR: Could not retrieve MC event");
      return;
    }

    AliHeader* header = mcEvent->Header();
    if (!header)
    {
      AliDebug(AliLog::kError, "Header not available");
      return;
    }

    if (fUseMCVertex)
    {
      Printf("WARNING: Replacing vertex by MC vertex. This is for systematical checks only.");
      // get the MC vertex
      AliGenEventHeader* genHeader = header->GenEventHeader();
      TArrayF vtxMC(3);
      genHeader->PrimaryVertex(vtxMC);

      vtx[2] = vtxMC[2];
    }

    if (fUseMCKine)
    {
      stack = mcEvent->Stack();
      if (!stack)
      {
        AliDebug(AliLog::kError, "Stack not available");
        return;
      }
    }
  }

  // create list of (label, eta, pt) tuples
  Int_t inputCount = 0;
  Int_t* labelArr = 0;
  Float_t* etaArr = 0;
  Float_t* ptArr = 0;
  if (fAnalysisMode == AliPWG0Helper::kSPD)
  {
    // get tracklets
    const AliMultiplicity* mult = fESD->GetMultiplicity();
    if (!mult)
    {
      AliDebug(AliLog::kError, "AliMultiplicity not available");
      return;
    }

    labelArr = new Int_t[mult->GetNumberOfTracklets()];
    etaArr = new Float_t[mult->GetNumberOfTracklets()];
    ptArr = new Float_t[mult->GetNumberOfTracklets()];

    if (fUseMCKine && stack)
      Printf("Processing only primaries (MC information used). This is for systematical checks only.");

    // get multiplicity from ITS tracklets
    for (Int_t i=0; i<mult->GetNumberOfTracklets(); ++i)
    {
      //printf("%d %f %f %f\n", i, mult->GetTheta(i), mult->GetPhi(i), mult->GetDeltaPhi(i));

      if (fUseMCKine && stack)
        if (mult->GetLabel(i, 0) < 0 || mult->GetLabel(i, 0) != mult->GetLabel(i, 1) || !stack->IsPhysicalPrimary(mult->GetLabel(i, 0)))
          continue;
      
      Float_t deltaPhi = mult->GetDeltaPhi(i);
      // prevent values to be shifted by 2 Pi()
      if (deltaPhi < -TMath::Pi())
        deltaPhi += TMath::Pi() * 2;
      if (deltaPhi > TMath::Pi())
        deltaPhi -= TMath::Pi() * 2;

      if (TMath::Abs(deltaPhi) > 1)
        printf("WARNING: Very high Delta Phi: %d %f %f %f\n", i, mult->GetTheta(i), mult->GetPhi(i), deltaPhi);

      fDeltaPhi->Fill(deltaPhi);

      etaArr[inputCount] = mult->GetEta(i);
      labelArr[inputCount] = mult->GetLabel(i, 0);
      ptArr[inputCount] = 0; // no pt for tracklets
      ++inputCount;
    }

    //Printf("Accepted %d tracklets", inputCount);
  }
  else if (fAnalysisMode == AliPWG0Helper::kTPC || fAnalysisMode == AliPWG0Helper::kTPCITS)
  {
    if (!fEsdTrackCuts)
    {
      AliDebug(AliLog::kError, "fESDTrackCuts not available");
      return;
    }

    // get multiplicity from ESD tracks
    TObjArray* list = fEsdTrackCuts->GetAcceptedTracks(fESD);
    Int_t nGoodTracks = list->GetEntries();

    labelArr = new Int_t[nGoodTracks];
    etaArr = new Float_t[nGoodTracks];
    ptArr = new Float_t[nGoodTracks];

    // loop over esd tracks
    for (Int_t i=0; i<nGoodTracks; i++)
    {
      AliESDtrack* esdTrack = dynamic_cast<AliESDtrack*> (list->At(i));
      if (!esdTrack)
      {
        AliDebug(AliLog::kError, Form("ERROR: Could not retrieve track %d.", i));
        continue;
      }

      etaArr[inputCount] = esdTrack->Eta();
      labelArr[inputCount] = TMath::Abs(esdTrack->GetLabel());
      ptArr[inputCount] = esdTrack->Pt();
      ++inputCount;
    }

    Printf("Accepted %d tracks", nGoodTracks);

    delete list;
  }
  else
    return;

  // Processing of ESD information (always)
  if (eventTriggered)
  {
    // control hist
    fMult->Fill(inputCount);
    fdNdEtaAnalysisESD->FillTriggeredEvent(inputCount);

    if (vtxESD)
    {
      // control hist
      fMultVtx->Fill(inputCount);

      for (Int_t i=0; i<inputCount; ++i)
      {
        Float_t eta = etaArr[i];
        Float_t pt  = ptArr[i];

        fdNdEtaAnalysisESD->FillTrack(vtx[2], eta, pt);

        if (TMath::Abs(vtx[2]) < 20)
        {
          fPartEta[0]->Fill(eta);

          if (vtx[2] < 0)
            fPartEta[1]->Fill(eta);
          else
            fPartEta[2]->Fill(eta);
        }
      }

      // for event count per vertex
      fdNdEtaAnalysisESD->FillEvent(vtx[2], inputCount);

      // control hist
      fEvents->Fill(vtx[2]);
    }
  }

  if (fReadMC)   // Processing of MC information (optional)
  {
    AliMCEventHandler* eventHandler = dynamic_cast<AliMCEventHandler*> (AliAnalysisManager::GetAnalysisManager()->GetMCtruthEventHandler());
    if (!eventHandler) {
      Printf("ERROR: Could not retrieve MC event handler");
      return;
    }

    AliMCEvent* mcEvent = eventHandler->MCEvent();
    if (!mcEvent) {
      Printf("ERROR: Could not retrieve MC event");
      return;
    }

    stack = mcEvent->Stack();
    if (!stack)
    {
      AliDebug(AliLog::kError, "Stack not available");
      return;
    }

    AliHeader* header = mcEvent->Header();
    if (!header)
    {
      AliDebug(AliLog::kError, "Header not available");
      return;
    }

    // get the MC vertex
    AliGenEventHeader* genHeader = header->GenEventHeader();
    if (!genHeader)
    {
      AliDebug(AliLog::kError, "Could not retrieve genHeader from Header");
      return;
    }

    TArrayF vtxMC(3);
    genHeader->PrimaryVertex(vtxMC);

    // get process type; NB: this only works for Pythia
    Int_t processType = AliPWG0Helper::GetPythiaEventProcessType(header);
    AliDebug(AliLog::kDebug+1, Form("Found pythia process type %d", processType));

    if (processType<0)
      AliDebug(AliLog::kError, Form("Unknown Pythia process type %d.", processType));

    // loop over mc particles
    Int_t nPrim  = stack->GetNprimary();

    Int_t nAcceptedParticles = 0;

    for (Int_t iMc = 0; iMc < nPrim; ++iMc)
    {
      //Printf("Getting particle %d", iMc);
      TParticle* particle = stack->Particle(iMc);

      if (!particle)
        continue;

      if (AliPWG0Helper::IsPrimaryCharged(particle, nPrim) == kFALSE)
        continue;

      AliDebug(AliLog::kDebug+1, Form("Accepted primary %d, unique ID: %d", iMc, particle->GetUniqueID()));
      Float_t eta = particle->Eta();
      Float_t pt = particle->Pt();

       // 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))
      if (TMath::Abs(eta) < 1.5 && pt > 0.3)
        nAcceptedParticles++;

      fdNdEtaAnalysis->FillTrack(vtxMC[2], eta, pt);
      fVertex->Fill(particle->Vx(), particle->Vy(), particle->Vz());

      if (processType != 92 && processType != 93)
        fdNdEtaAnalysisNSD->FillTrack(vtxMC[2], eta, pt);

      if (eventTriggered)
      {
        fdNdEtaAnalysisTr->FillTrack(vtxMC[2], eta, pt);
        if (vtxESD)
          fdNdEtaAnalysisTrVtx->FillTrack(vtxMC[2], eta, pt);
      }

      if (TMath::Abs(eta) < 0.8)
        fPartPt->Fill(pt);
    }

    fdNdEtaAnalysis->FillEvent(vtxMC[2], nAcceptedParticles);
    if (processType != 92 && processType != 93)
      fdNdEtaAnalysisNSD->FillEvent(vtxMC[2], nAcceptedParticles);

    if (eventTriggered)
    {
      fdNdEtaAnalysisTr->FillEvent(vtxMC[2], nAcceptedParticles);
      if (vtxESD)
        fdNdEtaAnalysisTrVtx->FillEvent(vtxMC[2], nAcceptedParticles);
    }

    if (eventTriggered && vtxESD)
    {
      // from tracks is only done for triggered and vertex reconstructed events

      for (Int_t i=0; i<inputCount; ++i)
      {
        Int_t label = labelArr[i];

        if (label < 0)
          continue;

        //Printf("Getting particle of track %d", label);
        TParticle* particle = stack->Particle(label);

        if (!particle)
        {
          AliDebug(AliLog::kError, Form("ERROR: Could not retrieve particle %d.", label));
          continue;
        }

        fdNdEtaAnalysisTracks->FillTrack(vtxMC[2], particle->Eta(), particle->Pt());
      } // end of track loop

      // for event count per vertex
      fdNdEtaAnalysisTracks->FillEvent(vtxMC[2], inputCount);
    }
  }

  delete[] etaArr;
  delete[] labelArr;
  delete[] ptArr;
}

void AlidNdEtaTask::Terminate(Option_t *)
{
  // The Terminate() function is the last function to be called during
  // a query. It always runs on the client, it can be used to present
  // the results graphically or save the results to file.

  fOutput = dynamic_cast<TList*> (GetOutputData(0));
  if (!fOutput) {
    Printf("ERROR: fOutput not available");
    return;
  }

  fdNdEtaAnalysisESD = dynamic_cast<dNdEtaAnalysis*> (fOutput->FindObject("fdNdEtaAnalysisESD"));
  fMult = dynamic_cast<TH1F*> (fOutput->FindObject("fMult"));
  fMultVtx = dynamic_cast<TH1F*> (fOutput->FindObject("fMultVtx"));
  fEvents = dynamic_cast<TH1F*> (fOutput->FindObject("dndeta_check_vertex"));
  fVertexResolution = dynamic_cast<TH1F*> (fOutput->FindObject("dndeta_vertex_resolution_z"));
  for (Int_t i=0; i<3; ++i)
    fPartEta[i] = dynamic_cast<TH1F*> (fOutput->FindObject(Form("dndeta_check_%d", i)));

  if (!fdNdEtaAnalysisESD)
  {
    AliDebug(AliLog::kError, "ERROR: fdNdEtaAnalysisESD not available");
    return;
  }

  if (fMult && fMultVtx)
  {
    new TCanvas;
    fMult->Draw();
    fMultVtx->SetLineColor(2);
    fMultVtx->Draw("SAME");
  }

  if (fPartEta[0])
  {
    Int_t events1 = (Int_t) fEvents->Integral(fEvents->GetXaxis()->FindBin(-19.9), fEvents->GetXaxis()->FindBin(-0.001));
    Int_t events2 = (Int_t) fEvents->Integral(fEvents->GetXaxis()->FindBin(0.001), fEvents->GetXaxis()->FindBin(19.9));

    Printf("%d events with vertex used", events1 + events2);

    if (events1 > 0 && events2 > 0)
    {
    fPartEta[0]->Scale(1.0 / (events1 + events2));
    fPartEta[1]->Scale(1.0 / events1);
    fPartEta[2]->Scale(1.0 / events2);

    for (Int_t i=0; i<3; ++i)
      fPartEta[i]->Scale(1.0 / fPartEta[i]->GetBinWidth(1));

    new TCanvas("control", "control", 500, 500);
    for (Int_t i=0; i<3; ++i)
    {
      fPartEta[i]->SetLineColor(i+1);
      fPartEta[i]->Draw((i==0) ? "" : "SAME");
    }
    }
  }

  if (fEvents)
  {
    new TCanvas("control3", "control3", 500, 500);
    fEvents->Draw();
  }

  TFile* fout = new TFile("analysis_esd_raw.root", "RECREATE");

  if (fdNdEtaAnalysisESD)
    fdNdEtaAnalysisESD->SaveHistograms();

  if (fEsdTrackCuts)
    fEsdTrackCuts->SaveHistograms("esd_tracks_cuts");

  if (fMult)
    fMult->Write();

  if (fMultVtx)
    fMultVtx->Write();

  for (Int_t i=0; i<3; ++i)
    if (fPartEta[i])
      fPartEta[i]->Write();

  if (fEvents)
    fEvents->Write();

  if (fVertexResolution)
    fVertexResolution->Write();

  if (fDeltaPhi)
    fDeltaPhi->Write();

  fout->Write();
  fout->Close();

  Printf("Writting result to analysis_esd_raw.root");

  if (fReadMC)
  {
    fdNdEtaAnalysis = dynamic_cast<dNdEtaAnalysis*> (fOutput->FindObject("dndeta"));
    fdNdEtaAnalysisNSD = dynamic_cast<dNdEtaAnalysis*> (fOutput->FindObject("dndetaNSD"));
    fdNdEtaAnalysisTr = dynamic_cast<dNdEtaAnalysis*> (fOutput->FindObject("dndetaTr"));
    fdNdEtaAnalysisTrVtx = dynamic_cast<dNdEtaAnalysis*> (fOutput->FindObject("dndetaTrVtx"));
    fdNdEtaAnalysisTracks = dynamic_cast<dNdEtaAnalysis*> (fOutput->FindObject("dndetaTracks"));
    fPartPt = dynamic_cast<TH1F*> (fOutput->FindObject("dndeta_check_pt"));

    if (!fdNdEtaAnalysis || !fdNdEtaAnalysisTr || !fdNdEtaAnalysisTrVtx || !fPartPt)
    {
      AliDebug(AliLog::kError, Form("ERROR: Histograms not available %p %p", (void*) fdNdEtaAnalysis, (void*) fPartPt));
      return;
    }

    fdNdEtaAnalysis->Finish(0, -1, AlidNdEtaCorrection::kNone);
    fdNdEtaAnalysisNSD->Finish(0, -1, AlidNdEtaCorrection::kNone);
    fdNdEtaAnalysisTr->Finish(0, -1, AlidNdEtaCorrection::kNone);
    fdNdEtaAnalysisTrVtx->Finish(0, -1, AlidNdEtaCorrection::kNone);
    fdNdEtaAnalysisTracks->Finish(0, -1, AlidNdEtaCorrection::kNone);

    Int_t events = (Int_t) fdNdEtaAnalysis->GetData()->GetEventCorrection()->GetMeasuredHistogram()->Integral();
    fPartPt->Scale(1.0/events);
    fPartPt->Scale(1.0/fPartPt->GetBinWidth(1));

    fout = new TFile("analysis_mc.root","RECREATE");

    fdNdEtaAnalysis->SaveHistograms();
    fdNdEtaAnalysisNSD->SaveHistograms();
    fdNdEtaAnalysisTr->SaveHistograms();
    fdNdEtaAnalysisTrVtx->SaveHistograms();
    fdNdEtaAnalysisTracks->SaveHistograms();
    fPartPt->Write();

    fout->Write();
    fout->Close();

    Printf("Writting result to analysis_mc.root");

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