[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 "esdTrackCuts/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(kTPC),
  fReadMC(kFALSE),
  fEsdTrackCuts(0),
  fdNdEtaAnalysisESD(0),
  fMult(0),
  fMultVtx(0),
  fEvents(0),
  fdNdEtaAnalysis(0),
  fdNdEtaAnalysisTr(0),
  fdNdEtaAnalysisTrVtx(0),
  fdNdEtaAnalysisTracks(0),
  fVertex(0),
  fPartPt(0)
{
  //
  // Constructor. Initialization of pointers
  //

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

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;
  }
}

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

  Printf("AlidNdEtaTask::ConnectInputData called");

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

    tree->SetBranchStatus("fTriggerMask", 1);
    tree->SetBranchStatus("fSPDVertex*", 1);

    if (fAnalysisMode == kSPD)
      tree->SetBranchStatus("fSPDMult*", 1);

    if (fAnalysisMode == kTPC) {
      AliESDtrackCuts::EnableNeededBranches(tree);
      tree->SetBranchStatus("fTracks.fLabel", 1);
    }

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

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

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

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

  TString detector;
  if (fAnalysisMode == kTPC)
  {
    detector = "TPC";
  }
  else if (fAnalysisMode == kSPD)
    detector = "SPD";

  fdNdEtaAnalysisESD = new dNdEtaAnalysis("fdNdEtaAnalysisESD", "fdNdEtaAnalysisESD", detector);
  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", 40, -20, 20);
  fOutput->Add(fEvents);

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

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

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

    fdNdEtaAnalysisTracks = new dNdEtaAnalysis("dndetaTracks", "dndetaTracks", detector);
    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);
  }
}

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(), AliPWG0Helper::kMB2);

  Bool_t eventVertex = AliPWG0Helper::IsVertexReconstructed(fESD->GetVertex());

  // get the ESD vertex
  const AliESDVertex* vtxESD = fESD->GetVertex();
  Double_t vtx[3];
  vtxESD->GetXYZ(vtx);

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

  // 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 == 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()];

    // 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));

      // This removes non-tracklets in PDC06 data. Very bad solution. New solution is implemented for newer data. Keeping this for compatibility.
      if (mult->GetDeltaPhi(i) < -1000)
        continue;

      etaArr[inputCount] = mult->GetEta(i);
      labelArr[inputCount] = mult->GetLabel(i);
      ptArr[inputCount] = 0; // no pt for tracklets
      ++inputCount;
    }
  }
  else if (fAnalysisMode == kTPC)
  {
    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("%f", esdTrack->Pt());
    }
  }
  else
    return;

  // Processing of ESD information (always)
  if (eventTriggered)
  {
    // control hists
    fMult->Fill(inputCount);
    if (eventVertex)
      fMultVtx->Fill(inputCount);

    // count all events
    if (!eventVertex)
      vtx[2] = 0;

    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 hists
    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;
    }

    AliStack* 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();
    TArrayF vtxMC(3);
    genHeader->PrimaryVertex(vtxMC);

    // 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()));
      ++nAcceptedParticles;

      Float_t eta = particle->Eta();
      Float_t pt = particle->Pt();

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

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

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

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

    if (eventTriggered && eventVertex)
    {
      // 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"));
  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->DrawCopy("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));

    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->DrawCopy();
  }

  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();

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

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

  if (fReadMC)
  {
    fdNdEtaAnalysis = dynamic_cast<dNdEtaAnalysis*> (fOutput->FindObject("dndeta"));
    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);
    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));

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

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