[u/mrichter/AliRoot.git] / PWG4 / AliGammaMCDataReader.cxx


/**************************************************************************
 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
 *                                                                        *
 * Author: The ALICE Off-line Project.                                    *
 * 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.                  *
 **************************************************************************/
/* $Id$ */

/* History of cvs commits:
 *
 * $Log$
 * Revision 1.1.2.1  2007/07/26 10:32:09  schutz
 * new analysis classes in the the new analysis framework
 *
 *
 */

//_________________________________________________________________________
// Class for reading data (Kinematics and ESDs) in order to do prompt gamma correlations
//  Class created from old AliPHOSGammaJet 
//  (see AliRoot versions previous Release 4-09)
//
//*-- Author: Gustavo Conesa (LNF-INFN) 
//////////////////////////////////////////////////////////////////////////////


// --- ROOT system ---
#include "Riostream.h"
#include <TParticle.h>

//---- ANALYSIS system ----
#include "AliGammaMCDataReader.h" 
#include "AliLog.h"
#include "AliESDEvent.h"
#include "AliESDVertex.h"
#include "AliESDCaloCluster.h"
#include "AliStack.h"

ClassImp(AliGammaMCDataReader)

//____________________________________________________________________________
  AliGammaMCDataReader::AliGammaMCDataReader() : 
    AliGammaReader(),  
    fEMCALPID(0),fPHOSPID(0),
    fEMCALPhotonWeight(0.), fEMCALPi0Weight(0.), fPHOSPhotonWeight(0.),  fPHOSPi0Weight(0.) 
{
  //Default Ctor
  
  //Initialize parameters
  fDataType=kMCData;
  InitParameters();
  
}

//____________________________________________________________________________
AliGammaMCDataReader::AliGammaMCDataReader(const AliGammaMCDataReader & g) :   
  AliGammaReader(g),
  fEMCALPID(g.fEMCALPID), 
  fPHOSPID(g.fPHOSPID),
  fEMCALPhotonWeight(g.fEMCALPhotonWeight), 
  fEMCALPi0Weight(g.fEMCALPi0Weight), 
  fPHOSPhotonWeight(g.fPHOSPhotonWeight),
  fPHOSPi0Weight(g.fPHOSPi0Weight)
{
  // cpy ctor
}

//_________________________________________________________________________
AliGammaMCDataReader & AliGammaMCDataReader::operator = (const AliGammaMCDataReader & source)
{
  // assignment operator
  
  if(&source == this) return *this;
  
  fEMCALPID = source.fEMCALPID ;
  fPHOSPID = source.fPHOSPID ;
  fEMCALPhotonWeight = source. fEMCALPhotonWeight ;
  fEMCALPi0Weight = source.fEMCALPi0Weight ;
  fPHOSPhotonWeight = source.fPHOSPhotonWeight ;
  fPHOSPi0Weight = source.fPHOSPi0Weight ;
  
  return *this;
  
}

