[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);
    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);
770a1f1d	332	if (vtxESD)
dd367a14	333	{
dd367a14	334	// control hist
0f67a57c	335	fMultVtx->Fill(inputCount);
0f67a57c	336
dd367a14	337	for (Int_t i=0; i<inputCount; ++i)
	338	{
	339	Float_t eta = etaArr[i];
	340	Float_t pt = ptArr[i];
0f67a57c	341
dd367a14	342	fdNdEtaAnalysisESD->FillTrack(vtx[2], eta, pt);
0f67a57c	343
dd367a14	344	if (TMath::Abs(vtx[2]) < 20)
	345	{
	346	fPartEta[0]->Fill(eta);
0f67a57c	347
dd367a14	348	if (vtx[2] < 0)
	349	fPartEta[1]->Fill(eta);
	350	else
	351	fPartEta[2]->Fill(eta);
	352	}
0f67a57c	353	}
0f67a57c	354
dd367a14	355	// for event count per vertex
dd367a14	356	fdNdEtaAnalysisESD->FillEvent(vtx[2], inputCount);
0f67a57c	357
54b096ef	358	// control hist
dd367a14	359	fEvents->Fill(vtx[2]);
dd367a14	360	}
0f67a57c	361	}
	362
	363	if (fReadMC) // Processing of MC information (optional)
	364	{
	365	AliMCEventHandler* eventHandler = dynamic_cast<AliMCEventHandler*> (AliAnalysisManager::GetAnalysisManager()->GetMCtruthEventHandler());
	366	if (!eventHandler) {
	367	Printf("ERROR: Could not retrieve MC event handler");
	368	return;
	369	}
	370
	371	AliMCEvent* mcEvent = eventHandler->MCEvent();
	372	if (!mcEvent) {
	373	Printf("ERROR: Could not retrieve MC event");
	374	return;
	375	}
	376
	377	AliStack* stack = mcEvent->Stack();
	378	if (!stack)
	379	{
	380	AliDebug(AliLog::kError, "Stack not available");
	381	return;
	382	}
	383
	384	AliHeader* header = mcEvent->Header();
	385	if (!header)
	386	{
	387	AliDebug(AliLog::kError, "Header not available");
	388	return;
	389	}
	390
	391	// get the MC vertex
	392	AliGenEventHeader* genHeader = header->GenEventHeader();
	393	TArrayF vtxMC(3);
	394	genHeader->PrimaryVertex(vtxMC);
	395
	396	// loop over mc particles
	397	Int_t nPrim = stack->GetNprimary();
	398
	399	Int_t nAcceptedParticles = 0;
	400
	401	for (Int_t iMc = 0; iMc < nPrim; ++iMc)
	402	{
	403	//Printf("Getting particle %d", iMc);
	404	TParticle* particle = stack->Particle(iMc);
	405
	406	if (!particle)
	407	continue;
	408
	409	if (AliPWG0Helper::IsPrimaryCharged(particle, nPrim) == kFALSE)
	410	continue;
	411
	412	AliDebug(AliLog::kDebug+1, Form("Accepted primary %d, unique ID: %d", iMc, particle->GetUniqueID()));
0f67a57c	413	Float_t eta = particle->Eta();
	414	Float_t pt = particle->Pt();
	415
54b096ef	416	// 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	417	if (TMath::Abs(eta) < 1.5 && pt > 0.3)
2b26296f	418	nAcceptedParticles++;
2b26296f	419
0f67a57c	420	fdNdEtaAnalysis->FillTrack(vtxMC[2], eta, pt);
	421	fVertex->Fill(particle->Vx(), particle->Vy(), particle->Vz());
	422
	423	if (eventTriggered)
	424	{
	425	fdNdEtaAnalysisTr->FillTrack(vtxMC[2], eta, pt);
770a1f1d	426	if (vtxESD)
0f67a57c	427	fdNdEtaAnalysisTrVtx->FillTrack(vtxMC[2], eta, pt);
	428	}
	429
	430	if (TMath::Abs(eta) < 0.8)
	431	fPartPt->Fill(pt);
	432	}
	433
	434	fdNdEtaAnalysis->FillEvent(vtxMC[2], nAcceptedParticles);
	435	if (eventTriggered)
	436	{
	437	fdNdEtaAnalysisTr->FillEvent(vtxMC[2], nAcceptedParticles);
770a1f1d	438	if (vtxESD)
0f67a57c	439	fdNdEtaAnalysisTrVtx->FillEvent(vtxMC[2], nAcceptedParticles);
	440	}
	441
770a1f1d	442	if (eventTriggered && vtxESD)
0f67a57c	443	{
	444	// from tracks is only done for triggered and vertex reconstructed events
	445
	446	for (Int_t i=0; i<inputCount; ++i)
	447	{
	448	Int_t label = labelArr[i];
	449
	450	if (label < 0)
	451	continue;
	452
	453	//Printf("Getting particle of track %d", label);
	454	TParticle* particle = stack->Particle(label);
	455
	456	if (!particle)
	457	{
	458	AliDebug(AliLog::kError, Form("ERROR: Could not retrieve particle %d.", label));
	459	continue;
	460	}
	461
	462	fdNdEtaAnalysisTracks->FillTrack(vtxMC[2], particle->Eta(), particle->Pt());
	463	} // end of track loop
	464
	465	// for event count per vertex
	466	fdNdEtaAnalysisTracks->FillEvent(vtxMC[2], inputCount);
	467	}
	468	}
	469
	470	delete[] etaArr;
	471	delete[] labelArr;
	472	delete[] ptArr;
	473	}
	474
	475	void AlidNdEtaTask::Terminate(Option_t *)
	476	{
	477	// The Terminate() function is the last function to be called during
	478	// a query. It always runs on the client, it can be used to present
	479	// the results graphically or save the results to file.
	480
	481	fOutput = dynamic_cast<TList*> (GetOutputData(0));
	482	if (!fOutput) {
	483	Printf("ERROR: fOutput not available");
	484	return;
	485	}
	486
	487	fdNdEtaAnalysisESD = dynamic_cast<dNdEtaAnalysis*> (fOutput->FindObject("fdNdEtaAnalysisESD"));
	488	fMult = dynamic_cast<TH1F*> (fOutput->FindObject("fMult"));
	489	fMultVtx = dynamic_cast<TH1F*> (fOutput->FindObject("fMultVtx"));
	490	fEvents = dynamic_cast<TH1F*> (fOutput->FindObject("dndeta_check_vertex"));
54b096ef	491	fVertexResolution = dynamic_cast<TH1F*> (fOutput->FindObject("dndeta_vertex_resolution_z"));
0f67a57c	492	for (Int_t i=0; i<3; ++i)
	493	fPartEta[i] = dynamic_cast<TH1F*> (fOutput->FindObject(Form("dndeta_check_%d", i)));
	494
	495	if (!fdNdEtaAnalysisESD)
	496	{
	497	AliDebug(AliLog::kError, "ERROR: fdNdEtaAnalysisESD not available");
	498	return;
	499	}
	500
	501	if (fMult && fMultVtx)
	502	{
	503	new TCanvas;
	504	fMult->Draw();
	505	fMultVtx->SetLineColor(2);
770a1f1d	506	fMultVtx->Draw("SAME");
0f67a57c	507	}
	508
	509	if (fPartEta[0])
	510	{
	511	Int_t events1 = (Int_t) fEvents->Integral(fEvents->GetXaxis()->FindBin(-19.9), fEvents->GetXaxis()->FindBin(-0.001));
	512	Int_t events2 = (Int_t) fEvents->Integral(fEvents->GetXaxis()->FindBin(0.001), fEvents->GetXaxis()->FindBin(19.9));
	513
54b096ef	514	Printf("%d events with vertex used", events1 + events2);
	515
	516	if (events1 > 0 && events2 > 0)
	517	{
0f67a57c	518	fPartEta[0]->Scale(1.0 / (events1 + events2));
	519	fPartEta[1]->Scale(1.0 / events1);
	520	fPartEta[2]->Scale(1.0 / events2);
	521
	522	for (Int_t i=0; i<3; ++i)
	523	fPartEta[i]->Scale(1.0 / fPartEta[i]->GetBinWidth(1));
	524
	525	new TCanvas("control", "control", 500, 500);
	526	for (Int_t i=0; i<3; ++i)
	527	{
	528	fPartEta[i]->SetLineColor(i+1);
	529	fPartEta[i]->Draw((i==0) ? "" : "SAME");
	530	}
54b096ef	531	}
0f67a57c	532	}
	533
	534	if (fEvents)
	535	{
	536	new TCanvas("control3", "control3", 500, 500);
770a1f1d	537	fEvents->Draw();
0f67a57c	538	}
	539
	540	TFile* fout = new TFile("analysis_esd_raw.root", "RECREATE");
	541
	542	if (fdNdEtaAnalysisESD)
	543	fdNdEtaAnalysisESD->SaveHistograms();
	544
	545	if (fEsdTrackCuts)
	546	fEsdTrackCuts->SaveHistograms("esd_tracks_cuts");
	547
	548	if (fMult)
	549	fMult->Write();
	550
	551	if (fMultVtx)
	552	fMultVtx->Write();
	553
	554	for (Int_t i=0; i<3; ++i)
	555	if (fPartEta[i])
	556	fPartEta[i]->Write();
	557
	558	if (fEvents)
	559	fEvents->Write();
	560
54b096ef	561	if (fVertexResolution)
	562	fVertexResolution->Write();
	563
	564	if (fDeltaPhi)
	565	fDeltaPhi->Write();
	566
0f67a57c	567	fout->Write();
	568	fout->Close();
	569
	570	Printf("Writting result to analysis_esd_raw.root");
	571
	572	if (fReadMC)
	573	{
	574	fdNdEtaAnalysis = dynamic_cast<dNdEtaAnalysis*> (fOutput->FindObject("dndeta"));
	575	fdNdEtaAnalysisTr = dynamic_cast<dNdEtaAnalysis*> (fOutput->FindObject("dndetaTr"));
	576	fdNdEtaAnalysisTrVtx = dynamic_cast<dNdEtaAnalysis*> (fOutput->FindObject("dndetaTrVtx"));
	577	fdNdEtaAnalysisTracks = dynamic_cast<dNdEtaAnalysis*> (fOutput->FindObject("dndetaTracks"));
	578	fPartPt = dynamic_cast<TH1F*> (fOutput->FindObject("dndeta_check_pt"));
	579
	580	if (!fdNdEtaAnalysis \|\| !fdNdEtaAnalysisTr \|\| !fdNdEtaAnalysisTrVtx \|\| !fPartPt)
	581	{
	582	AliDebug(AliLog::kError, Form("ERROR: Histograms not available %p %p", (void) fdNdEtaAnalysis, (void) fPartPt));
	583	return;
	584	}
	585
	586	fdNdEtaAnalysis->Finish(0, -1, AlidNdEtaCorrection::kNone);
	587	fdNdEtaAnalysisTr->Finish(0, -1, AlidNdEtaCorrection::kNone);
	588	fdNdEtaAnalysisTrVtx->Finish(0, -1, AlidNdEtaCorrection::kNone);
	589	fdNdEtaAnalysisTracks->Finish(0, -1, AlidNdEtaCorrection::kNone);
	590
	591	Int_t events = (Int_t) fdNdEtaAnalysis->GetData()->GetEventCorrection()->GetMeasuredHistogram()->Integral();
	592	fPartPt->Scale(1.0/events);
	593	fPartPt->Scale(1.0/fPartPt->GetBinWidth(1));
	594
	595	TFile* fout = new TFile("analysis_mc.root","RECREATE");
	596
	597	fdNdEtaAnalysis->SaveHistograms();
	598	fdNdEtaAnalysisTr->SaveHistograms();
	599	fdNdEtaAnalysisTrVtx->SaveHistograms();
	600	fdNdEtaAnalysisTracks->SaveHistograms();
	601	fPartPt->Write();
	602
	603	fout->Write();
	604	fout->Close();
	605
	606	Printf("Writting result to analysis_mc.root");
	607
	608	if (fPartPt)
	609	{
	610	new TCanvas("control2", "control2", 500, 500);
770a1f1d	611	fPartPt->Draw();
0f67a57c	612	}
	613	}
	614	}