[u/mrichter/AliRoot.git] / PWGCF / EBYE / PIDFluctuation / task / AliEbyEPidTaskFastGen.cxx

/**************************************************************************
 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
 *                                                                        *
 * Author: Satyajit Jena.                                                 *
 * Contributors are mentioned in the code where appropriate.              *
 *                                                                        *
 * Permission to use, copy, modify and distribute this software and its   *
 * documentation strictly for non-commercial purposes is hereby granted   *
 * without fee, provided that the above copyright notice appears in all   *
 * copies and that both the copyright notice and this permission notice   *
 * appear in the supporting documentation. The authors make no claims     *
 * about the suitability of this software for any purpose. It is          *
 * provided "as is" without express or implied warranty.                  *
 **************************************************************************/


//=========================================================================//
//                 AliEbyE Fluctuaion Analysis for PID                     //
//                   Deepika Jena | drathee@cern.ch                        //
//                   Satyajit Jena | sjena@cern.ch                         //
//=========================================================================//

#include "TChain.h"
#include "TList.h"
#include "TFile.h"
#include "TTree.h"
#include "TH1D.h"
#include "TH2F.h"
#include "TH3F.h"

#include "THnSparse.h"

#include "AliAnalysisTask.h"
#include "AliAnalysisManager.h"

#include "AliVEvent.h"
#include "AliESDEvent.h"
#include "AliMCEvent.h"
#include "AliAODEvent.h"

#include "AliStack.h"
#include "AliGenEventHeader.h"
#include "AliGenPythiaEventHeader.h"
#include "AliGenHijingEventHeader.h"
#include "AliGenDPMjetEventHeader.h"
#include "TDatabasePDG.h"

#include "AliEbyEPidTaskFastGen.h"

ClassImp(AliEbyEPidTaskFastGen)

//-----------------------------------------------------------------------
AliEbyEPidTaskFastGen::AliEbyEPidTaskFastGen( const char *name )
: AliAnalysisTaskSE( name ),
  fThnList(0),
  fVxMax(3.),
  fVyMax(3.),
  fVzMax(15.),
  fPtLowerLimit(0.2),
  fPtHigherLimit(5.),
  fEtaLowerLimit(-1.),
  fEtaHigherLimit(1.),
  fHistoCorrelationMC(NULL)
{ 
 
  DefineOutput(1, TList::Class()); 
}

AliEbyEPidTaskFastGen::~AliEbyEPidTaskFastGen()
{
  if (fThnList) delete fThnList;
}


//---------------------------------------------------------------------------------
void AliEbyEPidTaskFastGen::UserCreateOutputObjects() {
  fThnList = new TList();
  fThnList->SetOwner(kTRUE);
  Int_t fgSparseDataBins[kNSparseData]   = {50000,   500,     440,  30000,   15000,  15000,   10000,  5000,   5000,  5000,  2500,  2500,  2500,  1250, 1250};
  Double_t fgSparseDataMin[kNSparseData] = {0.,       0.,      0.,     0.,      0.,     0.,      0.,    0.,     0.,    0.,    0.,    0.,    0.,   0.,    0.};
  Double_t fgSparseDataMax[kNSparseData] = {50000.,  100.,    440., 30000.,  15000., 15000.,  10000., 5000.,  5000., 5000., 2500., 2500., 2500., 1250., 1250.};
  
  const Char_t *fgkSparseDataTitle[] = {
    "Nnumber of Partiles",
    "Impact Parameters",
    "N_{p}",
    "N_{ch}",
    "N_{+}",
    "N_{-}",
    "N_{#pi}",
    "N_{#pi^{+}}",
    "N_{#pi^{-}}",
    "N_{K}",
    "N_{K^{+}}",
    "N_{K^{-}}",
    "N_{pr}",
    "N_{p}",
    "N_{#bar{p}}"
  };
  
  fHistoCorrelationMC = new THnSparseI("hHistoPR", "", kNSparseData, fgSparseDataBins, fgSparseDataMin, fgSparseDataMax);
  for (Int_t iaxis = 0; iaxis < kNSparseData; iaxis++)
    fHistoCorrelationMC->GetAxis(iaxis)->SetTitle(fgkSparseDataTitle[iaxis]);
  fThnList->Add(fHistoCorrelationMC);
  
  PostData(1, fThnList);
}