//____________________________________________________________________________
void AliGammaMCDataReader::CreateParticleList(TObject * data, TObject * kine,
					      TClonesArray * plCTS, 
					      TClonesArray * plEMCAL,  
					      TClonesArray * plPHOS, TClonesArray *){
  
  //Create a list of particles from the ESD. These particles have been measured 
  //by the Central Tracking system (TPC+ITS), PHOS and EMCAL 

  AliESDEvent* esd = (AliESDEvent*) data;
  AliStack* stack = (AliStack*) kine;
  
  Int_t npar  = 0 ;
  Float_t *pid = new Float_t[AliPID::kSPECIESN];  
  AliDebug(3,"Fill particle lists");
  
  //Get vertex for momentum calculation  
  Double_t v[3] ; //vertex ;
  esd->GetVertex()->GetXYZ(v) ; 
  
  //########### PHOS ##############
  
  Int_t begphos = esd->GetFirstPHOSCluster();  
  Int_t endphos = esd->GetFirstPHOSCluster() + 
    esd->GetNumberOfPHOSClusters() ;  
  Int_t indexPH = plPHOS->GetEntries() ;
  AliDebug(3,Form("First PHOS particle %d, last particle %d", begphos,endphos));
  
  for (npar =  begphos; npar <  endphos; npar++) {//////////////PHOS track loop
    AliESDCaloCluster * clus = esd->GetCaloCluster(npar) ; // retrieve track from esd
    if(clus->GetTrackMatched()==-1 ){
      //Create a TParticle to fill the particle list
      TLorentzVector momentum ;
      clus->GetMomentum(momentum, v);
      pid=clus->GetPid();	
      Int_t pdg = 22;

      if(fPHOSPID){
	if( pid[AliPID::kPhoton] > fPHOSPhotonWeight) pdg=22;
	if( pid[AliPID::kPi0] > fPHOSPi0Weight) pdg=111;
	else pdg = 0 ;
      }
      
      TParticle * particle = new TParticle(pdg, 1, -1, -1, -1, -1, 
					   momentum.Px(), momentum.Py(), momentum.Pz(), momentum.E(), v[0], v[1], v[2], 0);

      AliDebug(4,Form("PHOS clusters: pt %f, phi %f, eta %f", particle->Pt(),particle->Phi(),particle->Eta()));
      
      //Look  if parent is prompt photon
      Int_t label = clus->GetLabel();
      if(label>=0){
	TParticle * pmother = stack->Particle(label);
	Int_t imother = pmother->GetFirstMother(); 
	if(imother == 6 || imother == 7)
	  pmother->SetFirstMother(22);
      }
      
      new((*plPHOS)[indexPH++])   TParticle(*particle) ;
    
    }//not charged
  }//cluster loop

  //########### CTS (TPC+ITS) #####################
  Int_t begtpc   = 0 ;  
  Int_t endtpc   = esd->GetNumberOfTracks() ;
  Int_t indexCh  = plCTS->GetEntries() ;
  AliDebug(3,Form("First CTS particle %d, last particle %d", begtpc,endtpc));
  
  for (npar =  begtpc; npar <  endtpc; npar++) {////////////// track loop
    AliESDtrack * track = esd->GetTrack(npar) ; // retrieve track from esd
    
    //We want tracks fitted in the detectors:
    ULong_t status=AliESDtrack::kTPCrefit;
    status|=AliESDtrack::kITSrefit;
    
    //We want tracks whose PID bit is set:
    //     ULong_t status =AliESDtrack::kITSpid;
    //     status|=AliESDtrack::kTPCpid;
    
    if ( (track->GetStatus() & status) == status) {//Check if the bits we want are set
      // Do something with the tracks which were successfully
      // re-fitted 
      Double_t en = 0; //track ->GetTPCsignal() ;
      Double_t mom[3];
      track->GetPxPyPz(mom) ;
      Double_t px = mom[0];
      Double_t py = mom[1];
      Double_t pz = mom[2]; //Check with TPC people if this is correct.
      Int_t pdg = 11; //Give any charged PDG code, in this case electron.
      //I just want to tag the particle as charged
       TParticle * particle = new TParticle(pdg, 1, -1, -1, -1, -1, 
					    px, py, pz, en, v[0], v[1], v[2], 0);
  
      //TParticle * particle = new TParticle() ;
      //particle->SetMomentum(px,py,pz,en) ;
 
      new((*plCTS)[indexCh++])       TParticle(*particle) ;    
    }
  }
  
  //################ EMCAL ##############
  
  Int_t indexEM  = plEMCAL->GetEntries() ; 
  Int_t begem = esd->GetFirstEMCALCluster();  
  Int_t endem = esd->GetFirstEMCALCluster() + 
    esd->GetNumberOfEMCALClusters() ;  
  
  AliDebug(3,Form("First EMCAL particle %d, last particle %d",begem,endem));
  
  for (npar =  begem; npar <  endem; npar++) {//////////////EMCAL track loop
    AliESDCaloCluster * clus = esd->GetCaloCluster(npar) ; // retrieve track from esd
    Int_t clustertype= clus->GetClusterType();
    if(clustertype == AliESDCaloCluster::kClusterv1 && clus->GetTrackMatched()==-1 ){
      
      TLorentzVector momentum ;
      clus->GetMomentum(momentum, v);            
      pid=clus->GetPid();	
      Int_t pdg = 22;

      if(fEMCALPID){
	if( pid[AliPID::kPhoton] > fEMCALPhotonWeight) pdg=22;
	else if( pid[AliPID::kPi0] > fEMCALPi0Weight) pdg=111;
	else pdg = 0;
      }
      
      TParticle * particle = new TParticle(pdg, 1, -1, -1, -1, -1, 
					   momentum.Px(), momentum.Py(), momentum.Pz(), momentum.E(), v[0], v[1], v[2], 0);
      AliDebug(4,Form("EMCAL clusters: pt %f, phi %f, eta %f", particle->Pt(),particle->Phi(),particle->Eta()));

      //Look  if parent is prompt photon
      Int_t label = clus->GetLabel();
      if(label>=0){
	TParticle * pmother = stack->Particle(label);
	Int_t imother = pmother->GetFirstMother(); 
	if(imother == 6 || imother == 7)
	  pmother->SetFirstMother(22);
      }
      
      new((*plEMCAL)[indexEM++])   TParticle(*particle) ;
      
    }//not charged, not pseudocluster
  }//Cluster loop
  
  AliDebug(3,"Particle lists filled");
  
}

  //____________________________________________________________________________
