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

/* $Id$ */

#include "AlidNdEtaCorrectionTask.h"

#include <TCanvas.h>
#include <TChain.h>
#include <TFile.h>
#include <TH1F.h>
#include <TParticle.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"
#include "dNdEta/AlidNdEtaCorrection.h"

ClassImp(AlidNdEtaCorrectionTask)

AlidNdEtaCorrectionTask::AlidNdEtaCorrectionTask(const char* opt) :
  AliAnalysisTask("AlidNdEtaCorrectionTask", ""),
  fESD(0),
  fOutput(0),
  fOption(opt),
  fAnalysisMode(kTPC),
  fSignMode(0),
  fEsdTrackCuts(0),
  fdNdEtaCorrection(0),
  fdNdEtaAnalysisMC(0),
  fdNdEtaAnalysisESD(0),
  fPIDParticles(0),
  fPIDTracks(0)
{
  //
  // Constructor. Initialization of pointers
  //

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

AlidNdEtaCorrectionTask::~AlidNdEtaCorrectionTask()
{
  //
  // 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 AlidNdEtaCorrectionTask::ConnectInputData(Option_t *)
{
  // Connect ESD
  // Called once

  Printf("AlidNdEtaCorrectionTask::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 AlidNdEtaCorrectionTask::CreateOutputObjects()
{
  // create result objects and add to output list

  if (fOption.Contains("only-positive"))
  {
    Printf("INFO: Processing only positive particles.");
    fSignMode = 1;
  }
  else if (fOption.Contains("only-negative"))
  {
    Printf("INFO: Processing only negative particles.");
    fSignMode = -1;
  }

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

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

  fdNdEtaCorrection = new AlidNdEtaCorrection("dndeta_correction", "dndeta_correction", detector);
  fOutput->Add(fdNdEtaCorrection);

  fPIDParticles = new TH1F("pid_particles", "PID of generated primary particles", 10001, -5000.5, 5000.5);
  fOutput->Add(fPIDParticles);

  fPIDTracks = new TH1F("pid_tracks", "MC PID of reconstructed tracks", 10001, -5000.5, 5000.5);
  fOutput->Add(fPIDTracks);

  fdNdEtaAnalysisMC = new dNdEtaAnalysis("dndetaMC", "dndetaMC", detector);
  fOutput->Add(fdNdEtaAnalysisMC);

  fdNdEtaAnalysisESD = new dNdEtaAnalysis("dndetaESD", "dndetaESD", detector);
  fOutput->Add(fdNdEtaAnalysisESD);
}

void AlidNdEtaCorrectionTask::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());

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

  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 process type
  Int_t processType = AliPWG0Helper::GetPythiaEventProcessType(header);
  AliDebug(AliLog::kDebug+1, Form("Found pythia procces type %d", processType));

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

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

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

  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;

    if (SignOK(particle->GetPDG()) == kFALSE)
      continue;

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

    fdNdEtaCorrection->FillMCParticle(vtxMC[2], eta, pt, eventTriggered, eventVertex, processType);

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

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

    if (label < 0)
    {
      AliDebug(AliLog::kWarning, Form("WARNING: cannot find corresponding mc part for track %d.", label));
      continue;
    }

    TParticle* particle = stack->Particle(label);
    if (!particle)
    {
      AliDebug(AliLog::kError, Form("UNEXPECTED: particle with label %d not found in stack (track loop).", label));
      continue;
    }

    if (SignOK(particle->GetPDG()) == kFALSE)
        continue;

    if (eventTriggered && eventVertex)
    {
      fdNdEtaCorrection->FillTrackedParticle(vtxMC[2], particle->Eta(), particle->Pt());
      fdNdEtaAnalysisESD->FillTrack(vtxMC[2], particle->Eta(), particle->Pt());
      if (particle->Pt() > 0.1 && particle->Pt() < 0.2)
      {
        fPIDTracks->Fill(particle->GetPdgCode());
      }
    }
  }

  if (eventTriggered && eventVertex)
    fdNdEtaAnalysisMC->FillEvent(vtxMC[2], inputCount);

  // stuff regarding the vertex reco correction and trigger bias correction
  fdNdEtaCorrection->FillEvent(vtxMC[2], inputCount, eventTriggered, eventVertex, processType);
  if (eventTriggered) {
    if (eventVertex)
    {
      fdNdEtaAnalysisESD->FillEvent(vtxMC[2], inputCount);
    }
  }

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

void AlidNdEtaCorrectionTask::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;
  }

  fdNdEtaCorrection = dynamic_cast<AlidNdEtaCorrection*> (fOutput->FindObject("dndeta_correction"));
  fdNdEtaAnalysisMC = dynamic_cast<dNdEtaAnalysis*> (fOutput->FindObject("dndetaMC"));
  fdNdEtaAnalysisESD = dynamic_cast<dNdEtaAnalysis*> (fOutput->FindObject("dndetaESD"));
  if (!fdNdEtaCorrection || !fdNdEtaAnalysisMC || !fdNdEtaAnalysisESD)
  {
    AliDebug(AliLog::kError, "Could not read object from output list");
    return;
  }

  fdNdEtaCorrection->Finish();

  TString fileName;
  fileName.Form("correction_map%s.root", fOption.Data());
  TFile* fout = new TFile(fileName, "RECREATE");

  if (fEsdTrackCuts)
    fEsdTrackCuts->SaveHistograms("esd_track_cuts");
  fdNdEtaCorrection->SaveHistograms();
  fdNdEtaAnalysisMC->SaveHistograms();
  fdNdEtaAnalysisESD->SaveHistograms();

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

  fdNdEtaCorrection->DrawHistograms();

  if (fPIDParticles && fPIDTracks)
  {
    new TCanvas("pidcanvas", "pidcanvas", 500, 500);

    fPIDParticles->Draw();
    fPIDTracks->SetLineColor(2);
    fPIDTracks->Draw("SAME");

    TDatabasePDG* pdgDB = new TDatabasePDG;

    for (Int_t i=0; i <= fPIDParticles->GetNbinsX()+1; ++i)
      if (fPIDParticles->GetBinContent(i) > 0)
        printf("PDG = %d (%s): generated: %d, reconstructed: %d, ratio: %f\n", (Int_t) fPIDParticles->GetBinCenter(i), pdgDB->GetParticle((Int_t) fPIDParticles->GetBinCenter(i))->GetName(), (Int_t) fPIDParticles->GetBinContent(i), (Int_t) fPIDTracks->GetBinContent(i), ((fPIDTracks->GetBinContent(i) > 0) ? fPIDParticles->GetBinContent(i) / fPIDTracks->GetBinContent(i) : -1));

    delete pdgDB;
    pdgDB = 0;
  }

  Printf("Writting result to %s", fileName.Data());
}