//----------------------------------------------------------------------------------
void AliEbyEPidTaskFastGen::UserExec( Option_t * ){
  AliMCEvent* mcEvent = MCEvent();
  if (!mcEvent) {
    Printf("ERROR: Could not retrieve MC event");
    return;
  }
  
  AliStack *stack = mcEvent->Stack();
 
  AliGenEventHeader* genHeader = mcEvent->GenEventHeader();
  if(!genHeader){
    printf("  Event generator header not available!!!\n");
    return;
  }
  
  Double_t imp = -1;
  Double_t np = -1;
  Double_t nt = -1;
  
  if(genHeader->InheritsFrom(AliGenHijingEventHeader::Class())){
    imp = ((AliGenHijingEventHeader*) genHeader)->ImpactParameter();
    nt  = ((AliGenHijingEventHeader*) genHeader)->TargetParticipants();
    np  = ((AliGenHijingEventHeader*) genHeader)->ProjectileParticipants();
  }  
  

  Int_t count[4][2];
  for (Int_t i = 0; i < 4; i++) {
    for (Int_t j = 0; j < 2; j++) count[i][j] = 0; 
  }
  
  Int_t nParticles = stack->GetNtrack();
  for (Int_t iParticle = 0; iParticle < nParticles; iParticle++) {
    TParticle* part = stack->Particle(iParticle);
    if (!part) {
      Printf(" No Particle Available ");
      continue;
    }
    
    //  TDatabasePDG* pdgDB   = TDatabasePDG::Instance();
    //  TParticlePDG* pdgPart = pdgDB->GetParticle(part->GetPdgCode());
    
    Float_t pt = part->Pt();
    if(pt < fPtLowerLimit || pt > fPtHigherLimit ) continue;
    Float_t gEta = part->Eta();
    if(gEta < fEtaLowerLimit || gEta > fEtaHigherLimit ) continue;
    
    
    //  Float_t gCharge = (pdgPart ? pdgPart->Charge() : 0);
    
    Int_t pid = part->GetPdgCode();
    if(pid == 211)         { count[1][0]++;  count[0][0]++; } 
    else if(pid ==  -211)  { count[1][1]++;  count[0][1]++; } 
    else if(pid ==   321)  { count[2][0]++;  count[0][0]++; } 
    else if(pid ==  -321)  { count[2][1]++;  count[0][1]++; } 
    else if(pid ==  2212)  { count[3][0]++;  count[0][0]++; } 
    else if(pid == -2212)  { count[3][1]++;  count[0][1]++; } 
  }
  
  Double_t vsparseMC[kNSparseData];
  vsparseMC[0]  = nParticles;
  vsparseMC[1]  = imp;
  vsparseMC[2]  = np + nt;
  vsparseMC[3]  = count[0][0] + count[0][1];
  vsparseMC[4]  = count[0][0];
  vsparseMC[5]  = count[0][1];

  vsparseMC[6]  = count[1][0] + count[1][1];
  vsparseMC[7]  = count[1][0];
  vsparseMC[8]  = count[1][1];

  vsparseMC[9]  = count[2][0] + count[2][1];
  vsparseMC[10]  = count[2][0];
  vsparseMC[11]  = count[2][1];

  vsparseMC[12] = count[3][0] + count[3][1];
  vsparseMC[13] = count[3][0];
  vsparseMC[14] = count[3][1];
  

   /* Printf(" %10.2f %10.2f %10.2f %10.2f %10.2f %10.2f %10.2f %10.2f %10.2f %10.2f %10.2f %10.2f %10.2f %10.2f %10.2f",
  	 vsparseMC[0], vsparseMC[1], vsparseMC[2], 
	 vsparseMC[3], vsparseMC[4], vsparseMC[5], 
	 vsparseMC[6], vsparseMC[7], vsparseMC[8], 
	 vsparseMC[9], vsparseMC[10], vsparseMC[11],
	 vsparseMC[12], vsparseMC[13], vsparseMC[13]);
   */
  fHistoCorrelationMC->Fill(vsparseMC);
  PostData(1, fThnList);
  
}

