[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(AliPWG0Helper::kTPC),
  fReadMC(kFALSE),
  fUseMCVertex(kFALSE),
  fEsdTrackCuts(0),
  fdNdEtaAnalysisESD(0),
  fMult(0),
  fMultVtx(0),
  fEvents(0),
  fVertexResolution(0),
  fdNdEtaAnalysis(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());
}

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);
    // PrimaryVertex also needed

    if (fAnalysisMode == AliPWG0Helper::kSPD) {
      tree->SetBranchStatus("fSPDMult*", 1);
      //AliPWG0Helper::SetBranchStatusRecursive(tree, "AliMultiplicity", kTRUE);
    }
    
    if (fAnalysisMode == AliPWG0Helper::kTPC || fAnalysisMode == AliPWG0Helper::kTPCITS) {
      AliESDtrackCuts::EnableNeededBranches(tree);
      tree->SetBranchStatus("fTracks.fLabel", 1);
      //AliPWG0Helper::SetBranchStatusRecursive(tree, "Tracks", kTRUE);
    }

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

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

  // 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", 40, -20, 20);
  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);

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

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

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

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

      fDeltaPhi->Fill(mult->GetDeltaPhi(i));

      etaArr[inputCount] = mult->GetEta(i);
      labelArr[inputCount] = mult->GetLabel(i);
      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;

  if (fUseMCVertex) {
    Printf("WARNING: Replacing vertex by MC vertex. This is for systematical checks only.");
    if (!fReadMC) {
      Printf("ERROR: fUseMCVertex 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;
    }

    // get the MC vertex
    AliGenEventHeader* genHeader = header->GenEventHeader();
    TArrayF vtxMC(3);
    genHeader->PrimaryVertex(vtxMC);

    vtx[2] = vtxMC[2];
  }

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

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