[u/mrichter/AliRoot.git] / PWG4 / AliGammaDataReader.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 (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 "AliGammaDataReader.h" 
#include "AliLog.h"
#include "AliESDEvent.h"
#include "AliESDVertex.h"
#include "AliESDCaloCluster.h"

ClassImp(AliGammaDataReader)

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

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

//_________________________________________________________________________
AliGammaDataReader & AliGammaDataReader::operator = (const AliGammaDataReader & 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 AliGammaDataReader::CreateParticleList(TObject * data, TObject *,
					    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;

  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 = 0;  
  Int_t endphos = esd->GetNumberOfCaloClusters() ;  
  Int_t indexPH = plPHOS->GetEntries() ;
  Int_t begem = 0;
  Int_t endem = endphos;

  AliDebug(3,Form("First Calorimeter 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

    //Check if it is PHOS cluster
    if(!clus->IsEMCAL())
      begem=npar+1;
    else
      continue;

    AliDebug(4,Form("PHOS clusters: E %f, match %d", clus->E(),clus->GetTrackMatched()));
    if(clus->GetTrackMatched()==-1){
      //Create a TParticle to fill the particle list
      TLorentzVector momentum ;
      clus->GetMomentum(momentum, v);      
      Double_t phi = momentum.Phi();
      if(phi<0) phi+=TMath::TwoPi() ;
      if(momentum.Pt() > fNeutralPtCut &&  TMath::Abs(momentum.Eta()) < fPHOSEtaCut &&
	 phi > fPhiPHOSCut[0] && phi < fPhiPHOSCut[1] ) {
      	
	pid=clus->GetPid();	
	Int_t pdg = 22;
	
	if(fPHOSPID){
	  AliDebug(4, Form("PID: photon %f, pi0 %f, electron %f",pid[AliPID::kPhoton],pid[AliPID::kPi0],pid[AliPID::kElectron]));
	  if( pid[AliPID::kPhoton] > fPHOSPhotonWeight) pdg=22;
	  else 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(3,Form("PHOS added: pt %f, phi %f, eta %f", particle->Pt(),particle->Phi(),particle->Eta()));
	
	new((*plPHOS)[indexPH++])   TParticle(*particle) ;
	
      }//pt, eta, phi cut
      else 	AliDebug(4,"Particle not added");
    }//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) ;
       if(particle->Pt() > fChargedPtCut && TMath::Abs(particle->Eta())<fCTSEtaCut)
	 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();
    AliDebug(4,Form("EMCAL clusters: E %f, match %d", clus->E(),clus->GetTrackMatched()));

    if(clustertype == AliESDCaloCluster::kClusterv1 && clus->GetTrackMatched()==-1 ){
      
      TLorentzVector momentum ;
      clus->GetMomentum(momentum, v); 
      Double_t phi = momentum.Phi();
      if(phi<0) phi+=TMath::TwoPi() ;
      if(momentum.Pt() > fNeutralPtCut &&  TMath::Abs(momentum.Eta()) < fEMCALEtaCut &&
	phi > fPhiEMCALCut[0] && phi < fPhiEMCALCut[1] ) {
      
	pid=clus->GetPid();	
	Int_t pdg = 22;
	if(fEMCALPID){
	  AliDebug(4, Form("PID: photon %f, pi0 %f, electron %f",pid[AliPID::kPhoton],pid[AliPID::kPi0],pid[AliPID::kElectron]));
	  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(3,Form("EMCAL added: pt %f, phi %f, eta %f", particle->Pt(),particle->Phi(),particle->Eta()));
	
	new((*plEMCAL)[indexEM++])   TParticle(*particle) ;
  
      }//pt, phi, eta cut
      else 	AliDebug(4,"Particle not added");
    }//not charged, not pseudocluster
  }//Cluster loop
  
  AliDebug(3,"Particle lists filled");
  
}

  //____________________________________________________________________________
void AliGammaDataReader::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 ;

}


void AliGammaDataReader::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\n",  fPHOSPhotonWeight) ; 
  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 (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 "AliGammaDataReader.h"
	42	#include "AliLog.h"
	43	#include "AliESDEvent.h"
	44	#include "AliESDVertex.h"
	45	#include "AliESDCaloCluster.h"
	46
	47	ClassImp(AliGammaDataReader)
	48
	49	//____________________________________________________________________________
	50	AliGammaDataReader::AliGammaDataReader() :
	51	AliGammaReader(),
	52	fEMCALPID(0),fPHOSPID(0),
	53	fEMCALPhotonWeight(0.), fEMCALPi0Weight(0.), fPHOSPhotonWeight(0.), fPHOSPi0Weight(0.)
	54	{
	55	//Default Ctor
	56
	57	//Initialize parameters
	58	fDataType=kData;
	59	InitParameters();
	60
	61	}
	62
	63	//____________________________________________________________________________
	64	AliGammaDataReader::AliGammaDataReader(const AliGammaDataReader & g) :
65	AliGammaReader(g),
66	fEMCALPID(g.fEMCALPID),
67	fPHOSPID(g.fPHOSPID),
68	fEMCALPhotonWeight(g.fEMCALPhotonWeight),
69	fEMCALPi0Weight(g.fEMCALPi0Weight),
70	fPHOSPhotonWeight(g.fPHOSPhotonWeight),
71	fPHOSPi0Weight(g.fPHOSPi0Weight)
72	{
73	// cpy ctor
74	}
75
76	//_________________________________________________________________________
77	AliGammaDataReader & AliGammaDataReader::operator = (const AliGammaDataReader & source)
78	{
79	// assignment operator
80
81	if(&source == this) return *this;
82
83	fEMCALPID = source.fEMCALPID ;
84	fPHOSPID = source.fPHOSPID ;
85	fEMCALPhotonWeight = source. fEMCALPhotonWeight ;
86	fEMCALPi0Weight = source.fEMCALPi0Weight ;
87	fPHOSPhotonWeight = source.fPHOSPhotonWeight ;
88	fPHOSPi0Weight = source.fPHOSPi0Weight ;
89
90	return *this;
91
92	}
93
94	//____________________________________________________________________________
95	void AliGammaDataReader::CreateParticleList(TObject * data, TObject *,
96	TClonesArray * plCTS,
97	TClonesArray * plEMCAL,
98	TClonesArray * plPHOS, TClonesArray*){
99
100	//Create a list of particles from the ESD. These particles have been measured
101	//by the Central Tracking system (TPC+ITS), PHOS and EMCAL
102
103	AliESDEvent* esd = (AliESDEvent*) data;
104
105	Int_t npar = 0 ;
106	Float_t *pid = new Float_t[AliPID::kSPECIESN];
107	AliDebug(3,"Fill particle lists");
108
109	//Get vertex for momentum calculation
110	Double_t v[3] ; //vertex ;
111	esd->GetVertex()->GetXYZ(v) ;
112
113	//########### PHOS ##############
114
115	Int_t begphos = 0;
116	Int_t endphos = esd->GetNumberOfCaloClusters() ;
117	Int_t indexPH = plPHOS->GetEntries() ;
118	Int_t begem = 0;
119	Int_t endem = endphos;
120
121	AliDebug(3,Form("First Calorimeter 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
126	//Check if it is PHOS cluster
127	if(!clus->IsEMCAL())
128	begem=npar+1;
129	else
130	continue;
131
132	AliDebug(4,Form("PHOS clusters: E %f, match %d", clus->E(),clus->GetTrackMatched()));
133	if(clus->GetTrackMatched()==-1){
134	//Create a TParticle to fill the particle list
135	TLorentzVector momentum ;
136	clus->GetMomentum(momentum, v);
137	Double_t phi = momentum.Phi();
138	if(phi<0) phi+=TMath::TwoPi() ;
139	if(momentum.Pt() > fNeutralPtCut && TMath::Abs(momentum.Eta()) < fPHOSEtaCut &&
140	phi > fPhiPHOSCut[0] && phi < fPhiPHOSCut[1] ) {
141
142	pid=clus->GetPid();
143	Int_t pdg = 22;
144
145	if(fPHOSPID){
146	AliDebug(4, Form("PID: photon %f, pi0 %f, electron %f",pid[AliPID::kPhoton],pid[AliPID::kPi0],pid[AliPID::kElectron]));
147	if( pid[AliPID::kPhoton] > fPHOSPhotonWeight) pdg=22;
148	else if( pid[AliPID::kPi0] > fPHOSPi0Weight) pdg=111;
149	else pdg = 0;
150	}
151
152	TParticle * particle = new TParticle(pdg, 1, -1, -1, -1, -1,
153	momentum.Px(), momentum.Py(), momentum.Pz(), momentum.E(), v[0], v[1], v[2], 0);
154
155	AliDebug(3,Form("PHOS added: pt %f, phi %f, eta %f", particle->Pt(),particle->Phi(),particle->Eta()));
156
157	new((plPHOS)[indexPH++]) TParticle(particle) ;
158
159	}//pt, eta, phi cut
160	else AliDebug(4,"Particle not added");
161	}//not charged
162	}//cluster loop
163
164	//########### CTS (TPC+ITS) #####################
165	Int_t begtpc = 0 ;
166	Int_t endtpc = esd->GetNumberOfTracks() ;
167	Int_t indexCh = plCTS->GetEntries() ;
168	AliDebug(3,Form("First CTS particle %d, last particle %d", begtpc,endtpc));
169
170	for (npar = begtpc; npar < endtpc; npar++) {////////////// track loop
171	AliESDtrack * track = esd->GetTrack(npar) ; // retrieve track from esd
172
173	//We want tracks fitted in the detectors:
174	ULong_t status=AliESDtrack::kTPCrefit;
175	status\|=AliESDtrack::kITSrefit;
176
177	//We want tracks whose PID bit is set:
178	// ULong_t status =AliESDtrack::kITSpid;
179	// status\|=AliESDtrack::kTPCpid;
180
181	if ( (track->GetStatus() & status) == status) {//Check if the bits we want are set
182	// Do something with the tracks which were successfully
183	// re-fitted
184	Double_t en = 0; //track ->GetTPCsignal() ;
185	Double_t mom[3];
186	track->GetPxPyPz(mom) ;
187	Double_t px = mom[0];
188	Double_t py = mom[1];
189	Double_t pz = mom[2]; //Check with TPC people if this is correct.
190	Int_t pdg = 11; //Give any charged PDG code, in this case electron.
191	//I just want to tag the particle as charged
192	TParticle * particle = new TParticle(pdg, 1, -1, -1, -1, -1,
193	px, py, pz, en, v[0], v[1], v[2], 0);
194
195	//TParticle * particle = new TParticle() ;
196	//particle->SetMomentum(px,py,pz,en) ;
197	if(particle->Pt() > fChargedPtCut && TMath::Abs(particle->Eta())<fCTSEtaCut)
198	new((plCTS)[indexCh++]) TParticle(particle) ;
199	}
200	}
201
202	//################ EMCAL ##############
203
204	Int_t indexEM = plEMCAL->GetEntries() ;
205	//Int_t begem = esd->GetFirstEMCALCluster();
206	//Int_t endem = esd->GetFirstEMCALCluster() +
207	//esd->GetNumberOfEMCALClusters() ;
208
209	//AliDebug(3,Form("First EMCAL particle %d, last particle %d",begem,endem));
210
211	for (npar = begem; npar < endem; npar++) {//////////////EMCAL track loop
212	AliESDCaloCluster * clus = esd->GetCaloCluster(npar) ; // retrieve track from esd
213	Int_t clustertype= clus->GetClusterType();
214	AliDebug(4,Form("EMCAL clusters: E %f, match %d", clus->E(),clus->GetTrackMatched()));
215
216	if(clustertype == AliESDCaloCluster::kClusterv1 && clus->GetTrackMatched()==-1 ){
217
218	TLorentzVector momentum ;
219	clus->GetMomentum(momentum, v);
220	Double_t phi = momentum.Phi();
221	if(phi<0) phi+=TMath::TwoPi() ;
222	if(momentum.Pt() > fNeutralPtCut && TMath::Abs(momentum.Eta()) < fEMCALEtaCut &&
223	phi > fPhiEMCALCut[0] && phi < fPhiEMCALCut[1] ) {
224
225	pid=clus->GetPid();
226	Int_t pdg = 22;
227	if(fEMCALPID){
228	AliDebug(4, Form("PID: photon %f, pi0 %f, electron %f",pid[AliPID::kPhoton],pid[AliPID::kPi0],pid[AliPID::kElectron]));
229	if( pid[AliPID::kPhoton] > fEMCALPhotonWeight) pdg=22;
230	else if( pid[AliPID::kPi0] > fEMCALPi0Weight) pdg=111;
231	else pdg = 0 ;
232	}
233
234
235	TParticle * particle = new TParticle(pdg, 1, -1, -1, -1, -1,
236	momentum.Px(), momentum.Py(), momentum.Pz(), momentum.E(), v[0], v[1], v[2], 0);
237	AliDebug(3,Form("EMCAL added: pt %f, phi %f, eta %f", particle->Pt(),particle->Phi(),particle->Eta()));
238
239	new((plEMCAL)[indexEM++]) TParticle(particle) ;
240
241	}//pt, phi, eta cut
242	else AliDebug(4,"Particle not added");
243	}//not charged, not pseudocluster
244	}//Cluster loop
245
246	AliDebug(3,"Particle lists filled");
247
248	}
249
250	//____________________________________________________________________________
251	void AliGammaDataReader::InitParameters()
252	{
253
254	//Initialize the parameters of the analysis.
255
256	//Fill particle lists when PID is ok
257	fEMCALPID = kFALSE;
258	fPHOSPID = kFALSE;
259	fEMCALPhotonWeight = 0.5 ;
260	fEMCALPi0Weight = 0.5 ;
261	fPHOSPhotonWeight = 0.8 ;
262
263	}
264
265
266	void AliGammaDataReader::Print(const Option_t * opt) const
267	{
268
269	//Print some relevant parameters set for the analysis
270	if(! opt)
271	return;
272
273	Info("Print", "%s %s", GetName(), GetTitle() ) ;
274	printf("PHOS PID on? = %d\n", fPHOSPID) ;
275	printf("EMCAL PID on? = %d\n", fEMCALPID) ;
276	printf("PHOS PID weight , photon %f\n", fPHOSPhotonWeight) ;
277	printf("EMCAL PID weight, photon %f, pi0 %f\n", fEMCALPhotonWeight, fEMCALPi0Weight) ;
278
279
280	}