void AliEbyEPidTaskFastGen::Terminate( Option_t * ){

  Info("AliEbyEPiDTask"," Task Successfully finished");
  
}
Commit	Line	Data
30409e1c	1	/**************************************************************************
	2	* Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
	3	* *
	4	* Author: Satyajit Jena. *
	5	* Contributors are mentioned in the code where appropriate. *
	6	* *
	7	* Permission to use, copy, modify and distribute this software and its *
	8	* documentation strictly for non-commercial purposes is hereby granted *
	9	* without fee, provided that the above copyright notice appears in all *
	10	* copies and that both the copyright notice and this permission notice *
	11	* appear in the supporting documentation. The authors make no claims *
	12	* about the suitability of this software for any purpose. It is *
	13	* provided "as is" without express or implied warranty. *
	14	**************************************************************************/
	15
	16
	17	//=========================================================================//
	18	// AliEbyE Fluctuaion Analysis for PID //
	19	// Deepika Jena \| drathee@cern.ch //
	20	// Satyajit Jena \| sjena@cern.ch //
	21	//=========================================================================//
	22
	23	#include "TChain.h"
	24	#include "TList.h"
	25	#include "TFile.h"
	26	#include "TTree.h"
	27	#include "TH1D.h"
	28	#include "TH2F.h"
	29	#include "TH3F.h"
	30
	31	#include "THnSparse.h"
	32
	33	#include "AliAnalysisTask.h"
	34	#include "AliAnalysisManager.h"
	35
	36	#include "AliVEvent.h"
	37	#include "AliESDEvent.h"
	38	#include "AliMCEvent.h"
	39	#include "AliAODEvent.h"
	40
	41	#include "AliStack.h"
	42	#include "AliGenEventHeader.h"
	43	#include "AliGenPythiaEventHeader.h"
	44	#include "AliGenHijingEventHeader.h"
	45	#include "AliGenDPMjetEventHeader.h"
	46	#include "TDatabasePDG.h"
	47
	48	#include "AliEbyEPidTaskFastGen.h"
	49
	50	ClassImp(AliEbyEPidTaskFastGen)
	51
	52	//-----------------------------------------------------------------------
	53	AliEbyEPidTaskFastGen::AliEbyEPidTaskFastGen( const char *name )
	54	: AliAnalysisTaskSE( name ),
	55	fThnList(0),
	56	fVxMax(3.),
	57	fVyMax(3.),
	58	fVzMax(15.),
	59	fPtLowerLimit(0.2),
	60	fPtHigherLimit(5.),
	61	fEtaLowerLimit(-1.),
	62	fEtaHigherLimit(1.),
	63	fHistoCorrelationMC(NULL)
	64	{
65
66	DefineOutput(1, TList::Class());
67	}
68
69	AliEbyEPidTaskFastGen::~AliEbyEPidTaskFastGen()
70	{
71	if (fThnList) delete fThnList;
72	}
73
74
75	//---------------------------------------------------------------------------------
76	void AliEbyEPidTaskFastGen::UserCreateOutputObjects() {
77	fThnList = new TList();
78	fThnList->SetOwner(kTRUE);
79	Int_t fgSparseDataBins[kNSparseData] = {50000, 500, 440, 30000, 15000, 15000, 10000, 5000, 5000, 5000, 2500, 2500, 2500, 1250, 1250};
80	Double_t fgSparseDataMin[kNSparseData] = {0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.};
81	Double_t fgSparseDataMax[kNSparseData] = {50000., 100., 440., 30000., 15000., 15000., 10000., 5000., 5000., 5000., 2500., 2500., 2500., 1250., 1250.};
82
83	const Char_t *fgkSparseDataTitle[] = {
84	"Nnumber of Partiles",
85	"Impact Parameters",
86	"N_{p}",
87	"N_{ch}",
88	"N_{+}",
89	"N_{-}",
90	"N_{#pi}",
91	"N_{#pi^{+}}",
92	"N_{#pi^{-}}",
93	"N_{K}",
94	"N_{K^{+}}",
95	"N_{K^{-}}",
96	"N_{pr}",
97	"N_{p}",
98	"N_{#bar{p}}"
99	};
100
101	fHistoCorrelationMC = new THnSparseI("hHistoPR", "", kNSparseData, fgSparseDataBins, fgSparseDataMin, fgSparseDataMax);
102	for (Int_t iaxis = 0; iaxis < kNSparseData; iaxis++)
103	fHistoCorrelationMC->GetAxis(iaxis)->SetTitle(fgkSparseDataTitle[iaxis]);
104	fThnList->Add(fHistoCorrelationMC);
105
106	PostData(1, fThnList);
107	}
108
109
110	//----------------------------------------------------------------------------------
111	void AliEbyEPidTaskFastGen::UserExec( Option_t * ){
112	AliMCEvent* mcEvent = MCEvent();
113	if (!mcEvent) {
114	Printf("ERROR: Could not retrieve MC event");
115	return;
116	}
117
118	AliStack *stack = mcEvent->Stack();
119
120	AliGenEventHeader* genHeader = mcEvent->GenEventHeader();
121	if(!genHeader){
122	printf(" Event generator header not available!!!\n");
123	return;
124	}
125
126	Double_t imp = -1;
127	Double_t np = -1;
128	Double_t nt = -1;
129
130	if(genHeader->InheritsFrom(AliGenHijingEventHeader::Class())){
131	imp = ((AliGenHijingEventHeader*) genHeader)->ImpactParameter();
132	nt = ((AliGenHijingEventHeader*) genHeader)->TargetParticipants();
133	np = ((AliGenHijingEventHeader*) genHeader)->ProjectileParticipants();
134	}
135
136
137	Int_t count[4][2];
138	for (Int_t i = 0; i < 4; i++) {
139	for (Int_t j = 0; j < 2; j++) count[i][j] = 0;
140	}
141
142	Int_t nParticles = stack->GetNtrack();
143	for (Int_t iParticle = 0; iParticle < nParticles; iParticle++) {
144	TParticle* part = stack->Particle(iParticle);
145	if (!part) {
146	Printf(" No Particle Available ");
147	continue;
148	}
149
150	// TDatabasePDG* pdgDB = TDatabasePDG::Instance();
151	// TParticlePDG* pdgPart = pdgDB->GetParticle(part->GetPdgCode());
152
153	Float_t pt = part->Pt();
154	if(pt < fPtLowerLimit \|\| pt > fPtHigherLimit ) continue;
155	Float_t gEta = part->Eta();
156	if(gEta < fEtaLowerLimit \|\| gEta > fEtaHigherLimit ) continue;
157
158
159	// Float_t gCharge = (pdgPart ? pdgPart->Charge() : 0);
160
161	Int_t pid = part->GetPdgCode();
162	if(pid == 211) { count[1][0]++; count[0][0]++; }
163	else if(pid == -211) { count[1][1]++; count[0][1]++; }
164	else if(pid == 321) { count[2][0]++; count[0][0]++; }
165	else if(pid == -321) { count[2][1]++; count[0][1]++; }
166	else if(pid == 2212) { count[3][0]++; count[0][0]++; }
167	else if(pid == -2212) { count[3][1]++; count[0][1]++; }
168	}
169
170	Double_t vsparseMC[kNSparseData];
171	vsparseMC[0] = nParticles;
172	vsparseMC[1] = imp;
173	vsparseMC[2] = np + nt;
174	vsparseMC[3] = count[0][0] + count[0][1];
175	vsparseMC[4] = count[0][0];
176	vsparseMC[5] = count[0][1];
177
178	vsparseMC[6] = count[1][0] + count[1][1];
179	vsparseMC[7] = count[1][0];
180	vsparseMC[8] = count[1][1];
181
182	vsparseMC[9] = count[2][0] + count[2][1];
183	vsparseMC[10] = count[2][0];
184	vsparseMC[11] = count[2][1];
185
186	vsparseMC[12] = count[3][0] + count[3][1];
187	vsparseMC[13] = count[3][0];
188	vsparseMC[14] = count[3][1];
189
190
191	/* Printf(" %10.2f %10.2f %10.2f %10.2f %10.2f %10.2f %10.2f %10.2f %10.2f %10.2f %10.2f %10.2f %10.2f %10.2f %10.2f",
192	vsparseMC[0], vsparseMC[1], vsparseMC[2],
193	vsparseMC[3], vsparseMC[4], vsparseMC[5],
194	vsparseMC[6], vsparseMC[7], vsparseMC[8],
195	vsparseMC[9], vsparseMC[10], vsparseMC[11],
196	vsparseMC[12], vsparseMC[13], vsparseMC[13]);
197	*/
198	fHistoCorrelationMC->Fill(vsparseMC);
199	PostData(1, fThnList);
200
201	}
202
203	void AliEbyEPidTaskFastGen::Terminate( Option_t * ){
204
205	Info("AliEbyEPiDTask"," Task Successfully finished");
206
207	}
208
209