void AliGammaMCDataReader::InitParameters()
{
 
  //Initialize the parameters of the analysis.

  //Fill particle lists when PID is ok
  fEMCALPID = kFALSE;
  fPHOSPID = kFALSE;
  fEMCALPhotonWeight = 0.5 ;
  fEMCALPi0Weight = 0.5 ;
  fPHOSPhotonWeight = 0.8 ;
  fPHOSPi0Weight = 0.5 ;

}


void AliGammaMCDataReader::Print(const Option_t * opt) const
{

  //Print some relevant parameters set for the analysis
  if(! opt)
    return;
  
  Info("Print", "%s %s", GetName(), GetTitle() ) ;
  printf("PHOS PID on?               =     %d\n",  fPHOSPID) ; 
  printf("EMCAL PID  on?         =     %d\n",  fEMCALPID) ;
  printf("PHOS PID weight , photon  %f, pi0 %f\n\n",  fPHOSPhotonWeight,  fPHOSPi0Weight) ; 
  printf("EMCAL PID weight, photon %f, pi0 %f\n",   fEMCALPhotonWeight,  fEMCALPi0Weight) ; 

}
Commit	Line	Data
bdcfac30	1
	2	/**************************************************************************
	3	* Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
	4	* *
	5	* Author: The ALICE Off-line Project. *
	6	* Contributors are mentioned in the code where appropriate. *
	7	* *
	8	* Permission to use, copy, modify and distribute this software and its *
	9	* documentation strictly for non-commercial purposes is hereby granted *
	10	* without fee, provided that the above copyright notice appears in all *
	11	* copies and that both the copyright notice and this permission notice *
	12	* appear in the supporting documentation. The authors make no claims *
	13	* about the suitability of this software for any purpose. It is *
	14	* provided "as is" without express or implied warranty. *
	15	**************************************************************************/
	16	/* $Id$ */
	17
	18	/* History of cvs commits:
	19	*
	20	* $Log$
	21	* Revision 1.1.2.1 2007/07/26 10:32:09 schutz
	22	* new analysis classes in the the new analysis framework
	23	*
	24	*
	25	*/
	26
	27	//_________________________________________________________________________
	28	// Class for reading data (Kinematics and ESDs) in order to do prompt gamma correlations
	29	// Class created from old AliPHOSGammaJet
	30	// (see AliRoot versions previous Release 4-09)
	31	//
	32	//*-- Author: Gustavo Conesa (LNF-INFN)
	33	//////////////////////////////////////////////////////////////////////////////
	34
	35
	36	// --- ROOT system ---
	37	#include "Riostream.h"
	38	#include <TParticle.h>
	39
	40	//---- ANALYSIS system ----
	41	#include "AliGammaMCDataReader.h"
	42	#include "AliLog.h"
	43	#include "AliESDEvent.h"
	44	#include "AliESDVertex.h"
	45	#include "AliESDCaloCluster.h"
	46	#include "AliStack.h"
	47
	48	ClassImp(AliGammaMCDataReader)
	49
	50	//____________________________________________________________________________
	51	AliGammaMCDataReader::AliGammaMCDataReader() :
	52	AliGammaReader(),
	53	fEMCALPID(0),fPHOSPID(0),
	54	fEMCALPhotonWeight(0.), fEMCALPi0Weight(0.), fPHOSPhotonWeight(0.), fPHOSPi0Weight(0.)
	55	{
	56	//Default Ctor
	57
	58	//Initialize parameters
	59	fDataType=kMCData;
	60	InitParameters();
	61
	62	}
	63
	64	//____________________________________________________________________________