Bool_t AlidNdEtaCorrectionTask::SignOK(TParticlePDG* particle)
{
  // returns if a particle with this sign should be counted
  // this is determined by the value of fSignMode, which should have the same sign
  // as the charge
  // if fSignMode is 0 all particles are counted

  if (fSignMode == 0)
    return kTRUE;

  if (!particle)
  {
    printf("WARNING: not counting a particle that does not have a pdg particle\n");
    return kFALSE;
  }

  Double_t charge = particle->Charge();

  if (fSignMode > 0)
    if (charge < 0)
      return kFALSE;

  if (fSignMode < 0)
    if (charge > 0)
      return kFALSE;

  return kTRUE;
}
Commit	Line	Data
0f67a57c	1	/* $Id$ */
	2
	3	#include "AlidNdEtaCorrectionTask.h"
	4
	5	#include <TCanvas.h>
	6	#include <TChain.h>
	7	#include <TFile.h>
	8	#include <TH1F.h>
	9	#include <TParticle.h>
	10
	11	#include <AliLog.h>
	12	#include <AliESDVertex.h>
	13	#include <AliESDEvent.h>
	14	#include <AliStack.h>
	15	#include <AliHeader.h>
	16	#include <AliGenEventHeader.h>
	17	#include <AliMultiplicity.h>
	18	#include <AliAnalysisManager.h>
	19	#include <AliMCEventHandler.h>
	20	#include <AliMCEvent.h>
	21	#include <AliESDInputHandler.h>
	22
	23	#include "esdTrackCuts/AliESDtrackCuts.h"
	24	#include "AliPWG0Helper.h"
	25	//#include "AliCorrection.h"
	26	//#include "AliCorrectionMatrix3D.h"
	27	#include "dNdEta/dNdEtaAnalysis.h"
	28	#include "dNdEta/AlidNdEtaCorrection.h"
	29
	30	ClassImp(AlidNdEtaCorrectionTask)
	31
	32	AlidNdEtaCorrectionTask::AlidNdEtaCorrectionTask(const char* opt) :
	33	AliAnalysisTask("AlidNdEtaCorrectionTask", ""),
	34	fESD(0),
	35	fOutput(0),
	36	fOption(opt),
	37	fAnalysisMode(kTPC),
	38	fSignMode(0),
	39	fEsdTrackCuts(0),
	40	fdNdEtaCorrection(0),
	41	fdNdEtaAnalysisMC(0),
	42	fdNdEtaAnalysisESD(0),
	43	fPIDParticles(0),
	44	fPIDTracks(0)
	45	{
	46	//
	47	// Constructor. Initialization of pointers
	48	//
	49
	50	// Define input and output slots here
	51	DefineInput(0, TChain::Class());
	52	DefineOutput(0, TList::Class());
	53	}
	54
	55	AlidNdEtaCorrectionTask::~AlidNdEtaCorrectionTask()
	56	{
	57	//
	58	// Destructor
	59	//
	60
	61	// histograms are in the output list and deleted when the output
	62	// list is deleted by the TSelector dtor
	63
	64	if (fOutput) {
65	delete fOutput;
66	fOutput = 0;
67	}
68	}
69
70	//________________________________________________________________________
71	void AlidNdEtaCorrectionTask::ConnectInputData(Option_t *)
72	{
73	// Connect ESD
74	// Called once
75
76	Printf("AlidNdEtaCorrectionTask::ConnectInputData called");
77
78	TTree* tree = dynamic_cast<TTree*> (GetInputData(0));
79	if (!tree) {
80	Printf("ERROR: Could not read tree from input slot 0");
81	} else {
82	// Disable all branches and enable only the needed ones
83	tree->SetBranchStatus("*", 0);
84
85	tree->SetBranchStatus("fTriggerMask", 1);
86	tree->SetBranchStatus("fSPDVertex*", 1);
87
88	if (fAnalysisMode == kSPD)
89	tree->SetBranchStatus("fSPDMult*", 1);
90
91	if (fAnalysisMode == kTPC) {
92	AliESDtrackCuts::EnableNeededBranches(tree);
93	tree->SetBranchStatus("fTracks.fLabel", 1);
94	}
95
0f67a57c	96	AliESDInputHandler esdH = dynamic_cast<AliESDInputHandler> (AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler());
	97
	98	if (!esdH) {
	99	Printf("ERROR: Could not get ESDInputHandler");
	100	} else
	101	fESD = esdH->GetEvent();
	102	}
	103	}
	104
	105	void AlidNdEtaCorrectionTask::CreateOutputObjects()
	106	{
	107	// create result objects and add to output list
	108
	109	if (fOption.Contains("only-positive"))
	110	{
	111	Printf("INFO: Processing only positive particles.");
	112	fSignMode = 1;
	113	}
	114	else if (fOption.Contains("only-negative"))
	115	{
	116	Printf("INFO: Processing only negative particles.");
	117	fSignMode = -1;
	118	}
	119
	120	TString detector;
	121	if (fAnalysisMode == kTPC)
	122	{
	123	detector = "TPC";
	124	}
	125	else if (fAnalysisMode == kSPD)
	126	detector = "SPD";
	127
	128	fOutput = new TList;
	129	fOutput->SetOwner();
	130
	131	fdNdEtaCorrection = new AlidNdEtaCorrection("dndeta_correction", "dndeta_correction", detector);
	132	fOutput->Add(fdNdEtaCorrection);
	133
	134	fPIDParticles = new TH1F("pid_particles", "PID of generated primary particles", 10001, -5000.5, 5000.5);
	135	fOutput->Add(fPIDParticles);
	136
	137	fPIDTracks = new TH1F("pid_tracks", "MC PID of reconstructed tracks", 10001, -5000.5, 5000.5);
	138	fOutput->Add(fPIDTracks);
	139
	140	fdNdEtaAnalysisMC = new dNdEtaAnalysis("dndetaMC", "dndetaMC", detector);
	141	fOutput->Add(fdNdEtaAnalysisMC);
	142
	143	fdNdEtaAnalysisESD = new dNdEtaAnalysis("dndetaESD", "dndetaESD", detector);
	144	fOutput->Add(fdNdEtaAnalysisESD);
	145	}
	146
	147	void AlidNdEtaCorrectionTask::Exec(Option_t*)
	148	{
	149	// process the event
	150
	151	// Check prerequisites
	152	if (!fESD)
	153	{
	154	AliDebug(AliLog::kError, "ESD branch not available");
	155	return;
	156	}
	157
	158	// trigger definition
	159	Bool_t eventTriggered = AliPWG0Helper::IsEventTriggered(fESD->GetTriggerMask(), AliPWG0Helper::kMB2);
160
161	Bool_t eventVertex = AliPWG0Helper::IsVertexReconstructed(fESD->GetVertex());
162
163	// post the data already here
164	PostData(0, fOutput);
165
166	// create list of (label, eta, pt) tuples
167	Int_t inputCount = 0;
168	Int_t* labelArr = 0;
169	Float_t* etaArr = 0;
170	Float_t* ptArr = 0;
171	if (fAnalysisMode == kSPD)
172	{
173	// get tracklets
174	const AliMultiplicity* mult = fESD->GetMultiplicity();
175	if (!mult)
176	{
177	AliDebug(AliLog::kError, "AliMultiplicity not available");
178	return;
179	}
180
181	labelArr = new Int_t[mult->GetNumberOfTracklets()];
182	etaArr = new Float_t[mult->GetNumberOfTracklets()];
183	ptArr = new Float_t[mult->GetNumberOfTracklets()];
184
185	// get multiplicity from ITS tracklets
186	for (Int_t i=0; i<mult->GetNumberOfTracklets(); ++i)
187	{
188	//printf("%d %f %f %f\n", i, mult->GetTheta(i), mult->GetPhi(i), mult->GetDeltaPhi(i));
189
190	// This removes non-tracklets in PDC06 data. Very bad solution. New solution is implemented for newer data. Keeping this for compatibility.
191	if (mult->GetDeltaPhi(i) < -1000)
192	continue;
193
194	etaArr[inputCount] = mult->GetEta(i);
195	labelArr[inputCount] = mult->GetLabel(i);
196	ptArr[inputCount] = 0; // no pt for tracklets
197	++inputCount;
198	}
199	}
200	else if (fAnalysisMode == kTPC)
201	{
202	if (!fEsdTrackCuts)
203	{
204	AliDebug(AliLog::kError, "fESDTrackCuts not available");
205	return;
206	}
207
208	// get multiplicity from ESD tracks
209	TObjArray* list = fEsdTrackCuts->GetAcceptedTracks(fESD);
210	Int_t nGoodTracks = list->GetEntries();
211
212	labelArr = new Int_t[nGoodTracks];
213	etaArr = new Float_t[nGoodTracks];
214	ptArr = new Float_t[nGoodTracks];
215
216	// loop over esd tracks
217	for (Int_t i=0; i<nGoodTracks; i++)
218	{
219	AliESDtrack* esdTrack = dynamic_cast<AliESDtrack*> (list->At(i));
220	if (!esdTrack)
221	{
222	AliDebug(AliLog::kError, Form("ERROR: Could not retrieve track %d.", i));
223	continue;
224	}
225
226	etaArr[inputCount] = esdTrack->Eta();
227	labelArr[inputCount] = TMath::Abs(esdTrack->GetLabel());
228	ptArr[inputCount] = esdTrack->Pt();
229	++inputCount;
230	}
231	}
232	else
233	return;
234
235	AliMCEventHandler* eventHandler = dynamic_cast<AliMCEventHandler*> (AliAnalysisManager::GetAnalysisManager()->GetMCtruthEventHandler());
236	if (!eventHandler) {
237	Printf("ERROR: Could not retrieve MC event handler");
238	return;
239	}
240
241	AliMCEvent* mcEvent = eventHandler->MCEvent();
242	if (!mcEvent) {
243	Printf("ERROR: Could not retrieve MC event");
244	return;
245	}
246
247	AliStack* stack = mcEvent->Stack();
248	if (!stack)
249	{
250	AliDebug(AliLog::kError, "Stack not available");
251	return;
252	}
253
254	AliHeader* header = mcEvent->Header();
255	if (!header)
256	{
257	AliDebug(AliLog::kError, "Header not available");
258	return;
259	}
260
261	// get process type
262	Int_t processType = AliPWG0Helper::GetPythiaEventProcessType(header);
263	AliDebug(AliLog::kDebug+1, Form("Found pythia procces type %d", processType));
264
265	if (processType<0)
266	AliDebug(AliLog::kError, Form("Unknown Pythia process type %d.", processType));
267
268	// get the MC vertex
269	AliGenEventHeader* genHeader = header->GenEventHeader();
270	TArrayF vtxMC(3);
271	genHeader->PrimaryVertex(vtxMC);
272
273	// loop over mc particles
274	Int_t nPrim = stack->GetNprimary();
275
276	for (Int_t iMc = 0; iMc < nPrim; ++iMc)
277	{
278	//Printf("Getting particle %d", iMc);
279	TParticle* particle = stack->Particle(iMc);
280
281	if (!particle)
282	continue;
283
284	if (AliPWG0Helper::IsPrimaryCharged(particle, nPrim) == kFALSE)
285	continue;
286
287	if (SignOK(particle->GetPDG()) == kFALSE)
288	continue;
289
290	Float_t eta = particle->Eta();
291	Float_t pt = particle->Pt();
292
293	fdNdEtaCorrection->FillMCParticle(vtxMC[2], eta, pt, eventTriggered, eventVertex, processType);
294
295	if (eventTriggered)
296	if (eventVertex)
297	fdNdEtaAnalysisMC->FillTrack(vtxMC[2], eta, pt);
298	}
299
300	for (Int_t i=0; i<inputCount; ++i)
301	{
302	Int_t label = labelArr[i];
303
304	if (label < 0)
305	{
306	AliDebug(AliLog::kWarning, Form("WARNING: cannot find corresponding mc part for track %d.", label));
307	continue;
308	}
309
310	TParticle* particle = stack->Particle(label);
311	if (!particle)
312	{
313	AliDebug(AliLog::kError, Form("UNEXPECTED: particle with label %d not found in stack (track loop).", label));
314	continue;
315	}
316
317	if (SignOK(particle->GetPDG()) == kFALSE)
318	continue;
319
320	if (eventTriggered && eventVertex)
321	{
322	fdNdEtaCorrection->FillTrackedParticle(vtxMC[2], particle->Eta(), particle->Pt());
323	fdNdEtaAnalysisESD->FillTrack(vtxMC[2], particle->Eta(), particle->Pt());
324	if (particle->Pt() > 0.1 && particle->Pt() < 0.2)
325	{
326	fPIDTracks->Fill(particle->GetPdgCode());
327	}
328	}
329	}
330
331	if (eventTriggered && eventVertex)
332	fdNdEtaAnalysisMC->FillEvent(vtxMC[2], inputCount);
333
334	// stuff regarding the vertex reco correction and trigger bias correction
335	fdNdEtaCorrection->FillEvent(vtxMC[2], inputCount, eventTriggered, eventVertex, processType);
336	if (eventTriggered) {
337	if (eventVertex)
338	{
339	fdNdEtaAnalysisESD->FillEvent(vtxMC[2], inputCount);
340	}
341	}
342
343	delete[] etaArr;
344	delete[] labelArr;
345	delete[] ptArr;
346	}
347
348	void AlidNdEtaCorrectionTask::Terminate(Option_t *)
349	{
350	// The Terminate() function is the last function to be called during
351	// a query. It always runs on the client, it can be used to present
352	// the results graphically or save the results to file.
353
354	fOutput = dynamic_cast<TList*> (GetOutputData(0));
355	if (!fOutput) {
356	Printf("ERROR: fOutput not available");
357	return;
358	}
359
360	fdNdEtaCorrection = dynamic_cast<AlidNdEtaCorrection*> (fOutput->FindObject("dndeta_correction"));
361	fdNdEtaAnalysisMC = dynamic_cast<dNdEtaAnalysis*> (fOutput->FindObject("dndetaMC"));
362	fdNdEtaAnalysisESD = dynamic_cast<dNdEtaAnalysis*> (fOutput->FindObject("dndetaESD"));
363	if (!fdNdEtaCorrection \|\| !fdNdEtaAnalysisMC \|\| !fdNdEtaAnalysisESD)
364	{
365	AliDebug(AliLog::kError, "Could not read object from output list");
366	return;
367	}
368
369	fdNdEtaCorrection->Finish();
370
371	TString fileName;
372	fileName.Form("correction_map%s.root", fOption.Data());
373	TFile* fout = new TFile(fileName, "RECREATE");
374
375	if (fEsdTrackCuts)
376	fEsdTrackCuts->SaveHistograms("esd_track_cuts");
377	fdNdEtaCorrection->SaveHistograms();
378	fdNdEtaAnalysisMC->SaveHistograms();
379	fdNdEtaAnalysisESD->SaveHistograms();
380
381	fout->Write();
382	fout->Close();
383
384	fdNdEtaCorrection->DrawHistograms();
385
386	if (fPIDParticles && fPIDTracks)
387	{
388	new TCanvas("pidcanvas", "pidcanvas", 500, 500);
389
390	fPIDParticles->Draw();
391	fPIDTracks->SetLineColor(2);
392	fPIDTracks->Draw("SAME");
393
394	TDatabasePDG* pdgDB = new TDatabasePDG;
395
396	for (Int_t i=0; i <= fPIDParticles->GetNbinsX()+1; ++i)
397	if (fPIDParticles->GetBinContent(i) > 0)
398	printf("PDG = %d (%s): generated: %d, reconstructed: %d, ratio: %f\n", (Int_t) fPIDParticles->GetBinCenter(i), pdgDB->GetParticle((Int_t) fPIDParticles->GetBinCenter(i))->GetName(), (Int_t) fPIDParticles->GetBinContent(i), (Int_t) fPIDTracks->GetBinContent(i), ((fPIDTracks->GetBinContent(i) > 0) ? fPIDParticles->GetBinContent(i) / fPIDTracks->GetBinContent(i) : -1));
399
400	delete pdgDB;
401	pdgDB = 0;
402	}
403
404	Printf("Writting result to %s", fileName.Data());
405	}
406
407	Bool_t AlidNdEtaCorrectionTask::SignOK(TParticlePDG* particle)
408	{
409	// returns if a particle with this sign should be counted
410	// this is determined by the value of fSignMode, which should have the same sign
411	// as the charge
412	// if fSignMode is 0 all particles are counted
413
414	if (fSignMode == 0)
415	return kTRUE;
416
417	if (!particle)
418	{
419	printf("WARNING: not counting a particle that does not have a pdg particle\n");
420	return kFALSE;
421	}
422
423	Double_t charge = particle->Charge();
424
425	if (fSignMode > 0)
426	if (charge < 0)
427	return kFALSE;
428
429	if (fSignMode < 0)
430	if (charge > 0)
431	return kFALSE;
432
433	return kTRUE;
434	}