65	AliGammaMCDataReader::AliGammaMCDataReader(const AliGammaMCDataReader & g) :
66	AliGammaReader(g),
67	fEMCALPID(g.fEMCALPID),
68	fPHOSPID(g.fPHOSPID),
69	fEMCALPhotonWeight(g.fEMCALPhotonWeight),
70	fEMCALPi0Weight(g.fEMCALPi0Weight),
71	fPHOSPhotonWeight(g.fPHOSPhotonWeight),
72	fPHOSPi0Weight(g.fPHOSPi0Weight)
73	{
74	// cpy ctor
75	}
76
77	//_________________________________________________________________________
78	AliGammaMCDataReader & AliGammaMCDataReader::operator = (const AliGammaMCDataReader & source)
79	{
80	// assignment operator
81
82	if(&source == this) return *this;
83
84	fEMCALPID = source.fEMCALPID ;
85	fPHOSPID = source.fPHOSPID ;
86	fEMCALPhotonWeight = source. fEMCALPhotonWeight ;
87	fEMCALPi0Weight = source.fEMCALPi0Weight ;
88	fPHOSPhotonWeight = source.fPHOSPhotonWeight ;
89	fPHOSPi0Weight = source.fPHOSPi0Weight ;
90
91	return *this;
92
93	}
94
95	//____________________________________________________________________________
96	void AliGammaMCDataReader::CreateParticleList(TObject * data, TObject * kine,
97	TClonesArray * plCTS,
98	TClonesArray * plEMCAL,
99	TClonesArray * plPHOS, TClonesArray *){
100
101	//Create a list of particles from the ESD. These particles have been measured
102	//by the Central Tracking system (TPC+ITS), PHOS and EMCAL
103
104	AliESDEvent* esd = (AliESDEvent*) data;
105	AliStack* stack = (AliStack*) kine;
106
107	Int_t npar = 0 ;
108	Float_t *pid = new Float_t[AliPID::kSPECIESN];
109	AliDebug(3,"Fill particle lists");
110
111	//Get vertex for momentum calculation
112	Double_t v[3] ; //vertex ;
113	esd->GetVertex()->GetXYZ(v) ;
114
115	//########### PHOS ##############
116
117	Int_t begphos = esd->GetFirstPHOSCluster();
118	Int_t endphos = esd->GetFirstPHOSCluster() +
119	esd->GetNumberOfPHOSClusters() ;
120	Int_t indexPH = plPHOS->GetEntries() ;
121	AliDebug(3,Form("First PHOS particle %d, last particle %d", begphos,endphos));
122
123	for (npar = begphos; npar < endphos; npar++) {//////////////PHOS track loop
124	AliESDCaloCluster * clus = esd->GetCaloCluster(npar) ; // retrieve track from esd
125	if(clus->GetTrackMatched()==-1 ){
126	//Create a TParticle to fill the particle list
127	TLorentzVector momentum ;
128	clus->GetMomentum(momentum, v);
129	pid=clus->GetPid();
130	Int_t pdg = 22;
131
132	if(fPHOSPID){
133	if( pid[AliPID::kPhoton] > fPHOSPhotonWeight) pdg=22;
134	if( pid[AliPID::kPi0] > fPHOSPi0Weight) pdg=111;
135	else pdg = 0 ;
136	}
137
138	TParticle * particle = new TParticle(pdg, 1, -1, -1, -1, -1,
139	momentum.Px(), momentum.Py(), momentum.Pz(), momentum.E(), v[0], v[1], v[2], 0);
140
141	AliDebug(4,Form("PHOS clusters: pt %f, phi %f, eta %f", particle->Pt(),particle->Phi(),particle->Eta()));
142
143	//Look if parent is prompt photon
144	Int_t label = clus->GetLabel();
145	if(label>=0){
146	TParticle * pmother = stack->Particle(label);
147	Int_t imother = pmother->GetFirstMother();
148	if(imother == 6 \|\| imother == 7)
149	pmother->SetFirstMother(22);
150	}
151
152	new((plPHOS)[indexPH++]) TParticle(particle) ;
153
154	}//not charged
155	}//cluster loop
156
157	//########### CTS (TPC+ITS) #####################
158	Int_t begtpc = 0 ;
159	Int_t endtpc = esd->GetNumberOfTracks() ;
160	Int_t indexCh = plCTS->GetEntries() ;
161	AliDebug(3,Form("First CTS particle %d, last particle %d", begtpc,endtpc));
162
163	for (npar = begtpc; npar < endtpc; npar++) {////////////// track loop
164	AliESDtrack * track = esd->GetTrack(npar) ; // retrieve track from esd
165
166	//We want tracks fitted in the detectors:
167	ULong_t status=AliESDtrack::kTPCrefit;
168	status\|=AliESDtrack::kITSrefit;
169
170	//We want tracks whose PID bit is set:
171	// ULong_t status =AliESDtrack::kITSpid;
172	// status\|=AliESDtrack::kTPCpid;
173
174	if ( (track->GetStatus() & status) == status) {//Check if the bits we want are set
175	// Do something with the tracks which were successfully
176	// re-fitted
177	Double_t en = 0; //track ->GetTPCsignal() ;
178	Double_t mom[3];
179	track->GetPxPyPz(mom) ;
180	Double_t px = mom[0];
181	Double_t py = mom[1];
182	Double_t pz = mom[2]; //Check with TPC people if this is correct.
183	Int_t pdg = 11; //Give any charged PDG code, in this case electron.
184	//I just want to tag the particle as charged
185	TParticle * particle = new TParticle(pdg, 1, -1, -1, -1, -1,
186	px, py, pz, en, v[0], v[1], v[2], 0);
187
188	//TParticle * particle = new TParticle() ;
189	//particle->SetMomentum(px,py,pz,en) ;
190
191	new((plCTS)[indexCh++]) TParticle(particle) ;
192	}
193	}
194
195	//################ EMCAL ##############
196
197	Int_t indexEM = plEMCAL->GetEntries() ;
198	Int_t begem = esd->GetFirstEMCALCluster();
199	Int_t endem = esd->GetFirstEMCALCluster() +
200	esd->GetNumberOfEMCALClusters() ;
201
202	AliDebug(3,Form("First EMCAL particle %d, last particle %d",begem,endem));
203
204	for (npar = begem; npar < endem; npar++) {//////////////EMCAL track loop
205	AliESDCaloCluster * clus = esd->GetCaloCluster(npar) ; // retrieve track from esd
206	Int_t clustertype= clus->GetClusterType();
207	if(clustertype == AliESDCaloCluster::kClusterv1 && clus->GetTrackMatched()==-1 ){
208
209	TLorentzVector momentum ;
210	clus->GetMomentum(momentum, v);
211	pid=clus->GetPid();
212	Int_t pdg = 22;
213
214	if(fEMCALPID){
215	if( pid[AliPID::kPhoton] > fEMCALPhotonWeight) pdg=22;
216	else if( pid[AliPID::kPi0] > fEMCALPi0Weight) pdg=111;
217	else pdg = 0;
218	}
219
220	TParticle * particle = new TParticle(pdg, 1, -1, -1, -1, -1,
221	momentum.Px(), momentum.Py(), momentum.Pz(), momentum.E(), v[0], v[1], v[2], 0);
222	AliDebug(4,Form("EMCAL clusters: pt %f, phi %f, eta %f", particle->Pt(),particle->Phi(),particle->Eta()));
223
224	//Look if parent is prompt photon
225	Int_t label = clus->GetLabel();
226	if(label>=0){
227	TParticle * pmother = stack->Particle(label);
228	Int_t imother = pmother->GetFirstMother();
229	if(imother == 6 \|\| imother == 7)
230	pmother->SetFirstMother(22);
231	}
232
233	new((plEMCAL)[indexEM++]) TParticle(particle) ;
234
235	}//not charged, not pseudocluster
236	}//Cluster loop
237
238	AliDebug(3,"Particle lists filled");
239
240	}
241
242	//____________________________________________________________________________
243	void AliGammaMCDataReader::InitParameters()
244	{
245
246	//Initialize the parameters of the analysis.
247
248	//Fill particle lists when PID is ok
249	fEMCALPID = kFALSE;
250	fPHOSPID = kFALSE;
251	fEMCALPhotonWeight = 0.5 ;
252	fEMCALPi0Weight = 0.5 ;
253	fPHOSPhotonWeight = 0.8 ;
254	fPHOSPi0Weight = 0.5 ;
255
256	}
257
258
259	void AliGammaMCDataReader::Print(const Option_t * opt) const
260	{
261
262	//Print some relevant parameters set for the analysis
263	if(! opt)
264	return;
265
266	Info("Print", "%s %s", GetName(), GetTitle() ) ;
267	printf("PHOS PID on? = %d\n", fPHOSPID) ;
268	printf("EMCAL PID on? = %d\n", fEMCALPID) ;
269	printf("PHOS PID weight , photon %f, pi0 %f\n\n", fPHOSPhotonWeight, fPHOSPi0Weight) ;
270	printf("EMCAL PID weight, photon %f, pi0 %f\n", fEMCALPhotonWeight, fEMCALPi0Weight) ;
271
272	}