1 /**************************************************************************
2 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
4 * Author: The ALICE Off-line Project. *
5 * Contributors are mentioned in the code where appropriate. *
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 **************************************************************************/
17 /* History of cvs commits:
20 * Revision 1.14 2006/08/28 10:01:56 kharlov
21 * Effective C++ warnings fixed (Timur Pocheptsov)
23 * Revision 1.13 2006/04/26 07:32:37 hristov
24 * Coding conventions, clean-up and related changes
26 * Revision 1.12 2006/03/10 13:23:36 hristov
27 * Using AliESDCaloCluster instead of AliESDtrack
29 * Revision 1.11 2006/01/31 20:30:52 hristov
32 * Revision 1.10 2006/01/23 18:04:08 hristov
33 * Removing meaningless const
35 * Revision 1.9 2006/01/12 16:23:26 schutz
36 * ESD is properly read with methods of macros/AliReadESD.C copied in it
38 * Revision 1.8 2005/12/20 07:08:32 schutz
39 * corrected error in call AliReadESD
41 * Revision 1.6 2005/05/28 14:19:04 schutz
42 * Compilation warnings fixed by T.P.
46 //_________________________________________________________________________
47 // Class for the analysis of gamma-jet correlations
48 // Basically it seaches for a prompt photon in the PHOS acceptance,
49 // if so we construct a jet around the highest pt particle in the opposite
50 // side in azimuth, inside the TPC and EMCAL acceptances. First the leading
51 // particle and then the jet have to fullfill several conditions
52 // (energy, direction ..) to be accepted. Then the fragmentation function
53 // of this jet is constructed
55 //*-- Author: Gustavo Conesa & Yves Schutz (IFIC, CERN)
56 //////////////////////////////////////////////////////////////////////////////
59 // --- ROOT system ---
62 #include <TParticle.h>
64 #include <TPaveLabel.h>
68 #include "AliPHOSGammaJet.h"
69 #include "AliPHOSGetter.h"
70 #include "AliPHOSGeometry.h"
71 #include "AliPHOSFastGlobalReconstruction.h"
73 #include "AliESDtrack.h"
74 #include "AliESDCaloCluster.h"
75 #include "Riostream.h"
78 ClassImp(AliPHOSGammaJet)
80 //____________________________________________________________________________
81 AliPHOSGammaJet::AliPHOSGammaJet() :
82 TTask(), fAnyConeOrPt(0), fOptionGJ(""),
83 fOutputFile(new TFile(gDirectory->GetName())),
84 fOutputFileName(gDirectory->GetName()),
85 fInputFileName(gDirectory->GetName()),
86 fHIJINGFileName(gDirectory->GetName()),
87 fHIJING(0), fESDdata(0), fEtaCut(0.),
88 fOnlyCharged(0), fPhiMaxCut(0.),
89 fPhiMinCut(0.), fPtCut(0.),
90 fNeutralPtCut(0.), fChargedPtCut(0.),
91 fInvMassMaxCut(0.), fInvMassMinCut(0.),
92 fMinDistance(0.), fRatioMaxCut(0.), fRatioMinCut(0.),
93 fTPCCutsLikeEMCAL(0), fDirName(""), fESDTree(""),
94 fPattern(""), fJetTPCRatioMaxCut(0.),
95 fJetTPCRatioMinCut(0.), fJetRatioMaxCut(0.),
96 fJetRatioMinCut(0.), fNEvent(0), fNCone(0),
97 fNPt(0), fCone(0), fPtThreshold(0),
98 fPtJetSelectionCut(0.0),
99 fListHistos(new TObjArray(100)),
100 fFastRec(0), fOptFast(0),
101 fRan(0), fResPara1(0.), fResPara2(0.), fResPara3(0.),
102 fPosParaA(0.), fPosParaB(0.), fAngleMaxParam(), fSelect(0)
105 fAngleMaxParam.Set(4) ;
106 fAngleMaxParam.Reset(0.);
108 for(Int_t i = 0; i<10; i++){
112 fNamePtThres[i] = "" ;
123 fJetSigma1[i] = 0.0 ;
124 fJetSigma2[i] = 0.0 ;
125 fPhiEMCALCut[i] = 0.0 ;
130 TList * list = gDirectory->GetListOfKeys() ;
134 for (index = 0 ; index < list->GetSize()-1 ; index++) {
135 //-1 to avoid GammaJet Task
136 h = dynamic_cast<TH2F*>(gDirectory->Get(list->At(index)->GetName())) ;
137 fListHistos->Add(h) ;
142 //____________________________________________________________________________
143 AliPHOSGammaJet::AliPHOSGammaJet(const TString inputfilename) :
144 TTask("GammaJet","Analysis of gamma-jet correlations"),
145 fAnyConeOrPt(0), fOptionGJ(),
150 fHIJING(0), fESDdata(0), fEtaCut(0.),
151 fOnlyCharged(0), fPhiMaxCut(0.),
152 fPhiMinCut(0.), fPtCut(0.),
153 fNeutralPtCut(0.), fChargedPtCut(0.),
154 fInvMassMaxCut(0.), fInvMassMinCut(0.),
155 fMinDistance(0.), fRatioMaxCut(0.), fRatioMinCut(0.),
156 fTPCCutsLikeEMCAL(0), fDirName(), fESDTree(),
157 fPattern(), fJetTPCRatioMaxCut(0.),
158 fJetTPCRatioMinCut(0.), fJetRatioMaxCut(0.),
159 fJetRatioMinCut(0.), fNEvent(0), fNCone(0),
160 fNPt(0), fCone(0), fPtThreshold(0),
161 fPtJetSelectionCut(0.0),
163 fFastRec(0), fOptFast(0),
164 fRan(0), fResPara1(0.), fResPara2(0.), fResPara3(0.),
165 fPosParaA(0.), fPosParaB(0.), fAngleMaxParam(), fSelect(0)
169 fInputFileName = inputfilename;
170 fFastRec = new AliPHOSFastGlobalReconstruction(fInputFileName);
171 AliPHOSGetter * gime = AliPHOSGetter::Instance(fInputFileName) ;
172 fNEvent = gime->MaxEvent();
176 //____________________________________________________________________________
177 AliPHOSGammaJet::AliPHOSGammaJet(const AliPHOSGammaJet & gj) :
179 fAnyConeOrPt(gj.fAnyConeOrPt), fOptionGJ(gj.fOptionGJ),
180 fOutputFile(gj.fOutputFile),
181 fOutputFileName(gj.fOutputFileName),
182 fInputFileName(gj.fInputFileName),
183 fHIJINGFileName(gj.fHIJINGFileName),
184 fHIJING(gj.fHIJING), fESDdata(gj.fESDdata), fEtaCut(gj.fEtaCut),
185 fOnlyCharged(gj.fOnlyCharged), fPhiMaxCut(gj.fPhiMaxCut),
186 fPhiMinCut(gj.fPhiMinCut), fPtCut(gj.fPtCut),
187 fNeutralPtCut(gj.fNeutralPtCut), fChargedPtCut(gj.fChargedPtCut),
188 fInvMassMaxCut(gj.fInvMassMaxCut), fInvMassMinCut(gj.fInvMassMinCut),
189 fMinDistance(gj.fMinDistance), fRatioMaxCut(gj.fRatioMaxCut),
190 fRatioMinCut(gj.fRatioMinCut), fTPCCutsLikeEMCAL(gj.fTPCCutsLikeEMCAL),
191 fDirName(gj.fDirName), fESDTree(gj.fESDTree),
192 fPattern(gj.fPattern), fJetTPCRatioMaxCut(gj.fJetTPCRatioMaxCut),
193 fJetTPCRatioMinCut(gj.fJetTPCRatioMinCut), fJetRatioMaxCut(gj.fJetRatioMaxCut),
194 fJetRatioMinCut(gj.fJetRatioMinCut), fNEvent(gj.fNEvent), fNCone(gj.fNCone),
195 fNPt(gj.fNPt), fCone(gj.fCone), fPtThreshold(gj.fPtThreshold),
196 fPtJetSelectionCut(gj.fPtJetSelectionCut),
197 fListHistos(0),//?????
198 fFastRec(gj.fFastRec), fOptFast(gj.fOptFast),
200 fResPara1(gj.fResPara1), fResPara2(gj.fResPara2), fResPara3(gj.fResPara3),
201 fPosParaA(gj.fPosParaA), fPosParaB(gj.fPosParaB),
202 fAngleMaxParam(gj.fAngleMaxParam), fSelect(gj.fSelect)
205 SetName (gj.GetName()) ;
206 SetTitle(gj.GetTitle()) ;
208 for(Int_t i = 0; i<10; i++){
209 fCones[i] = gj.fCones[i] ;
210 fNameCones[i] = gj.fNameCones[i] ;
211 fPtThres[i] = gj.fPtThres[i] ;
212 fNamePtThres[i] = gj.fNamePtThres[i] ;
214 fJetXMin1[i] = gj.fJetXMin1[i] ;
215 fJetXMin2[i] = gj.fJetXMin2[i] ;
216 fJetXMax1[i] = gj.fJetXMax1[i] ;
217 fJetXMax2[i] = gj.fJetXMax2[i] ;
218 fBkgMean[i] = gj.fBkgMean[i] ;
219 fBkgRMS[i] = gj.fBkgRMS[i] ;
221 fJetE1[i] = gj.fJetE1[i] ;
222 fJetE2[i] = gj.fJetE2[i] ;
223 fJetSigma1[i] = gj.fJetSigma1[i] ;
224 fJetSigma2[i] = gj.fJetSigma2[i] ;
225 fPhiEMCALCut[i] = gj.fPhiEMCALCut[i] ;
231 //____________________________________________________________________________
232 AliPHOSGammaJet::~AliPHOSGammaJet()
234 fOutputFile->Close() ;
237 //____________________________________________________________________________
238 void AliPHOSGammaJet::AddHIJINGToList(Int_t iEvent, TClonesArray * particleList,
241 TClonesArray * plNePHOS)
244 // List of particles copied to a root file.
246 // Char_t sb[] = "bgrd/";
247 // // cout<<sb<<endl;
249 // Char_t sf[] = "/list.root";
250 // // cout<<sf<<endl;
251 // sprintf(si,"%d",iEvent);
253 // // cout<<si<<endl;
255 // // cout<<si<<endl;
256 // TFile * f = TFile::Open(sb) ;
257 // //cout<<f->GetName()<<endl;
260 sprintf(fi,"bgrd/%d/list.root",iEvent);
261 TFile * f = TFile::Open(fi) ;
262 //cout<<f->GetName()<<endl;
264 TParticle *particle = new TParticle();
265 TTree *t = (TTree*) f->Get("T");
266 TBranch *branch = t->GetBranch("primaries");
267 branch->SetAddress(&particle);
269 Int_t index = particleList->GetEntries() ;
270 Int_t indexNe = plNe->GetEntries() ;
271 Int_t indexCh = plCh->GetEntries() ;
272 Int_t indexNePHOS = plNePHOS->GetEntries() ;
273 Double_t charge = 0.;
274 Int_t iParticle = 0 ;
276 Double_t x = 0., z = 0.;
277 // cout<<"bkg entries "<<t->GetEntries()<<endl;
279 AliPHOSGetter * gime = AliPHOSGetter::Instance() ;
280 const AliPHOSGeometry * geom = gime->PHOSGeometry() ;
282 //DP. To be fixed: extract vertex position
283 Double_t vtx[3]={0.,0.,0.} ;
286 for (iParticle=0 ; iParticle < t->GetEntries() ; iParticle++) {
287 t->GetEvent(iParticle) ;
293 charge = TDatabasePDG::Instance()
294 ->GetParticle(particle->GetPdgCode())->Charge();
296 if((charge != 0) && (particle->Pt() > fChargedPtCut)){
297 if(TMath::Abs(particle->Eta())<fEtaCut){
298 new((*particleList)[index]) TParticle(*particle) ;
299 (dynamic_cast<TParticle*>(particleList->At(index)))
303 new((*plCh)[indexCh]) TParticle(*particle) ;
304 (dynamic_cast<TParticle*>(plCh->At(indexCh)))->SetStatusCode(0) ;
307 if(strstr(fOptionGJ,"deb all")||strstr(fOptionGJ,"deb")){
308 dynamic_cast<TH1F*>(fListHistos->FindObject("PtSpectra"))->
309 Fill(particle->Pt());
312 else if((charge == 0) && (particle->Pt() > fNeutralPtCut) ){
313 geom->ImpactOnEmc(vtx,particle->Theta(),particle->Phi(), m,z,x);
317 if(strstr(fOptionGJ,"deb all")|| strstr(fOptionGJ,"deb"))
318 dynamic_cast<TH1F*>(fListHistos->FindObject("PtSpectra"))->
319 Fill(particle->Pt());
321 new((*plNePHOS)[indexNePHOS]) TParticle(*particle) ;
322 (dynamic_cast<TParticle*>(plNePHOS->At(indexNePHOS)))->SetStatusCode(0) ;
326 if((particle->Phi()>fPhiEMCALCut[0]) &&
327 (particle->Phi()<fPhiEMCALCut[1]) && m == 0)
329 if(strstr(fOptionGJ,"deb all")|| strstr(fOptionGJ,"deb"))
330 dynamic_cast<TH1F*>(fListHistos->FindObject("PtSpectra"))->
331 Fill(particle->Pt());
333 new((*particleList)[index]) TParticle(*particle) ;
334 (dynamic_cast<TParticle*>(particleList->At(index)))
337 new((*plNe)[indexNe]) TParticle(*particle) ;
338 (dynamic_cast<TParticle*>(plNe->At(indexNe)))->SetStatusCode(0) ;
347 Double_t mass = TDatabasePDG::Instance()->GetParticle(111)->Mass();
348 TLorentzVector pPi0, pGamma1, pGamma2 ;
349 Double_t angle = 0, cellDistance = 0.;
352 // fFastRec = new AliPHOSFastGlobalReconstruction(fHIJINGFileName);
353 // fFastRec->FastReconstruction(iiEvent);
355 for (iParticle=0 ; iParticle < t->GetEntries() ; iParticle++) {
356 t->GetEvent(iParticle) ;
360 charge = TDatabasePDG::Instance()
361 ->GetParticle(particle->GetPdgCode())->Charge();
363 if((charge != 0) && (particle->Pt() > fChargedPtCut)
364 && (TMath::Abs(particle->Eta())<fEtaCut)){
366 new((*particleList)[index]) TParticle(*particle) ;
367 (dynamic_cast<TParticle*>(particleList->At(index)))
371 new((*plCh)[indexCh]) TParticle(*particle) ;
372 (dynamic_cast<TParticle*>(plCh->At(indexCh)))->SetStatusCode(0) ;
375 else if(charge == 0){
376 geom->ImpactOnEmc(vtx,particle->Theta(),particle->Phi(), m,z,x);
377 if((particle->GetPdgCode() != 111) && (particle->Pt() > fNeutralPtCut) &&
378 (TMath::Abs(particle->Eta())<fEtaCut) ){
379 TLorentzVector part(particle->Px(),particle->Py(),
380 particle->Pz(),particle->Energy());
382 if(part.Pt() > fNeutralPtCut){
384 if(particle->Phi()>fPhiEMCALCut[0] &&
385 particle->Phi()<fPhiEMCALCut[1] && m == 0)
387 new((*particleList)[index]) TParticle(*particle) ;
388 (dynamic_cast<TParticle*>(particleList->At(index)))
390 (dynamic_cast<TParticle*>(particleList->At(index)))
391 ->SetMomentum(part.Px(),part.Py(),part.Pz(),part.E());
395 new((*plNe)[indexNe]) TParticle(*particle) ;
396 (dynamic_cast<TParticle*>(plNe->At(indexNe)))
397 ->SetMomentum(part.Px(),part.Py(),part.Pz(),part.E());
398 (dynamic_cast<TParticle*>(plNe->At(indexNe)))->SetStatusCode(0) ;
403 new((*plNePHOS)[indexNePHOS]) TParticle(*particle) ;
404 (dynamic_cast<TParticle*>(plNePHOS->At(indexNePHOS)))
405 ->SetMomentum(part.Px(),part.Py(),part.Pz(),part.E());
406 (dynamic_cast<TParticle*>(plNePHOS->At(indexNePHOS)))->SetStatusCode(0) ;
411 if((particle->GetPdgCode() == 111) && (particle->Pt() > fNeutralPtCut) &&
412 (TMath::Abs(particle->Eta())<fEtaCut+1))
415 pPi0.SetPxPyPzE(particle->Px(),particle->Py(),particle->Pz(),
420 Pi0Decay(mass,pPi0,pGamma1,pGamma2,angle);
421 //Check if decay photons are too close for PHOS
422 cellDistance = angle*460; //cm
423 if (cellDistance < fMinDistance) {
424 if(strstr(fOptionGJ,"deb all")|| strstr(fOptionGJ,"deb"))
425 dynamic_cast<TH1F*>(fListHistos->FindObject("PtSpectra"))->
426 Fill(particle->Pt());
428 //Pi0 inside phi EMCAL acceptance
430 TLorentzVector part(particle->Px(),particle->Py(),
431 particle->Pz(),particle->Energy());
433 if(part.Pt() > fNeutralPtCut){
434 if(particle->Phi()>fPhiEMCALCut[0] &&
435 particle->Phi()<fPhiEMCALCut[1] && m == 0){
437 new((*particleList)[index]) TParticle(*particle) ;
438 (dynamic_cast<TParticle*>(particleList->At(index)))->SetStatusCode(0) ;
439 (dynamic_cast<TParticle*>(particleList->At(index)))
440 ->SetMomentum(part.Px(),part.Py(),part.Pz(),part.E());
443 new((*plNe)[indexNe]) TParticle(*particle) ;
444 (dynamic_cast<TParticle*>(plNe->At(indexNe))) ->SetStatusCode(0) ;
445 (dynamic_cast<TParticle*>(plNe->At(indexNe)))
446 ->SetMomentum(part.Px(),part.Py(),part.Pz(),part.E());
450 if(strstr(fOptionGJ,"deb all")|| strstr(fOptionGJ,"deb"))
451 dynamic_cast<TH1F*>(fListHistos->FindObject("PtSpectra"))->
452 Fill(particle->Pt());
453 new((*plNePHOS)[indexNePHOS]) TParticle(*particle) ;
454 (dynamic_cast<TParticle*>(plNePHOS->At(indexNePHOS))) ->SetStatusCode(0) ;
455 (dynamic_cast<TParticle*>(plNePHOS->At(indexNePHOS)))
456 ->SetMomentum(part.Px(),part.Py(),part.Pz(),part.E());
463 dynamic_cast<TH2F*>(fListHistos->FindObject("AnglePair"))
464 ->Fill(pPi0.E(),angle);
467 if(pGamma1.Pt() > 0. && TMath::Abs(pGamma1.Eta())<fEtaCut){
471 if(pGamma1.Pt() > fNeutralPtCut ){
473 TParticle * photon1 =
474 new TParticle(22,1,0,0,0,0,pGamma1.Px(),pGamma1.Py(),
475 pGamma1.Pz(),pGamma1.E(),0,0,0,0);
476 geom->ImpactOnEmc(vtx,photon1->Theta(),photon1->Phi(), m,z,x);
477 if( photon1->Phi()>fPhiEMCALCut[0] && photon1->Phi()<fPhiEMCALCut[1]
479 if(strstr(fOptionGJ,"deb all") || strstr(fOptionGJ,"deb"))
480 dynamic_cast<TH1F*>(fListHistos->FindObject("PtSpectra"))->
482 new((*particleList)[index]) TParticle(*photon1) ;
483 (dynamic_cast<TParticle*>(particleList->At(index)))->SetStatusCode(0) ;
486 new((*plNe)[indexNe]) TParticle(*photon1) ;
487 (dynamic_cast<TParticle*>(plNe->At(indexNe))) ->SetStatusCode(0) ;
492 if(strstr(fOptionGJ,"deb all") || strstr(fOptionGJ,"deb"))
493 dynamic_cast<TH1F*>(fListHistos->FindObject("PtSpectra"))->
495 new((*plNePHOS)[indexNePHOS]) TParticle(*photon1) ;
496 (dynamic_cast<TParticle*>(plNePHOS->At(indexNePHOS)))->SetStatusCode(0) ;
502 if(pGamma2.Pt() > 0. && TMath::Abs(pGamma2.Eta())<fEtaCut){
505 if(pGamma2.Pt() > fNeutralPtCut){
507 TParticle * photon2 =
508 new TParticle(22,1,0,0,0,0,pGamma2.Px(), pGamma2.Py(),
509 pGamma2.Pz(),pGamma2.E(),0,0,0,0);
510 geom->ImpactOnEmc(vtx,photon2->Theta(),photon2->Phi(), m,z,x);
511 if(photon2->Phi()>fPhiEMCALCut[0] &&
512 photon2->Phi()<fPhiEMCALCut[1] && m == 0){
513 if(strstr(fOptionGJ,"deb all") || strstr(fOptionGJ,"deb"))
514 dynamic_cast<TH1F*>(fListHistos->FindObject("PtSpectra"))->
516 new((*particleList)[index]) TParticle(*photon2) ;
517 (dynamic_cast<TParticle*>(particleList->At(index)))->SetStatusCode(0) ;
520 new((*plNe)[indexNe]) TParticle(*photon2) ;
521 (dynamic_cast<TParticle*>(plNe->At(indexNe))) ->SetStatusCode(0) ;
525 if(strstr(fOptionGJ,"deb all") || strstr(fOptionGJ,"deb"))
526 dynamic_cast<TH1F*>(fListHistos->FindObject("PtSpectra"))->
528 new((*plNePHOS)[indexNePHOS]) TParticle(*photon2) ;
529 (dynamic_cast<TParticle*>(plNePHOS->At(indexNePHOS))) ->SetStatusCode(0) ;
533 // e = (pGamma1+pGamma2).E();
534 // if(IsAngleInWindow(angle,e))
536 (fListHistos->FindObject("AnglePairAccepted"))->
537 Fill(pPi0.E(),angle);
540 }//photon2 in acceptance
541 }//if angle > mindist
544 }//for (iParticle<nParticle)
547 //Info("AddHIJINGToList","End HIJING");
550 //____________________________________________________________________________
551 Double_t AliPHOSGammaJet::CalculateJetRatioLimit(const Double_t ptg,
555 //Info("CalculateLimit","x1 %f, x2%f",x[0],x[1]);
556 Double_t epp = par[0] + par[1] * ptg ;
557 Double_t spp = par[2] + par[3] * ptg ;
558 Double_t f = x[0] + x[1] * ptg ;
559 Double_t epb = epp + par[4] ;
560 Double_t spb = TMath::Sqrt(spp*spp+ par[5]*par[5]) ;
561 Double_t rat = (epb - spb * f) / ptg ;
562 //Info("CalculateLimit","epp %f, spp %f, f %f", epp, spp, f);
563 //Info("CalculateLimit","epb %f, spb %f, rat %f", epb, spb, rat);
567 //____________________________________________________________________________
568 void AliPHOSGammaJet::CreateParticleList(Int_t iEvent,
569 TClonesArray * particleList,
572 TClonesArray * plNePHOS)
574 //Info("CreateParticleList","Inside");
575 AliPHOSGetter * gime = AliPHOSGetter::Instance(fInputFileName) ;
576 const AliPHOSGeometry * geom = gime->PHOSGeometry() ;
577 gime->Event(iEvent, "X") ;
579 //DP to be fixed: extract true vertex here
580 Double_t vtx[3]={0.,0.,0.} ;
582 Int_t index = particleList->GetEntries() ;
583 Int_t indexCh = plCh->GetEntries() ;
584 Int_t indexNe = plNe->GetEntries() ;
585 Int_t indexNePHOS = plNePHOS->GetEntries() ;
586 Int_t iParticle = 0 ;
587 Double_t charge = 0.;
589 Double_t x = 0., z = 0.;
592 for (iParticle=0 ; iParticle < gime->NPrimaries() ; iParticle++) {
593 const TParticle * particle = gime->Primary(iParticle) ;
599 //Keep Stable particles within eta range
600 if((particle->GetStatusCode() == 1) &&
601 (particle->Pt() > 0)){
602 if(TMath::Abs(particle->Eta())<fEtaCut){
605 charge = TDatabasePDG::Instance()
606 ->GetParticle(particle->GetPdgCode())->Charge();
607 if((charge != 0) && (particle->Pt() > fChargedPtCut)){
609 if(strstr(fOptionGJ,"deb all")|| strstr(fOptionGJ,"deb"))
610 dynamic_cast<TH1F*>(fListHistos->FindObject("PtSpectra"))->
611 Fill(particle->Pt());
612 new((*plCh)[indexCh++]) TParticle(*particle) ;
613 new((*particleList)[index++]) TParticle(*particle) ;
615 else if((charge == 0) && (particle->Pt() > fNeutralPtCut)){
616 geom->ImpactOnEmc(vtx,particle->Theta(),particle->Phi(), m,z,x);
619 if(strstr(fOptionGJ,"deb all")|| strstr(fOptionGJ,"deb"))
620 dynamic_cast<TH1F*>(fListHistos->FindObject("PtSpectra"))->
621 Fill(particle->Pt());
623 new((*plNePHOS)[indexNePHOS++]) TParticle(*particle) ;
625 if((particle->Phi()>fPhiEMCALCut[0]) &&
626 (particle->Phi()<fPhiEMCALCut[1]) && m == 0)
628 if(strstr(fOptionGJ,"deb all")|| strstr(fOptionGJ,"deb"))
629 dynamic_cast<TH1F*>(fListHistos->FindObject("PtSpectra"))->
630 Fill(particle->Pt());
631 new((*plNe)[indexNe++]) TParticle(*particle) ;
632 new((*particleList)[index++]) TParticle(*particle) ;
636 }//final particle, etacut
637 }//for (iParticle<nParticle)
641 Double_t mass = TDatabasePDG::Instance()->GetParticle(111)->Mass();
642 TLorentzVector pPi0, pGamma1, pGamma2 ;
643 Double_t angle = 0, cellDistance = 0.;
645 fFastRec = new AliPHOSFastGlobalReconstruction(fInputFileName);
646 fFastRec->FastReconstruction(iEvent);
648 for (iParticle=0 ; iParticle < gime->NPrimaries() ; iParticle++) {
649 const TParticle * particle = gime->Primary(iParticle) ;
653 //Keep Stable particles within eta range
654 if((particle->GetStatusCode() == 1) && (particle->Pt() > 0)){
658 charge = TDatabasePDG::Instance()
659 ->GetParticle(particle->GetPdgCode())->Charge();
660 if((charge != 0) && (particle->Pt() > fChargedPtCut) && (TMath::Abs(particle->Eta())<fEtaCut)){
661 new((*plCh)[indexCh++]) TParticle(*particle) ;
662 new((*particleList)[index++]) TParticle(*particle) ;
664 else if(charge == 0) {
665 geom->ImpactOnEmc(vtx,particle->Theta(),particle->Phi(), m,z,x);
666 if((particle->GetPdgCode() != 111) && particle->Pt() > 0 &&
667 (TMath::Abs(particle->Eta())<fEtaCut))
670 TLorentzVector part(particle->Px(),particle->Py(),
671 particle->Pz(),particle->Energy());
675 if(part.Pt() > fNeutralPtCut){
676 if(particle->Phi()>fPhiEMCALCut[0] &&
677 particle->Phi()<fPhiEMCALCut[1] && m == 0)
679 new((*particleList)[index]) TParticle(*particle) ;
680 (dynamic_cast<TParticle*>(particleList->At(index)))
681 ->SetMomentum(part.Px(),part.Py(),part.Pz(),part.E());
684 new((*plNe)[indexNe]) TParticle(*particle) ;
685 (dynamic_cast<TParticle*>(plNe->At(indexNe)))
686 ->SetMomentum(part.Px(),part.Py(),part.Pz(),part.E());
691 new((*plNePHOS)[indexNePHOS]) TParticle(*particle) ;
692 (dynamic_cast<TParticle*>(plNePHOS->At(indexNePHOS)))
693 ->SetMomentum(part.Px(),part.Py(),part.Pz(),part.E());
698 if((particle->GetPdgCode() == 111) && (particle->Pt() > fNeutralPtCut) &&
699 (TMath::Abs(particle->Eta())<fEtaCut+1))
702 pPi0.SetPxPyPzE(particle->Px(),particle->Py(),particle->Pz(),
707 Pi0Decay(mass,pPi0,pGamma1,pGamma2,angle);
708 //Check if decay photons are too close for PHOS
709 cellDistance = angle*460; //cm
711 if (cellDistance < fMinDistance) {
713 //Pi0 inside phi EMCAL acceptance
716 TLorentzVector part(particle->Px(),particle->Py(),
717 particle->Pz(),particle->Energy());
720 if(part.Pt() > fNeutralPtCut){
721 if(particle->Phi()>fPhiEMCALCut[0] &&
722 particle->Phi()<fPhiEMCALCut[1] && m == 0){
723 if(strstr(fOptionGJ,"deb all")|| strstr(fOptionGJ,"deb"))
724 dynamic_cast<TH1F*>(fListHistos->FindObject("PtSpectra"))->
725 Fill(particle->Pt());
727 new((*plNe)[indexNe]) TParticle(*particle) ;
728 (dynamic_cast<TParticle*>(plNe->At(indexNe)))
729 ->SetMomentum(part.Px(),part.Py(),part.Pz(),part.E());
730 new((*particleList)[index]) TParticle(*particle) ;
731 (dynamic_cast<TParticle*>(particleList->At(index)))
732 ->SetMomentum(part.Px(),part.Py(),part.Pz(),part.E());
737 if(strstr(fOptionGJ,"deb all")|| strstr(fOptionGJ,"deb"))
738 dynamic_cast<TH1F*>(fListHistos->FindObject("PtSpectra"))->
739 Fill(particle->Pt());
740 new((*plNePHOS)[indexNePHOS]) TParticle(*particle) ;
741 (dynamic_cast<TParticle*>(plNePHOS->At(indexNePHOS)))
742 ->SetMomentum(part.Px(),part.Py(),part.Pz(),part.E());
749 dynamic_cast<TH2F*>(fListHistos->FindObject("AnglePair"))
750 ->Fill(pPi0.E(),angle);
754 if(pGamma1.Pt() > 0 && TMath::Abs(pGamma1.Eta())<fEtaCut){
757 if(pGamma1.Pt() > fNeutralPtCut){
758 TParticle * photon1 =
759 new TParticle(22,1,0,0,0,0,pGamma1.Px(),pGamma1.Py(),
760 pGamma1.Pz(),pGamma1.E(),0,0,0,0);
761 geom->ImpactOnEmc(vtx,photon1->Theta(),photon1->Phi(), m,z,x);
762 if( photon1->Phi()>fPhiEMCALCut[0] && photon1->Phi()<fPhiEMCALCut[1]
764 if(strstr(fOptionGJ,"deb all") || strstr(fOptionGJ,"deb"))
765 dynamic_cast<TH1F*>(fListHistos->FindObject("PtSpectra"))->
767 new((*plNe)[indexNe++]) TParticle(*photon1) ;
768 new((*particleList)[index++]) TParticle(*photon1) ;
773 if(strstr(fOptionGJ,"deb all") || strstr(fOptionGJ,"deb"))
774 dynamic_cast<TH1F*>(fListHistos->FindObject("PtSpectra"))->
776 new((*plNePHOS)[indexNePHOS++]) TParticle(*photon1) ;
781 if(pGamma2.Pt() > 0 && TMath::Abs(pGamma2.Eta())<fEtaCut){
784 if(pGamma2.Pt() > fNeutralPtCut ){
785 TParticle * photon2 =
786 new TParticle(22,1,0,0,0,0,pGamma2.Px(), pGamma2.Py(),
787 pGamma2.Pz(),pGamma2.E(),0,0,0,0);
788 geom->ImpactOnEmc(vtx,photon2->Theta(),photon2->Phi(), m,z,x);
789 if(photon2->Phi()>fPhiEMCALCut[0] &&
790 photon2->Phi()<fPhiEMCALCut[1] && m == 0){
791 if(strstr(fOptionGJ,"deb all") || strstr(fOptionGJ,"deb"))
792 dynamic_cast<TH1F*>(fListHistos->FindObject("PtSpectra"))->
794 new((*plNe)[indexNe++]) TParticle(*photon2) ;
795 new((*particleList)[index++]) TParticle(*photon2) ;
798 if(strstr(fOptionGJ,"deb all") || strstr(fOptionGJ,"deb"))
799 dynamic_cast<TH1F*>(fListHistos->FindObject("PtSpectra"))->
801 new((*plNePHOS)[indexNePHOS++]) TParticle(*photon2) ;
805 // Float_t e = (pGamma1+pGamma2).E();
806 // if(IsAngleInWindow(angle,e))
808 (fListHistos->FindObject("AnglePairAccepted"))->
809 Fill(pPi0.E(),angle);
812 }//photon2 in acceptance
813 }//if angle > mindist
816 }//final particle, etacut
817 }//for (iParticle<nParticle)
821 //____________________________________________________________________________
822 void AliPHOSGammaJet::CreateParticleListFromESD(TClonesArray * pl,
825 TClonesArray * plNePHOS,
828 //Create a list of particles from the ESD. These particles have been measured
829 //by the Central Tracking system (TPC), PHOS and EMCAL
830 //(EMCAL not available for the moment).
831 //Info("CreateParticleListFromESD","Inside");
833 Int_t index = pl->GetEntries() ;
835 Float_t *pid = new Float_t[AliPID::kSPECIESN];
837 //########### PHOS ##############
838 //Info("CreateParticleListFromESD","Fill ESD PHOS list");
839 Int_t begphos = esd->GetFirstPHOSCluster();
840 Int_t endphos = esd->GetFirstPHOSCluster() +
841 esd->GetNumberOfPHOSClusters() ;
842 Int_t indexNePHOS = plNePHOS->GetEntries() ;
843 if(strstr(fOptionGJ,"deb all"))
844 Info("CreateParticleListFromESD","PHOS: first particle %d, last particle %d",
847 for (npar = begphos; npar < endphos; npar++) {//////////////PHOS track loop
848 AliESDCaloCluster * clus = esd->GetCaloCluster(npar) ; // retrieve track from esd
850 //Create a TParticle to fill the particle list
852 Float_t en = clus->GetClusterEnergy() ;
853 Float_t *p = new Float_t();
854 clus->GetGlobalPosition(p) ;
855 TVector3 pos(p[0],p[1],p[2]) ;
856 Double_t phi = pos.Phi();
857 Double_t theta= pos.Theta();
858 Double_t px = en*TMath::Cos(phi)*TMath::Sin(theta);;
859 Double_t py = en*TMath::Sin(phi)*TMath::Sin(theta);
860 Double_t pz = en*TMath::Cos(theta);
862 TParticle * particle = new TParticle() ;
863 particle->SetMomentum(px,py,pz,en) ;
865 //Select only photons
868 //cout<<"pid "<<pid[AliPID::kPhoton]<<endl ;
869 if( pid[AliPID::kPhoton] > 0.75)
870 new((*plNePHOS)[indexNePHOS++]) TParticle(*particle) ;
873 //########### TPC #####################
874 //Info("CreateParticleListFromESD","Fill ESD TPC list");
876 Int_t endtpc = esd->GetNumberOfTracks() ;
877 Int_t indexCh = plCh->GetEntries() ;
878 if(strstr(fOptionGJ,"deb all"))
879 Info("CreateParticleListFromESD","TPC: first particle %d, last particle %d",
882 for (npar = begtpc; npar < endtpc; npar++) {////////////// track loop
883 AliESDtrack * track = esd->GetTrack(npar) ; // retrieve track from esd
885 Double_t en = track ->GetTPCsignal() ;
887 track->GetPxPyPz(mom) ;
888 Double_t px = mom[0];
889 Double_t py = mom[1];
890 Double_t pz = mom[2]; //Check with TPC people if this is correct.
892 //cout<<"TPC signal "<<en<<endl;
893 //cout<<"px "<<px<<"; py "<<py<<"; pz "<<pz<<endl;
894 TParticle * particle = new TParticle() ;
895 particle->SetMomentum(px,py,pz,en) ;
897 new((*plCh)[indexCh++]) TParticle(*particle) ;
898 new((*pl)[index++]) TParticle(*particle) ;
902 //################ EMCAL ##############
903 Double_t v[3] ; //vertex ;
904 esd->GetVertex()->GetXYZ(v) ;
905 //##########Uncomment when ESD for EMCAL works ##########
906 //Info("CreateParticleListFromESD","Fill ESD EMCAL list");
908 Int_t begem = esd->GetFirstEMCALCluster();
909 Int_t endem = esd->GetFirstEMCALCluster() +
910 esd->GetNumberOfEMCALClusters() ;
911 Int_t indexNe = plNe->GetEntries() ;
912 if(strstr(fOptionGJ,"deb all"))
913 Info("CreateParticleListFromESD","EMCAL: first particle %d, last particle %d",
916 for (npar = begem; npar < endem; npar++) {//////////////EMCAL track loop
917 AliESDCaloCluster * clus = esd->GetCaloCluster(npar) ; // retrieve track from esd
919 Float_t en = clus->GetClusterEnergy() ;
920 Float_t *p = new Float_t();
921 clus->GetGlobalPosition(p) ;
922 TVector3 pos(p[0],p[1],p[2]) ;
923 Double_t phi = pos.Phi();
924 Double_t theta= pos.Theta();
925 Double_t px = en*TMath::Cos(phi)*TMath::Sin(theta);;
926 Double_t py = en*TMath::Sin(phi)*TMath::Sin(theta);
927 Double_t pz = en*TMath::Cos(theta);
928 //cout<<"EMCAL signal "<<en<<endl;
929 //cout<<"px "<<px<<"; py "<<py<<"; pz "<<pz<<endl;
930 //TParticle * particle = new TParticle() ;
931 //particle->SetMomentum(px,py,pz,en) ;
934 // //Uncomment if PID IS WORKING, photon and pi0 idenitification.
935 // //if( pid[AliPID::kPhoton] > 0.75) //This has to be fixen.
937 // //else if( pid[AliPID::kPi0] > 0.75)
939 pdg = 22; //No PID, assume all photons
940 TParticle * particle = new TParticle(pdg, 1, -1, -1, -1, -1,
941 px, py, pz, en, v[0], v[1], v[2], 0);
943 new((*plNe)[indexNe++]) TParticle(*particle) ;
944 new((*pl)[index++]) TParticle(*particle) ;
947 // Info("CreateParticleListFromESD","End Inside");
953 //____________________________________________________________________________
954 void AliPHOSGammaJet::Exec(Option_t *option)
964 // Create chain of esd trees
965 const UInt_t kNevent = static_cast<UInt_t>(GetNEvent()) ;
966 t = ReadESD(kNevent, fDirName, fESDTree, fPattern) ;
968 AliError("Could not create the TChain") ;
973 t->SetBranchAddress("ESD",&esd); // point to the container esd where to put the event from the esdTree
978 // AliGenPythia* pyth = (AliGenPythia*) gAlice->Generator();
981 TClonesArray * particleList = new TClonesArray("TParticle",1000);
982 TClonesArray * plCh = new TClonesArray("TParticle",1000);
983 TClonesArray * plNe = new TClonesArray("TParticle",1000);
984 TClonesArray * plNePHOS = new TClonesArray("TParticle",1000);
986 for (Int_t iEvent = 0 ; iEvent < fNEvent ; iEvent++) {
987 if(strstr(fOptionGJ,"deb")||strstr(fOptionGJ,"deb all"))
988 Info("Exec", "Event %d", iEvent) ;
992 Double_t phig = 0., phil = 0., phich = 0 , phipi = 0;
993 Double_t etag = 0., etal = 0., etach = 0., etapi = 0. ;
994 Double_t ptg = 0., ptl = 0., ptch = 0., ptpi = 0.;
996 TLorentzVector jet (0,0,0,0);
997 TLorentzVector jettpc(0,0,0,0);
1001 Int_t iNbytes = t->GetEntry(iEvent); // store event in esd
1002 //cout<<"nbytes "<<iNbytes<<endl;
1003 if ( iNbytes == 0 ) {
1004 AliError("Empty TChain") ;
1007 CreateParticleListFromESD(particleList, plCh,plNe,plNePHOS, esd); //,iEvent);
1010 CreateParticleList(iEvent, particleList, plCh,plNe,plNePHOS);
1012 // TLorentzVector pyjet(0,0,0,0);
1015 // Float_t jets[4][10];
1016 // pyth->SetJetReconstructionMode(1);
1017 // pyth->LoadEvent();
1018 // pyth->GetJets(nJ, nJT, jets);
1020 // Float_t pxJ = jets[0][0];
1021 // Float_t pyJ = jets[1][0];
1022 // Float_t pzJ = jets[2][0];
1023 // Float_t eJ = jets[3][0];
1024 // pyjet.SetPxPyPzE(pxJ,pyJ,pzJ,eJ ) ;
1027 // //Info("Exec",">>>>>>>>>>Number of jets !!!! %d",nJT);
1028 // for (Int_t iJ = 1; iJ < nJT; iJ++) {
1029 // Float_t pxJ = jets[0][iJ];
1030 // Float_t pyJ = jets[1][iJ];
1031 // Float_t pzJ = jets[2][iJ];
1032 // Float_t eJ = jets[3][iJ];
1033 // pyjet.SetPxPyPzE(pxJ,pyJ,pzJ,eJ ) ;
1034 // //Info("Exec",">>>>>Pythia Jet: %d, Phi %f, Eta %f, Pt %f",
1035 // // iJ,pyjet.Phi(),pyjet.Eta(),pyjet.Pt());
1041 AddHIJINGToList(iEvent, particleList, plCh,plNe, plNePHOS);
1044 Bool_t iIsInPHOS = kFALSE ;
1045 GetGammaJet(plNePHOS, ptg, phig, etag, iIsInPHOS) ;
1049 //Info("Exec"," In PHOS") ;
1050 dynamic_cast<TH1F*>(fListHistos->FindObject("NGamma"))->Fill(ptg);
1051 dynamic_cast<TH2F*>(fListHistos->FindObject("PhiGamma"))
1053 dynamic_cast<TH2F*>(fListHistos->FindObject("EtaGamma"))
1055 if(strstr(fOptionGJ,"deb")||strstr(fOptionGJ,"deb all"))
1056 Info("Exec", "Gamma: pt %f, phi %f, eta %f", ptg,
1059 // cout<<"n charged "<<plCh->GetEntries()<<endl;
1060 // cout<<"n neutral "<<plNe->GetEntries()<<endl;
1061 // cout<<"n All "<<particleList->GetEntries()<<endl;
1063 GetLeadingCharge(plCh, ptg, phig, ptch, etach, phich) ;
1064 GetLeadingPi0 (plNe, ptg, phig, ptpi, etapi, phipi) ;
1066 // cout<<"n2 charged "<<plCh->GetEntries()<<endl;
1067 // cout<<"n2 neutral "<<plNe->GetEntries()<<endl;
1068 // cout<<"n2 All "<<particleList->GetEntries()<<endl;
1073 //Is the leading cone inside EMCAL?
1074 Bool_t insidech = kFALSE ;
1075 if((phich - fCone) > fPhiEMCALCut[0] &&
1076 (phich + fCone) < fPhiEMCALCut[1]){
1079 Bool_t insidepi = kFALSE ;
1080 if((phipi - fCone) > fPhiEMCALCut[0] &&
1081 (phipi + fCone) < fPhiEMCALCut[1]){
1085 if ((ptch > 0 || ptpi > 0)){
1086 if((ptch > ptpi) && insidech){
1090 dynamic_cast<TH2F*>(fListHistos->FindObject("ChargeRatio"))
1091 ->Fill(ptg,ptch/ptg);
1092 dynamic_cast<TH2F*>(fListHistos->FindObject("DeltaPhiCharge"))
1093 ->Fill(ptg,phig-phich);
1094 dynamic_cast<TH2F*>(fListHistos->FindObject("DeltaEtaCharge"))
1095 ->Fill(ptg,etag-etach);
1096 if(strstr(fOptionGJ,"deb"))
1097 Info("Exec"," Charged Leading") ;
1099 if((ptpi > ptch) && insidepi){
1104 dynamic_cast<TH2F*>(fListHistos->FindObject("Pi0Ratio"))
1105 ->Fill(ptg,ptpi/ptg);
1106 dynamic_cast<TH2F*>(fListHistos->FindObject("DeltaPhiPi0"))
1107 ->Fill(ptg,phig-phipi);
1108 dynamic_cast<TH2F*>(fListHistos->FindObject("DeltaEtaPi0"))
1109 ->Fill(ptg,etag-etapi);
1111 if(ptpi > 0. && strstr(fOptionGJ,"deb"))
1112 Info("Exec"," Pi0 Leading") ;
1115 if(strstr(fOptionGJ,"deb"))
1116 Info("Exec","Leading pt %f, phi %f",ptl,phil);
1117 if(insidech || insidepi){
1120 MakeJet(particleList, ptg, phig, ptl, phil, etal, "", jet);
1122 if(strstr(fOptionGJ,"deb")){
1123 // Info("Exec","Pythia Jet: Phi %f, Eta %f, Pt %f",
1124 // pyjet.Phi(),pyjet.Eta(),pyjet.Pt());
1125 Info("Exec","TPC+EMCAL Jet: Phi %f, Eta %f, Pt %f",
1126 jet.Phi(),jet.Eta(),jet.Pt());
1128 // dynamic_cast<TH2F*>(fListHistos->FindObject("DeltaPhiJet"))
1129 // ->Fill(ptg,pyjet.Phi()-jet.Phi());
1130 // dynamic_cast<TH2F*>(fListHistos->FindObject("DeltaEtaJet"))
1131 // ->Fill(ptg,pyjet.Eta()-jet.Eta());
1132 // dynamic_cast<TH2F*>(fListHistos->FindObject("DeltaPtJet"))
1133 // ->Fill(ptg,pyjet.Pt()-jet.Pt());
1136 MakeJetAnyConeOrPt(particleList, ptg, phig, ptl, phil, etal, "");
1140 if(fOnlyCharged && ptch > 0.)
1142 if(strstr(fOptionGJ,"deb"))
1143 Info("Exec","Leading TPC pt %f, phi %f",ptch,phich);
1145 dynamic_cast<TH2F*>(fListHistos->FindObject("TPCRatio"))
1146 ->Fill(ptg,ptch/ptg);
1147 dynamic_cast<TH2F*>(fListHistos->FindObject("DeltaPhiTPC"))
1148 ->Fill(ptg,phig-phich);
1149 dynamic_cast<TH2F*>(fListHistos->FindObject("DeltaEtaTPC"))
1150 ->Fill(ptg,etag-etach);
1154 MakeJet(plCh, ptg, phig, ptch, phich, etach, "TPC",jettpc);
1156 if(strstr(fOptionGJ,"deb")){
1157 // Info("Exec","Pythia Jet: Phi %f, Eta %f, Pt %f",
1158 // pyjet.Phi(),pyjet.Eta(),pyjet.Pt());
1159 Info("Exec","TPC Jet: Phi %f, Eta %f, Pt %f",
1160 jettpc.Phi(),jettpc.Eta(),jettpc.Pt());
1162 // dynamic_cast<TH2F*>(fListHistos->FindObject("DeltaPhiTPCJet"))
1163 // ->Fill(ptg,pyjet.Phi()-jettpc.Phi());
1164 // dynamic_cast<TH2F*>(fListHistos->FindObject("DeltaEtaTPCJet"))
1165 // ->Fill(ptg,pyjet.Eta()-jettpc.Eta());
1166 // dynamic_cast<TH2F*>(fListHistos->FindObject("DeltaPtTPCJet"))
1167 // ->Fill(ptg,pyjet.Pt()-jettpc.Pt());
1170 MakeJetAnyConeOrPt(plCh, ptg, phig, ptch, phich, etach, "TPC");
1176 particleList->Delete() ;
1179 plNePHOS->Delete() ;
1184 delete particleList ;
1186 fOutputFile->Write() ;
1191 //____________________________________________________________________________
1192 void AliPHOSGammaJet::FillJetHistos(TClonesArray * pl, Double_t ptg,
1193 TString conf, TString type)
1195 //Fill jet fragmentation histograms if !fAnyCone,
1196 //only for fCone and fPtThres
1197 TParticle * particle = 0 ;
1202 while ( (particle = dynamic_cast<TParticle*>(next())) ) {
1204 Double_t pt = particle->Pt();
1206 charge = TDatabasePDG::Instance()
1207 ->GetParticle(particle->GetPdgCode())->Charge();
1208 if(charge != 0){//Only jet Charged particles
1210 (fListHistos->FindObject(type+conf+"Fragment"))
1213 (fListHistos->FindObject(type+conf+"PtDist"))
1219 (fListHistos->FindObject("NBkg"+conf))->Fill(ipr);
1222 //____________________________________________________________________________
1223 void AliPHOSGammaJet::FillJetHistosAnyConeOrPt(TClonesArray * pl, Double_t ptg,
1224 TString conf, TString type,
1225 TString cone, TString ptcut)
1227 //Fill jet fragmentation histograms if fAnyCone,
1228 //for several cones and pt thresholds
1229 TParticle *particle = 0;
1234 while ( (particle = dynamic_cast<TParticle*>(next())) ) {
1236 Double_t pt = particle->Pt();
1237 charge = TDatabasePDG::Instance()
1238 ->GetParticle(particle->GetPdgCode())->Charge();
1239 if(charge != 0){//Only jet Charged particles
1241 (fListHistos->FindObject(type+conf+"FragmentCone"+cone+"Pt"+ptcut))
1244 (fListHistos->FindObject(type+conf+"PtDistCone"+cone+"Pt"+ptcut))
1251 (fListHistos->FindObject("NBkg"+conf+"Cone"+cone+"Pt"+ptcut))
1256 //____________________________________________________________________________
1257 void AliPHOSGammaJet::GetGammaJet(TClonesArray * pl, Double_t &pt,
1258 Double_t &phi, Double_t &eta, Bool_t &Is) const
1260 //Search for the prompt photon in PHOS with pt > fPtCut
1265 for(Int_t ipr = 0;ipr < pl->GetEntries() ; ipr ++ ){
1266 TParticle * particle = dynamic_cast<TParticle *>(pl->At(ipr)) ;
1268 if((particle->Pt() > fPtCut) && (particle->Pt() > pt)){
1270 pt = particle->Pt();
1271 phi = particle->Phi() ;
1272 eta = particle->Eta() ;
1278 //____________________________________________________________________________
1279 void AliPHOSGammaJet::GetLeadingCharge(TClonesArray * pl,
1280 Double_t ptg, Double_t phig,
1281 Double_t &pt, Double_t &eta, Double_t &phi) const
1283 //Search for the charged particle with highest with
1284 //Phi=Phi_gamma-Pi and pT=0.1E_gamma
1289 for(Int_t ipr = 0;ipr < pl->GetEntries() ; ipr ++ ){
1291 TParticle * particle = dynamic_cast<TParticle *>(pl->At(ipr)) ;
1293 Double_t ptl = particle->Pt();
1294 Double_t rat = ptl/ptg ;
1295 Double_t phil = particle->Phi() ;
1297 if(((phig-phil)> fPhiMinCut) && ((phig-phil)<fPhiMaxCut) &&
1298 (rat > fRatioMinCut) && (rat < fRatioMaxCut) && (ptl > pt)) {
1299 eta = particle->Eta() ;
1302 //printf("GetLeadingCharge: %f %f %f %f \n", pt, eta, phi,rat) ;
1305 //printf("GetLeadingCharge: %f %f %f \n", pt, eta, phi) ;
1310 //____________________________________________________________________________
1311 void AliPHOSGammaJet::GetLeadingPi0(TClonesArray * pl,
1312 Double_t ptg, Double_t phig,
1313 Double_t &pt, Double_t &eta, Double_t &phi)
1316 //Search for the neutral pion with highest with
1317 //Phi=Phi_gamma-Pi and pT=0.1E_gamma
1321 Double_t ptl = -100.;
1322 Double_t rat = -100.;
1323 Double_t phil = -100. ;
1326 TParticle * particle = 0;
1330 while ( (particle = (TParticle*)next()) ) {
1331 if( particle->GetPdgCode() == 111){
1332 ptl = particle->Pt();
1334 phil = particle->Phi() ;
1335 e = particle->Energy();
1337 (fListHistos->FindObject("AnglePairNoCut"))->
1340 (fListHistos->FindObject("InvMassPairNoCut"))->
1343 if(((phig-phil)> fPhiMinCut) && ((phig-phil)<fPhiMaxCut) &&
1344 (rat > fRatioMinCut) && (rat < fRatioMaxCut)) {
1347 (fListHistos->FindObject("AnglePairLeadingCut"))->
1350 (fListHistos->FindObject("InvMassPairLeadingCut"))->
1354 (fListHistos->FindObject("AnglePairAngleCut"))->
1357 (fListHistos->FindObject("InvMassPairAngleCut"))->
1361 (fListHistos->FindObject("InvMassPairAllCut"))->
1364 (fListHistos->FindObject("AnglePairAllCut"))->
1369 eta = particle->Eta() ;
1373 //printf("GetLeadingPi0: %f %f %f %f %f \n", pt, eta, phi, rat, ptg) ;
1379 (fListHistos->FindObject("InvMassPairLeading"))->
1382 (fListHistos->FindObject("AnglePairLeading"))->
1387 Int_t iPrimary = -1;
1388 TLorentzVector gammai,gammaj;
1389 Double_t angle = 0., e = 0., invmass = 0.;
1390 Double_t anglef = 0., ef = 0., invmassf = 0.;
1394 while ( (particle = (TParticle*)next()) ) {
1397 ksPdg = particle->GetPdgCode();
1398 ptl = particle->Pt();
1399 if(ksPdg == 111){ //2 gamma
1401 phil = particle->Phi() ;
1402 if((ptl> pt)&& (rat > fRatioMinCut) && (rat < fRatioMaxCut) &&
1403 ((phig-phil)>fPhiMinCut)&&((phig-phil)<fPhiMaxCut)){
1404 eta = particle->Eta() ;
1409 if(ksPdg == 22){//1 gamma
1410 particle->Momentum(gammai);
1413 while ( (particle = (TParticle*)next2()) ) {
1415 if(jPrimary>iPrimary){
1416 ksPdg = particle->GetPdgCode();
1418 particle->Momentum(gammaj);
1419 //Info("GetLeadingPi0","Egammai %f, Egammaj %f",
1420 //gammai.Pt(),gammaj.Pt());
1422 ptl = (gammai+gammaj).Pt();
1423 phil = (gammai+gammaj).Phi();
1425 phil+=TMath::TwoPi();
1427 invmass = (gammai+gammaj).M();
1428 angle = gammaj.Angle(gammai.Vect());
1429 //Info("GetLeadingPi0","Angle %f", angle);
1430 e = (gammai+gammaj).E();
1433 (fListHistos->FindObject("AnglePairNoCut"))->
1436 (fListHistos->FindObject("InvMassPairNoCut"))->
1439 if((rat > fRatioMinCut) && (rat < fRatioMaxCut) &&
1440 ((phig-phil)>fPhiMinCut)&&((phig-phil)<fPhiMaxCut)){
1443 (fListHistos->FindObject("AnglePairLeadingCut"))->
1446 (fListHistos->FindObject("InvMassPairLeadingCut"))->
1449 if(IsAngleInWindow(angle,e)){
1451 (fListHistos->FindObject("AnglePairAngleCut"))->
1454 (fListHistos->FindObject("InvMassPairAngleCut"))->
1457 //Info("GetLeadingPi0","InvMass %f", invmass);
1458 if((invmass>fInvMassMinCut) && (invmass<fInvMassMaxCut)){
1460 (fListHistos->FindObject("InvMassPairAllCut"))->
1463 (fListHistos->FindObject("AnglePairAllCut"))->
1467 eta = particle->Eta() ;
1471 invmassf = invmass ;
1475 }//(invmass>0.125) && (invmass<0.145)
1476 }//gammaj.Angle(gammai.Vect())<0.04
1478 }//iprimary<jprimary
1486 (fListHistos->FindObject("InvMassPairLeading"))->
1489 (fListHistos->FindObject("AnglePairLeading"))->
1493 // printf("GetLeadingPi0: %f %f %f \n", pt, eta, phi) ;
1497 //____________________________________________________________________________
1498 void AliPHOSGammaJet::InitParameters()
1501 //Initialize the parameters of the analysis.
1503 fAngleMaxParam.Set(4) ;
1504 fAngleMaxParam.AddAt(0.4,0);//={0.4,-0.25,0.025,-2e-4};
1505 fAngleMaxParam.AddAt(-0.25,1) ;
1506 fAngleMaxParam.AddAt(0.025,2) ;
1507 fAngleMaxParam.AddAt(-2e-4,3) ;
1508 fAnyConeOrPt = kFALSE ;
1509 fOutputFileName = "GammaJet.root" ;
1511 fHIJINGFileName = "galice.root" ;
1513 fMinDistance = 3.6 ;
1515 fInvMassMaxCut = 0.15 ;
1516 fInvMassMinCut = 0.12 ;
1517 fOnlyCharged = kFALSE ;
1520 fPhiEMCALCut[0] = 60. *TMath::Pi()/180.;
1521 fPhiEMCALCut[1] = 180.*TMath::Pi()/180.;
1525 fNeutralPtCut = 0.5 ;
1526 fChargedPtCut = 0.5 ;
1527 fTPCCutsLikeEMCAL = kFALSE ;
1528 //Jet selection parameters
1530 fRatioMaxCut = 1.0 ;
1531 fRatioMinCut = 0.1 ;
1532 fJetRatioMaxCut = 1.2 ;
1533 fJetRatioMinCut = 0.8 ;
1534 fJetTPCRatioMaxCut = 1.2 ;
1535 fJetTPCRatioMinCut = 0.3 ;
1539 fESDTree = "esdTree" ;
1542 //Cut depending on gamma energy
1544 fPtJetSelectionCut = 20.; //For Low pt jets+BKG, another limits applyed
1545 //Reconstructed jet energy dependence parameters
1546 //e_jet = a1+e_gamma b2.
1547 //Index 0-> Pt>2 GeV r = 0.3; Index 1-> Pt>0.5 GeV r = 0.3
1548 fJetE1[0] = -5.75; fJetE1[1] = -4.1;
1549 fJetE2[0] = 1.005; fJetE2[1] = 1.05;
1551 //Reconstructed sigma of jet energy dependence parameters
1552 //s_jet = a1+e_gamma b2.
1553 //Index 0-> Pt>2 GeV r = 0.3; Index 1-> Pt>0.5 GeV r = 0.3
1554 fJetSigma1[0] = 2.65; fJetSigma1[1] = 2.75;
1555 fJetSigma2[0] = 0.0018; fJetSigma2[1] = 0.033;
1557 //Background mean energy and RMS
1558 //Index 0-> No BKG; Index 1-> BKG > 2 GeV;
1559 //Index 2-> (low pt jets)BKG > 0.5 GeV;
1560 //Index > 2, same for TPC conf
1561 fBkgMean[0] = 0.; fBkgMean[1] = 8.8 ; fBkgMean[2] = 69.5;
1562 fBkgMean[3] = 0.; fBkgMean[4] = 6.4; fBkgMean[5] = 48.6;
1563 fBkgRMS[0] = 0.; fBkgRMS[1] = 7.5; fBkgRMS[2] = 22.0;
1564 fBkgRMS[3] = 0.; fBkgRMS[4] = 5.4; fBkgRMS[5] = 13.2;
1566 //Factor x of min/max = E -+ x * sigma. Obtained after selecting the
1567 //limits for monoenergetic jets.
1568 //Index 0-> No BKG; Index 1-> BKG > 2 GeV;
1569 //Index 2-> (low pt jets) BKG > 0.5 GeV;
1570 //Index > 2, same for TPC conf
1572 fJetXMin1[0] =-0.69 ; fJetXMin1[1] = 0.39 ; fJetXMin1[2] =-0.88 ;
1573 fJetXMin1[3] =-2.0 ; fJetXMin1[4] =-0.442 ; fJetXMin1[5] =-1.1 ;
1574 fJetXMin2[0] = 0.066; fJetXMin2[1] = 0.038; fJetXMin2[2] = 0.034;
1575 fJetXMin2[3] = 0.25 ; fJetXMin2[4] = 0.113; fJetXMin2[5] = 0.077 ;
1576 fJetXMax1[0] =-3.8 ; fJetXMax1[1] =-0.76 ; fJetXMax1[2] =-3.6 ;
1577 fJetXMax1[3] =-2.7 ; fJetXMax1[4] =-1.21 ; fJetXMax1[5] =-3.7 ;
1578 fJetXMax2[0] =-0.012; fJetXMax2[1] =-0.022; fJetXMax2[2] = 0.016;
1579 fJetXMax2[3] =-0.024; fJetXMax2[4] =-0.008; fJetXMax2[5] = 0.027;
1582 //Photon fast reconstruction
1583 fResPara1 = 0.0255 ; // GeV
1584 fResPara2 = 0.0272 ;
1585 fResPara3 = 0.0129 ;
1587 fPosParaA = 0.096 ; // cm
1590 //Different cones and pt thresholds to construct the jet
1596 fCones[0] = 0.3 ; fNameCones[0] = "03" ;
1597 fPtThres[0] = 0.5 ; fNamePtThres[0] = "05" ;
1601 //__________________________________________________________________________-
1602 Bool_t AliPHOSGammaJet::IsAngleInWindow(const Float_t angle,const Float_t e) {
1603 //Check if the opening angle of the candidate pairs is inside
1604 //our selection windowd
1605 Bool_t result = kFALSE;
1606 Double_t mpi0 = 0.1349766;
1607 Double_t max = fAngleMaxParam.At(0)*TMath::Exp(fAngleMaxParam.At(1)*e)
1608 +fAngleMaxParam.At(2)+fAngleMaxParam.At(3)*e;
1609 Double_t arg = (e*e-2*mpi0*mpi0)/(e*e);
1610 Double_t min = 100. ;
1612 min = TMath::ACos(arg);
1614 if((angle<max)&&(angle>=min))
1620 //__________________________________________________________________________-
1621 Bool_t AliPHOSGammaJet::IsJetSelected(const Double_t ptg, const Double_t ptj,
1622 const TString type ){
1623 //Check if the energy of the reconstructed jet is within an energy window
1631 if(type == "TPC" && !fTPCCutsLikeEMCAL){
1632 iTPC = 3 ;//If(fTPCCutsLikeEMCAL) take jet energy cuts like EMCAL
1637 //Phythia alone, jets with pt_th > 0.2, r = 0.3
1638 par[0] = fJetE1[0]; par[1] = fJetE2[0];
1639 //Energy of the jet peak
1640 //e_jet = fJetE1[0]+fJetE2[0]*e_gamma, simulation fit
1641 par[2] = fJetSigma1[0]; par[3] = fJetSigma2[0];
1642 //Sigma of the jet peak
1643 //sigma_jet = fJetSigma1[0]+fJetSigma2[0]*e_gamma, simulation fit
1644 par[4] = fBkgMean[0 + iTPC]; par[5] = fBkgRMS[0 + iTPC];
1645 //Parameters reserved for HIJING bkg.
1646 xmax[0] = fJetXMax1[0 + iTPC]; xmax[1] = fJetXMax2[0 + iTPC];
1647 xmin[0] = fJetXMin1[0 + iTPC]; xmin[1] = fJetXMin2[0 + iTPC];
1648 //Factor that multiplies sigma to obtain the best limits,
1649 //by observation, of mono jet ratios (ptjet/ptg)
1650 //X_jet = fJetX1[0]+fJetX2[0]*e_gamma
1654 if(ptg > fPtJetSelectionCut){
1655 //Phythia +HIJING with pt_th > 2 GeV/c, r = 0.3
1656 par[0] = fJetE1[0]; par[1] = fJetE2[0];
1657 //Energy of the jet peak, same as in pp
1658 //e_jet = fJetE1[0]+fJetE2[0]*e_gamma, simulation fit
1659 par[2] = fJetSigma1[0]; par[3] = fJetSigma2[0];
1660 //Sigma of the jet peak, same as in pp
1661 //sigma_jet = fJetSigma1[0]+fJetSigma2[0]*e_gamma, simulation fit
1662 par[4] = fBkgMean[1 + iTPC]; par[5] = fBkgRMS[1 + iTPC];
1663 //Mean value and RMS of HIJING Bkg
1664 xmax[0] = fJetXMax1[1 + iTPC]; xmax[1] = fJetXMax2[1 + iTPC];
1665 xmin[0] = fJetXMin1[1 + iTPC]; xmin[1] = fJetXMin2[1 + iTPC];
1666 //Factor that multiplies sigma to obtain the best limits,
1667 //by observation, of mono jet ratios (ptjet/ptg) mixed with HIJING Bkg,
1668 //pt_th > 2 GeV, r = 0.3
1669 //X_jet = fJetX1[0]+fJetX2[0]*e_gamma
1673 //Phythia + HIJING with pt_th > 0.5 GeV/c, r = 0.3
1674 par[0] = fJetE1[1]; par[1] = fJetE2[1];
1675 //Energy of the jet peak, pt_th > 2 GeV/c, r = 0.3
1676 //e_jet = fJetE1[0]+fJetE2[0]*e_gamma, simulation fit
1677 par[2] = fJetSigma1[1]; par[3] = fJetSigma2[1];
1678 //Sigma of the jet peak, pt_th > 2 GeV/c, r = 0.3
1679 //sigma_jet = fJetSigma1[0]+fJetSigma2[0]*e_gamma, simulation fit
1680 par[4] = fBkgMean[2 + iTPC]; par[5] = fBkgRMS[2 + iTPC];
1681 //Mean value and RMS of HIJING Bkg in a 0.3 cone, pt > 2 GeV.
1682 xmax[0] = fJetXMax1[2 + iTPC]; xmax[1] = fJetXMax2[2 + iTPC];
1683 xmin[0] = fJetXMin1[2 + iTPC]; xmin[1] = fJetXMin2[2 + iTPC];
1684 //Factor that multiplies sigma to obtain the best limits,
1685 //by observation, of mono jet ratios (ptjet/ptg) mixed with HIJING Bkg,
1686 //pt_th > 2 GeV, r = 0.3
1687 //X_jet = fJetX1[0]+fJetX2[0]*e_gamma
1689 }//If low pt jet in bkg
1692 //Calculate minimum and maximum limits of the jet ratio.
1693 Double_t min = CalculateJetRatioLimit(ptg, par, xmin);
1694 Double_t max = CalculateJetRatioLimit(ptg, par, xmax);
1696 //Info("IsJetSeleted","%s : Limits min %f, max %f, ptg / ptj %f",
1697 // type.Data(),min,max,ptj/ptg);
1698 if(( min < ptj/ptg ) && ( max > ptj/ptg))
1705 //____________________________________________________________________________
1706 void AliPHOSGammaJet::List() const
1710 Info("List", "%d histograms found", fListHistos->GetEntries() ) ;
1711 TIter next(fListHistos) ;
1713 while ( (h = dynamic_cast<TH2F*>(next())) )
1714 Info("List", "%s", h->GetName()) ;
1717 //____________________________________________________________________________
1718 Double_t AliPHOSGammaJet::MakeEnergy(const Double_t energy)
1720 // Smears the energy according to the energy dependent energy resolution.
1721 // A gaussian distribution is assumed
1723 Double_t sigma = SigmaE(energy) ;
1724 return fRan.Gaus(energy, sigma) ;
1728 //____________________________________________________________________________
1729 void AliPHOSGammaJet::MakeHistos()
1731 // Create histograms to be saved in output file and
1732 // stores them in a TObjectArray
1734 fOutputFile = new TFile(fOutputFileName, "recreate") ;
1736 fListHistos = new TObjArray(10000) ;
1738 // Histos gamma pt vs leading pt
1740 TH1F * hPtSpectra = new TH1F
1741 ("PtSpectra","p_{T i} vs p_{T #gamma}",200,0,200);
1742 hPtSpectra->SetXTitle("p_{T} (GeV/c)");
1743 fListHistos->Add(hPtSpectra) ;
1745 //Histos ratio charged leading pt / gamma pt vs pt
1747 TH2F * hChargeRatio = new TH2F
1748 ("ChargeRatio","p_{T leading charge} /p_{T #gamma} vs p_{T #gamma}",
1750 hChargeRatio->SetYTitle("p_{T lead charge} /p_{T #gamma}");
1751 hChargeRatio->SetXTitle("p_{T #gamma} (GeV/c)");
1752 fListHistos->Add(hChargeRatio) ;
1754 TH2F * hTPCRatio = new TH2F
1755 ("TPCRatio","p_{T leading charge} /p_{T #gamma} vs p_{T #gamma}",
1757 hTPCRatio->SetYTitle("p_{T lead charge} /p_{T #gamma}");
1758 hTPCRatio->SetXTitle("p_{T #gamma} (GeV/c)");
1759 fListHistos->Add(hTPCRatio) ;
1762 TH2F * hPi0Ratio = new TH2F
1763 ("Pi0Ratio","p_{T leading #pi^{0}} /p_{T #gamma} vs p_{T #gamma}",
1765 hPi0Ratio->SetYTitle("p_{T lead #pi^{0}} /p_{T #gamma}");
1766 hPi0Ratio->SetXTitle("p_{T #gamma} (GeV/c)");
1767 fListHistos->Add(hPi0Ratio) ;
1769 TH2F * hPhiGamma = new TH2F
1770 ("PhiGamma","#phi_{#gamma}",200,0,120,200,0,7);
1771 hPhiGamma->SetYTitle("#phi");
1772 hPhiGamma->SetXTitle("p_{T #gamma} (GeV/c)");
1773 fListHistos->Add(hPhiGamma) ;
1775 TH2F * hEtaGamma = new TH2F
1776 ("EtaGamma","#phi_{#gamma}",200,0,120,200,-0.8,0.8);
1777 hEtaGamma->SetYTitle("#eta");
1778 hEtaGamma->SetXTitle("p_{T #gamma} (GeV/c)");
1779 fListHistos->Add(hEtaGamma) ;
1781 // //Jet reconstruction check
1782 // TH2F * hDeltaPhiJet = new TH2F
1783 // ("DeltaPhiJet","#phi_{jet} - #phi_{pyth jet} vs p_{T #gamma}",
1784 // 200,0,120,200,-1.5,1.5);
1785 // hDeltaPhiJet->SetYTitle("#Delta #phi");
1786 // hDeltaPhiJet->SetXTitle("p_{T #gamma} (GeV/c)");
1787 // fListHistos->Add(hDeltaPhiJet) ;
1789 // TH2F * hDeltaPhiTPCJet = new TH2F
1790 // ("DeltaPhiTPCJet","#phi_{jet TPC} - #phi_{pyth jet} vs p_{T #gamma}",
1791 // 200,0,120,200,-1.5,1.5);
1792 // hDeltaPhiTPCJet->SetYTitle("#Delta #phi");
1793 // hDeltaPhiTPCJet->SetXTitle("p_{T #gamma} (GeV/c)");
1794 // fListHistos->Add(hDeltaPhiTPCJet) ;
1796 // TH2F * hDeltaEtaJet = new TH2F
1797 // ("DeltaEtaJet","#phi_{jet} - #phi_{pyth jet} vs p_{T #gamma}",
1798 // 200,0,120,200,-1.5,1.5);
1799 // hDeltaEtaJet->SetYTitle("#Delta #phi");
1800 // hDeltaEtaJet->SetXTitle("p_{T #gamma} (GeV/c)");
1801 // fListHistos->Add(hDeltaEtaJet) ;
1803 // TH2F * hDeltaEtaTPCJet = new TH2F
1804 // ("DeltaEtaTPCJet","#phi_{jet TPC} - #phi_{pyth jet} vs p_{T #gamma}",
1805 // 200,0,120,200,-1.5,1.5);
1806 // hDeltaEtaTPCJet->SetYTitle("#Delta #phi");
1807 // hDeltaEtaTPCJet->SetXTitle("p_{T #gamma} (GeV/c)");
1808 // fListHistos->Add(hDeltaEtaTPCJet) ;
1810 // TH2F * hDeltaPtJet = new TH2F
1811 // ("DeltaPtJet","#phi_{jet} - #phi_{pyth jet} vs p_{T #gamma}",
1812 // 200,0,120,200,0.,100.);
1813 // hDeltaPtJet->SetYTitle("#Delta #phi");
1814 // hDeltaPtJet->SetXTitle("p_{T #gamma} (GeV/c)");
1815 // fListHistos->Add(hDeltaPtJet) ;
1817 // TH2F * hDeltaPtTPCJet = new TH2F
1818 // ("DeltaPtTPCJet","#phi_{jet TPC} - #phi_{pyth jet} vs p_{T #gamma}",
1819 // 200,0,120,200,0.,100.);
1820 // hDeltaPtTPCJet->SetYTitle("#Delta #phi");
1821 // hDeltaPtTPCJet->SetXTitle("p_{T #gamma} (GeV/c)");
1822 // fListHistos->Add(hDeltaPtTPCJet) ;
1825 TH2F * hDeltaPhiCharge = new TH2F
1826 ("DeltaPhiCharge","#phi_{#gamma} - #phi_{charge} vs p_{T #gamma}",
1827 200,0,120,200,0,6.4);
1828 hDeltaPhiCharge->SetYTitle("#Delta #phi");
1829 hDeltaPhiCharge->SetXTitle("p_{T #gamma} (GeV/c)");
1830 fListHistos->Add(hDeltaPhiCharge) ;
1832 TH2F * hDeltaPhiTPC = new TH2F
1833 ("DeltaPhiTPC","#phi_{#gamma} - #phi_{charge} vs p_{T #gamma}",
1834 200,0,120,200,0,6.4);
1835 hDeltaPhiTPC->SetYTitle("#Delta #phi");
1836 hDeltaPhiTPC->SetXTitle("p_{T #gamma} (GeV/c)");
1837 fListHistos->Add(hDeltaPhiTPC) ;
1838 TH2F * hDeltaPhiPi0 = new TH2F
1839 ("DeltaPhiPi0","#phi_{#gamma} - #phi_{ #pi^{0}} vs p_{T #gamma}",
1840 200,0,120,200,0,6.4);
1841 hDeltaPhiPi0->SetYTitle("#Delta #phi");
1842 hDeltaPhiPi0->SetXTitle("p_{T #gamma} (GeV/c)");
1843 fListHistos->Add(hDeltaPhiPi0) ;
1845 TH2F * hDeltaEtaCharge = new TH2F
1846 ("DeltaEtaCharge","#eta_{#gamma} - #eta_{charge} vs p_{T #gamma}",
1847 200,0,120,200,-2,2);
1848 hDeltaEtaCharge->SetYTitle("#Delta #eta");
1849 hDeltaEtaCharge->SetXTitle("p_{T #gamma} (GeV/c)");
1850 fListHistos->Add(hDeltaEtaCharge) ;
1852 TH2F * hDeltaEtaTPC = new TH2F
1853 ("DeltaEtaTPC","#eta_{#gamma} - #eta_{charge} vs p_{T #gamma}",
1854 200,0,120,200,-2,2);
1855 hDeltaEtaTPC->SetYTitle("#Delta #eta");
1856 hDeltaEtaTPC->SetXTitle("p_{T #gamma} (GeV/c)");
1857 fListHistos->Add(hDeltaEtaTPC) ;
1859 TH2F * hDeltaEtaPi0 = new TH2F
1860 ("DeltaEtaPi0","#eta_{#gamma} - #eta_{ #pi^{0}} vs p_{T #gamma}",
1861 200,0,120,200,-2,2);
1862 hDeltaEtaPi0->SetYTitle("#Delta #eta");
1863 hDeltaEtaPi0->SetXTitle("p_{T #gamma} (GeV/c)");
1864 fListHistos->Add(hDeltaEtaPi0) ;
1868 TH2F * hAnglePair = new TH2F
1870 "Angle between #pi^{0} #gamma pair vs p_{T #pi^{0}}",
1871 200,0,50,200,0,0.2);
1872 hAnglePair->SetYTitle("Angle (rad)");
1873 hAnglePair->SetXTitle("E_{ #pi^{0}} (GeV/c)");
1874 fListHistos->Add(hAnglePair) ;
1876 TH2F * hAnglePairAccepted = new TH2F
1877 ("AnglePairAccepted",
1878 "Angle between #pi^{0} #gamma pair vs p_{T #pi^{0}}, both #gamma in eta<0.7, inside window",
1879 200,0,50,200,0,0.2);
1880 hAnglePairAccepted->SetYTitle("Angle (rad)");
1881 hAnglePairAccepted->SetXTitle("E_{ #pi^{0}} (GeV/c)");
1882 fListHistos->Add(hAnglePairAccepted) ;
1884 TH2F * hAnglePairNoCut = new TH2F
1886 "Angle between all #gamma pair vs p_{T #pi^{0}}",200,0,50,200,0,0.2);
1887 hAnglePairNoCut->SetYTitle("Angle (rad)");
1888 hAnglePairNoCut->SetXTitle("E_{ #pi^{0}} (GeV/c)");
1889 fListHistos->Add(hAnglePairNoCut) ;
1891 TH2F * hAnglePairLeadingCut = new TH2F
1892 ("AnglePairLeadingCut",
1893 "Angle between all #gamma pair that have a good phi and pt vs p_{T #pi^{0}}",
1894 200,0,50,200,0,0.2);
1895 hAnglePairLeadingCut->SetYTitle("Angle (rad)");
1896 hAnglePairLeadingCut->SetXTitle("E_{ #pi^{0}} (GeV/c)");
1897 fListHistos->Add(hAnglePairLeadingCut) ;
1899 TH2F * hAnglePairAngleCut = new TH2F
1900 ("AnglePairAngleCut",
1901 "Angle between all #gamma pair (angle + leading cut) vs p_{T #pi^{0}}"
1902 ,200,0,50,200,0,0.2);
1903 hAnglePairAngleCut->SetYTitle("Angle (rad)");
1904 hAnglePairAngleCut->SetXTitle("E_{ #pi^{0}} (GeV/c)");
1905 fListHistos->Add(hAnglePairAngleCut) ;
1907 TH2F * hAnglePairAllCut = new TH2F
1909 "Angle between all #gamma pair (angle + inv mass cut+leading) vs p_{T #pi^{0}}"
1910 ,200,0,50,200,0,0.2);
1911 hAnglePairAllCut->SetYTitle("Angle (rad)");
1912 hAnglePairAllCut->SetXTitle("E_{ #pi^{0}} (GeV/c)");
1913 fListHistos->Add(hAnglePairAllCut) ;
1915 TH2F * hAnglePairLeading = new TH2F
1916 ("AnglePairLeading",
1917 "Angle between all #gamma pair finally selected vs p_{T #pi^{0}}",
1918 200,0,50,200,0,0.2);
1919 hAnglePairLeading->SetYTitle("Angle (rad)");
1920 hAnglePairLeading->SetXTitle("E_{ #pi^{0}} (GeV/c)");
1921 fListHistos->Add(hAnglePairLeading) ;
1924 TH2F * hInvMassPairNoCut = new TH2F
1925 ("InvMassPairNoCut","Invariant Mass of all #gamma pair vs p_{T #gamma}",
1926 120,0,120,360,0,0.5);
1927 hInvMassPairNoCut->SetYTitle("Invariant Mass (GeV/c^{2})");
1928 hInvMassPairNoCut->SetXTitle("p_{T #gamma} (GeV/c)");
1929 fListHistos->Add(hInvMassPairNoCut) ;
1931 TH2F * hInvMassPairLeadingCut = new TH2F
1932 ("InvMassPairLeadingCut",
1933 "Invariant Mass of #gamma pair (leading cuts) vs p_{T #gamma}",
1934 120,0,120,360,0,0.5);
1935 hInvMassPairLeadingCut->SetYTitle("Invariant Mass (GeV/c^{2})");
1936 hInvMassPairLeadingCut->SetXTitle("p_{T #gamma} (GeV/c)");
1937 fListHistos->Add(hInvMassPairLeadingCut) ;
1939 TH2F * hInvMassPairAngleCut = new TH2F
1940 ("InvMassPairAngleCut",
1941 "Invariant Mass of #gamma pair (angle cut) vs p_{T #gamma}",
1942 120,0,120,360,0,0.5);
1943 hInvMassPairAngleCut->SetYTitle("Invariant Mass (GeV/c^{2})");
1944 hInvMassPairAngleCut->SetXTitle("p_{T #gamma} (GeV/c)");
1945 fListHistos->Add(hInvMassPairAngleCut) ;
1948 TH2F * hInvMassPairAllCut = new TH2F
1949 ("InvMassPairAllCut",
1950 "Invariant Mass of #gamma pair (angle+invmass cut+leading) vs p_{T #gamma}",
1951 120,0,120,360,0,0.5);
1952 hInvMassPairAllCut->SetYTitle("Invariant Mass (GeV/c^{2})");
1953 hInvMassPairAllCut->SetXTitle("p_{T #gamma} (GeV/c)");
1954 fListHistos->Add(hInvMassPairAllCut) ;
1956 TH2F * hInvMassPairLeading = new TH2F
1957 ("InvMassPairLeading",
1958 "Invariant Mass of #gamma pair selected vs p_{T #gamma}",
1959 120,0,120,360,0,0.5);
1960 hInvMassPairLeading->SetYTitle("Invariant Mass (GeV/c^{2})");
1961 hInvMassPairLeading->SetXTitle("p_{T #gamma} (GeV/c)");
1962 fListHistos->Add(hInvMassPairLeading) ;
1967 TH1F * hNGamma = new TH1F("NGamma","Number of #gamma over PHOS",240,0,120);
1968 hNGamma->SetYTitle("N");
1969 hNGamma->SetXTitle("p_{T #gamma}(GeV/c)");
1970 fListHistos->Add(hNGamma) ;
1972 TH1F * hNBkg = new TH1F("NBkg","bkg multiplicity",9000,0,9000);
1973 hNBkg->SetYTitle("counts");
1974 hNBkg->SetXTitle("N");
1975 fListHistos->Add(hNBkg) ;
1977 TH2F * hNLeading = new TH2F
1978 ("NLeading","Accepted Jet Leading", 240,0,120,240,0,120);
1979 hNLeading->SetYTitle("p_{T charge} (GeV/c)");
1980 hNLeading->SetXTitle("p_{T #gamma}(GeV/c)");
1981 fListHistos->Add(hNLeading) ;
1984 TH1F * hN = new TH1F("NJet","Accepted jets",240,0,120);
1986 hN->SetXTitle("p_{T #gamma}(GeV/c)");
1987 fListHistos->Add(hN) ;
1990 //Ratios and Pt dist of reconstructed (not selected) jets
1992 TH2F * hJetRatio = new TH2F
1993 ("JetRatio","p_{T jet lead}/p_{T #gamma} vs p_{T #gamma}",
1994 240,0,120,200,0,10);
1995 hJetRatio->SetYTitle("p_{T jet lead #pi^{0}}/p_{T #gamma}");
1996 hJetRatio->SetXTitle("p_{T #gamma} (GeV/c)");
1997 fListHistos->Add(hJetRatio) ;
2000 TH2F * hJetPt = new TH2F
2001 ("JetPt", "p_{T jet lead} vs p_{T #gamma}",240,0,120,400,0,200);
2002 hJetPt->SetYTitle("p_{T jet lead #pi^{0}}/p_{T #gamma}");
2003 hJetPt->SetXTitle("p_{T #gamma} (GeV/c)");
2004 fListHistos->Add(hJetPt) ;
2009 TH2F * hBkgRatio = new TH2F
2010 ("BkgRatio","p_{T bkg lead}/p_{T #gamma} vs p_{T #gamma}",
2011 240,0,120,200,0,10);
2012 hBkgRatio->SetYTitle("p_{T bkg lead charge}/p_{T #gamma}");
2013 hBkgRatio->SetXTitle("p_{T #gamma} (GeV/c)");
2014 fListHistos->Add(hBkgRatio) ;
2017 TH2F * hBkgPt = new TH2F
2018 ("BkgPt","p_{T jet lead} vs p_{T #gamma}",240,0,120,400,0,200);
2019 hBkgPt->SetYTitle("p_{T jet lead charge}/p_{T #gamma}");
2020 hBkgPt->SetXTitle("p_{T #gamma} (GeV/c)");
2021 fListHistos->Add(hBkgPt) ;
2026 TH2F * hJetFragment =
2027 new TH2F("JetFragment","x = p_{T i charged}/p_{T #gamma}",
2028 240,0.,120.,1000,0.,1.2);
2029 hJetFragment->SetYTitle("x_{T}");
2030 hJetFragment->SetXTitle("p_{T #gamma}");
2031 fListHistos->Add(hJetFragment) ;
2033 TH2F * hBkgFragment = new TH2F
2034 ("BkgFragment","x = p_{T i charged}/p_{T #gamma}",
2035 240,0.,120.,1000,0.,1.2);
2036 hBkgFragment->SetYTitle("x_{T}");
2037 hBkgFragment->SetXTitle("p_{T #gamma}");
2038 fListHistos->Add(hBkgFragment) ;
2041 new TH2F("JetPtDist","x = p_{T i charged}",240,0.,120.,400,0.,200.);
2042 hJetPtDist->SetXTitle("p_{T #gamma} (GeV/c)");
2043 fListHistos->Add(hJetPtDist) ;
2045 TH2F * hBkgPtDist = new TH2F
2046 ("BkgPtDist","x = p_{T i charged}",240,0.,120.,400,0.,200.);
2047 hBkgPtDist->SetXTitle("p_{T #gamma} (GeV/c)");
2048 fListHistos->Add(hBkgPtDist) ;
2053 TH1F * hNBkgTPC = new TH1F
2054 ("NBkgTPC","TPC bkg multiplicity ",9000,0,9000);
2055 hNBkgTPC->SetYTitle("counts");
2056 hNBkgTPC->SetXTitle("N");
2057 fListHistos->Add(hNBkgTPC) ;
2059 TH2F * hNTPCLeading = new TH2F
2060 ("NTPCLeading","Accepted TPC jet leading",240,0,120,240,0,120);
2061 hNTPCLeading->SetYTitle("p_{T charge} (GeV/c)");
2062 hNTPCLeading->SetXTitle("p_{T #gamma}(GeV/c)");
2063 fListHistos->Add(hNTPCLeading) ;
2065 TH1F * hNTPC = new TH1F("NTPCJet","Number of TPC jets",240,0,120);
2066 hNTPC->SetYTitle("N");
2067 hNTPC->SetXTitle("p_{T #gamma}(GeV/c)");
2068 fListHistos->Add(hNTPC) ;
2070 TH2F * hJetTPCRatio = new TH2F
2071 ("JetTPCRatio", "p_{T jet lead TPC}/p_{T #gamma} vs p_{T #gamma}",
2072 240,0,120,200,0,10);
2073 hJetTPCRatio->SetYTitle("p_{T jet lead TPC}/p_{T #gamma}");
2074 hJetTPCRatio->SetXTitle("p_{T #gamma} (GeV/c)");
2075 fListHistos->Add(hJetTPCRatio) ;
2077 TH2F * hBkgTPCRatio = new TH2F
2078 ("BkgTPCRatio","p_{T bkg lead TPC}/p_{T #gamma} vs p_{T #gamma}",
2079 240,0,120,200,0,10);
2080 hBkgTPCRatio->SetYTitle("p_{T bkg lead TPC}/p_{T #gamma}");
2081 hBkgTPCRatio->SetXTitle("p_{T #gamma} (GeV/c)");
2082 fListHistos->Add(hBkgTPCRatio) ;
2084 TH2F * hJetTPCPt = new TH2F
2085 ("JetTPCPt", "p_{T jet lead TPC} vs p_{T #gamma}",240,0,120,400,0,200);
2086 hJetTPCPt->SetYTitle("p_{T jet lead TPC}/p_{T #gamma}");
2087 hJetTPCPt->SetXTitle("p_{T #gamma} (GeV/c)");
2088 fListHistos->Add(hJetTPCPt) ;
2090 TH2F * hBkgTPCPt = new TH2F
2091 ("BkgTPCPt", "p_{T bkg lead TPC} vs p_{T #gamma}",240,0,120,400,0,200);
2092 hBkgTPCPt->SetYTitle("p_{T bkg lead TPC}/p_{T #gamma}");
2093 hBkgTPCPt->SetXTitle("p_{T #gamma} (GeV/c)");
2094 fListHistos->Add(hBkgTPCPt) ;
2098 TH2F * hJetTPCFragment =
2099 new TH2F("JetTPCFragment","x = p_{T i charged}/p_{T #gamma}",
2100 240,0.,120.,1000,0.,1.2);
2101 hJetTPCFragment->SetYTitle("x_{T}");
2102 hJetTPCFragment->SetXTitle("p_{T #gamma}");
2103 fListHistos->Add(hJetTPCFragment) ;
2105 TH2F * hBkgTPCFragment = new TH2F
2106 ("BkgTPCFragment","x = p_{T i charged}/p_{T #gamma}",
2107 240,0.,120.,1000,0.,1.2);
2108 hBkgTPCFragment->SetYTitle("x_{T}");
2109 hBkgTPCFragment->SetXTitle("p_{T #gamma}");
2110 fListHistos->Add(hBkgTPCFragment) ;
2113 TH2F * hJetTPCPtDist = new TH2F("JetTPCPtDist",
2114 "x = p_{T i charged}",240,0.,120.,400,0.,200.);
2115 hJetTPCPtDist->SetXTitle("p_{T #gamma} (GeV/c)");
2116 fListHistos->Add(hJetTPCPtDist) ;
2118 TH2F * hBkgTPCPtDist = new TH2F
2119 ("BkgTPCPtDist","x = p_{T i charged}",240,0.,120.,400,0.,200.);
2120 hBkgTPCPtDist->SetXTitle("p_{T #gamma} (GeV/c)");
2121 fListHistos->Add(hBkgTPCPtDist) ;
2127 //If we want to study the jet for different cones and pt. Old version
2129 TH2F * hJetRatios[5][5];
2130 TH2F * hJetTPCRatios[5][5];
2132 TH2F * hJetPts[5][5];
2133 TH2F * hJetTPCPts[5][5];
2135 TH2F * hBkgRatios[5][5];
2136 TH2F * hBkgTPCRatios[5][5];
2138 TH2F * hBkgPts[5][5];
2139 TH2F * hBkgTPCPts[5][5];
2141 TH2F * hNLeadings[5][5];
2142 TH2F * hNTPCLeadings[5][5];
2145 TH1F * hNTPCs[5][5];
2147 TH1F * hNBkgs[5][5];
2148 TH1F * hNBkgTPCs[5][5];
2150 TH2F * hJetFragments[5][5];
2151 TH2F * hBkgFragments[5][5];
2152 TH2F * hJetPtDists[5][5];
2153 TH2F * hBkgPtDists[5][5];
2155 TH2F * hJetTPCFragments[5][5];
2156 TH2F * hBkgTPCFragments[5][5];
2157 TH2F * hJetTPCPtDists[5][5];
2158 TH2F * hBkgTPCPtDists[5][5];
2161 for(Int_t icone = 0; icone<fNCone; icone++){
2162 for(Int_t ipt = 0; ipt<fNPt;ipt++){
2167 hJetRatios[icone][ipt] = new TH2F
2168 ("JetRatioCone"+fNameCones[icone]+"Pt"+fNamePtThres[ipt],
2169 "p_{T jet lead #pi^{0}}/p_{T #gamma} vs p_{T #gamma}, cone ="
2170 +fNameCones[icone]+", pt>" +fNamePtThres[ipt]+" GeV/c",
2171 240,0,120,200,0,10);
2172 hJetRatios[icone][ipt]->
2173 SetYTitle("p_{T jet lead #pi^{0}}/p_{T #gamma}");
2174 hJetRatios[icone][ipt]->SetXTitle("p_{T #gamma} (GeV/c)");
2175 fListHistos->Add(hJetRatios[icone][ipt]) ;
2178 hJetPts[icone][ipt] = new TH2F
2179 ("JetPtCone"+fNameCones[icone]+"Pt"+fNamePtThres[ipt],
2180 "p_{T jet lead #pi^{0}}/p_{T #gamma} vs p_{T #gamma}, cone ="
2181 +fNameCones[icone]+", pt>" +fNamePtThres[ipt]+" GeV/c",
2182 240,0,120,400,0,200);
2183 hJetPts[icone][ipt]->
2184 SetYTitle("p_{T jet lead #pi^{0}}/p_{T #gamma}");
2185 hJetPts[icone][ipt]->SetXTitle("p_{T #gamma} (GeV/c)");
2186 fListHistos->Add(hJetPts[icone][ipt]) ;
2190 hBkgRatios[icone][ipt] = new TH2F
2191 ("BkgRatioCone"+fNameCones[icone]+"Pt"+fNamePtThres[ipt],
2192 "p_{T bkg lead #pi^{0}}/p_{T #gamma} vs p_{T #gamma}, cone ="
2193 +fNameCones[icone]+", pt>" +fNamePtThres[ipt]+" GeV/c",
2194 240,0,120,200,0,10);
2195 hBkgRatios[icone][ipt]->
2196 SetYTitle("p_{T bkg lead #pi^{0}}/p_{T #gamma}");
2197 hBkgRatios[icone][ipt]->SetXTitle("p_{T #gamma} (GeV/c)");
2198 fListHistos->Add(hBkgRatios[icone][ipt]) ;
2202 hBkgPts[icone][ipt] = new TH2F
2203 ("BkgPtCone"+fNameCones[icone]+"Pt"+fNamePtThres[ipt],
2204 "p_{T jet lead #pi^{0}}/p_{T #gamma} vs p_{T #gamma}, cone ="
2205 +fNameCones[icone]+", pt>" +fNamePtThres[ipt]+" GeV/c",
2206 240,0,120,400,0,200);
2207 hBkgPts[icone][ipt]->
2208 SetYTitle("p_{T jet lead #pi^{0}}/p_{T #gamma}");
2209 hBkgPts[icone][ipt]->SetXTitle("p_{T #gamma} (GeV/c)");
2210 fListHistos->Add(hBkgPts[icone][ipt]) ;
2215 hJetTPCRatios[icone][ipt] = new TH2F
2216 ("JetTPCRatioCone"+fNameCones[icone]+"Pt"+fNamePtThres[ipt],
2217 "p_{T jet lead TPC}/p_{T #gamma} vs p_{T #gamma}, cone ="
2218 +fNameCones[icone]+", pt>" +fNamePtThres[ipt]+" GeV/c",
2219 240,0,120,200,0,10);
2220 hJetTPCRatios[icone][ipt]->SetYTitle("p_{T jet lead TPC}/p_{T #gamma}");
2221 hJetTPCRatios[icone][ipt]->SetXTitle("p_{T #gamma} (GeV/c)");
2222 fListHistos->Add(hJetTPCRatios[icone][ipt]) ;
2224 hBkgTPCRatios[icone][ipt] = new TH2F
2225 ("BkgTPCRatioCone"+fNameCones[icone]+"Pt"+fNamePtThres[ipt],
2226 "p_{T bkg lead TPC}/p_{T #gamma} vs p_{T #gamma}, cone ="
2227 +fNameCones[icone]+", pt>" +fNamePtThres[ipt]+" GeV/c",
2228 240,0,120,200,0,10);
2229 hBkgTPCRatios[icone][ipt]->SetYTitle("p_{T bkg lead TPC}/p_{T #gamma}");
2230 hBkgTPCRatios[icone][ipt]->SetXTitle("p_{T #gamma} (GeV/c)");
2231 fListHistos->Add(hBkgTPCRatios[icone][ipt]) ;
2233 hJetTPCPts[icone][ipt] = new TH2F
2234 ("JetTPCPtCone"+fNameCones[icone]+"Pt"+fNamePtThres[ipt],
2235 "p_{T jet lead TPC}/p_{T #gamma} vs p_{T #gamma}, cone ="
2236 +fNameCones[icone]+", pt>" +fNamePtThres[ipt]+" GeV/c",
2237 240,0,120,400,0,200);
2238 hJetTPCPts[icone][ipt]->SetYTitle("p_{T jet lead TPC}/p_{T #gamma}");
2239 hJetTPCPts[icone][ipt]->SetXTitle("p_{T #gamma} (GeV/c)");
2240 fListHistos->Add(hJetTPCPts[icone][ipt]) ;
2242 hBkgTPCPts[icone][ipt] = new TH2F
2243 ("BkgTPCPtCone"+fNameCones[icone]+"Pt"+fNamePtThres[ipt],
2244 "p_{T bkg lead TPC}/p_{T #gamma} vs p_{T #gamma}, cone ="
2245 +fNameCones[icone]+", pt>" +fNamePtThres[ipt]+" GeV/c",
2246 240,0,120,400,0,200);
2247 hBkgTPCPts[icone][ipt]->SetYTitle("p_{T bkg lead TPC}/p_{T #gamma}");
2248 hBkgTPCPts[icone][ipt]->SetXTitle("p_{T #gamma} (GeV/c)");
2249 fListHistos->Add(hBkgTPCPts[icone][ipt]) ;
2254 hNBkgs[icone][ipt] = new TH1F
2255 ("NBkgCone"+fNameCones[icone]+"Pt"+fNamePtThres[ipt],
2256 "bkg multiplicity cone ="+fNameCones[icone]+", pt>"
2257 +fNamePtThres[ipt]+" GeV/c",9000,0,9000);
2258 hNBkgs[icone][ipt]->SetYTitle("counts");
2259 hNBkgs[icone][ipt]->SetXTitle("N");
2260 fListHistos->Add(hNBkgs[icone][ipt]) ;
2262 hNBkgTPCs[icone][ipt] = new TH1F
2263 ("NBkgTPCCone"+fNameCones[icone]+"Pt"+fNamePtThres[ipt],
2264 "bkg multiplicity cone ="+fNameCones[icone]+", pt>"
2265 +fNamePtThres[ipt]+" GeV/c",9000,0,9000);
2266 hNBkgTPCs[icone][ipt]->SetYTitle("counts");
2267 hNBkgTPCs[icone][ipt]->SetXTitle("N");
2268 fListHistos->Add(hNBkgTPCs[icone][ipt]) ;
2271 hNLeadings[icone][ipt] = new TH2F
2272 ("NLeadingCone"+fNameCones[icone]+"Pt"+fNamePtThres[ipt],
2273 "p_{T #gamma} vs p_{T #pi^{0}} cone ="+fNameCones[icone]+", pt>"
2274 +fNamePtThres[ipt]+" GeV/c",120,0,120,120,0,120);
2275 hNLeadings[icone][ipt]->SetYTitle("p_{T #pi^{0}}(GeV/c)");
2276 hNLeadings[icone][ipt]->SetXTitle("p_{T #gamma}(GeV/c)");
2277 fListHistos->Add(hNLeadings[icone][ipt]) ;
2279 hNs[icone][ipt] = new TH1F
2280 ("NJetCone"+fNameCones[icone]+"Pt"+fNamePtThres[ipt],
2281 "Number of neutral jets, cone ="+fNameCones[icone]+", pt>"
2282 +fNamePtThres[ipt]+" GeV/c",120,0,120);
2283 hNs[icone][ipt]->SetYTitle("N");
2284 hNs[icone][ipt]->SetXTitle("p_{T #gamma}(GeV/c)");
2285 fListHistos->Add(hNs[icone][ipt]) ;
2287 //Fragmentation Function
2288 hJetFragments[icone][ipt] = new TH2F
2289 ("JetFragmentCone"+fNameCones[icone]+"Pt"+fNamePtThres[ipt],
2290 "x_{T} = p_{T i}/p_{T #gamma}, cone ="+fNameCones[icone]+", pt>"
2291 +fNamePtThres[ipt]+" GeV/c",120,0.,120.,240,0.,1.2);
2292 hJetFragments[icone][ipt]->SetYTitle("x_{T}");
2293 hJetFragments[icone][ipt]->SetXTitle("p_{T #gamma}");
2294 fListHistos->Add(hJetFragments[icone][ipt]) ;
2296 hBkgFragments[icone][ipt] = new TH2F
2297 ("BkgFragmentCone"+fNameCones[icone]+"Pt"+fNamePtThres[ipt],
2298 "x_{T} = p_{T i}/p_{T #gamma}, cone ="+fNameCones[icone]+", pt>"
2299 +fNamePtThres[ipt]+" GeV/c",120,0.,120.,240,0.,1.2);
2300 hBkgFragments[icone][ipt]->SetYTitle("x_{T}");
2301 hBkgFragments[icone][ipt]->SetXTitle("p_{T #gamma}");
2302 fListHistos->Add(hBkgFragments[icone][ipt]) ;
2305 //Jet particle distribution
2307 hJetPtDists[icone][ipt] = new TH2F
2308 ("JetPtDistCone"+fNameCones[icone]+"Pt"+fNamePtThres[ipt],
2309 "p_{T i}, cone ="+fNameCones[icone]+", pt>" +fNamePtThres[ipt]+
2310 " GeV/c",120,0.,120.,120,0.,120.);
2311 hJetPtDists[icone][ipt]->SetXTitle("p_{T #gamma} (GeV/c)");
2312 fListHistos->Add(hJetPtDists[icone][ipt]) ;
2314 hBkgPtDists[icone][ipt] = new TH2F
2315 ("BkgPtDistCone"+fNameCones[icone]+"Pt"+fNamePtThres[ipt],
2316 "p_{T i}, cone ="+fNameCones[icone]+", pt>" +fNamePtThres[ipt]+
2317 " GeV/c",120,0.,120.,120,0.,120.);
2318 hBkgPtDists[icone][ipt]->SetXTitle("p_{T #gamma} (GeV/c)");
2319 fListHistos->Add(hBkgPtDists[icone][ipt]) ;
2323 hNTPCLeadings[icone][ipt] = new TH2F
2324 ("NTPCLeadingCone"+fNameCones[icone]+"Pt"+fNamePtThres[ipt],
2325 "p_{T #gamma} vs p_{T charge} cone ="+fNameCones[icone]+", pt>"
2326 +fNamePtThres[ipt]+" GeV/c",120,0,120,120,0,120);
2327 hNTPCLeadings[icone][ipt]->SetYTitle("p_{T charge} (GeV/c)");
2328 hNTPCLeadings[icone][ipt]->SetXTitle("p_{T #gamma}(GeV/c)");
2329 fListHistos->Add(hNTPCLeadings[icone][ipt]) ;
2331 hNTPCs[icone][ipt] = new TH1F
2332 ("NTPCJetCone"+fNameCones[icone]+"Pt"+fNamePtThres[ipt],
2333 "Number of charged jets, cone ="+fNameCones[icone]+", pt>"
2334 +fNamePtThres[ipt]+" GeV/c",120,0,120);
2335 hNTPCs[icone][ipt]->SetYTitle("N");
2336 hNTPCs[icone][ipt]->SetXTitle("p_{T #gamma}(GeV/c)");
2337 fListHistos->Add(hNTPCs[icone][ipt]) ;
2339 hJetTPCFragments[icone][ipt] = new TH2F
2340 ("JetTPCFragmentCone"+fNameCones[icone]+"Pt"+fNamePtThres[ipt],
2341 "x = p_{T i charged}/p_{T #gamma}, cone ="+fNameCones[icone]+", pt>"
2342 +fNamePtThres[ipt]+" GeV/c",120,0.,120.,240,0.,1.2);
2343 hJetTPCFragments[icone][ipt]->SetYTitle("x_{T}");
2344 hJetTPCFragments[icone][ipt]->SetXTitle("p_{T #gamma}");
2345 fListHistos->Add(hJetTPCFragments[icone][ipt]) ;
2347 hBkgTPCFragments[icone][ipt] = new TH2F
2348 ("BkgTPCFragmentCone"+fNameCones[icone]+"Pt"+fNamePtThres[ipt],
2349 "x = p_{T i charged}/p_{T #gamma}, cone ="+fNameCones[icone]+", pt>"
2350 +fNamePtThres[ipt]+" GeV/c",120,0.,120.,240,0.,1.2);
2351 hBkgTPCFragments[icone][ipt]->SetYTitle("x_{T}");
2352 hBkgTPCFragments[icone][ipt]->SetXTitle("p_{T #gamma}");
2353 fListHistos->Add(hBkgTPCFragments[icone][ipt]) ;
2355 hJetTPCPtDists[icone][ipt] = new TH2F
2356 ("JetTPCPtDistCone"+fNameCones[icone]+"Pt"+fNamePtThres[ipt],
2357 "x = p_{T i charged}, cone ="+fNameCones[icone]+", pt>"
2358 +fNamePtThres[ipt]+" GeV/c",120,0.,120.,120,0.,120.);
2359 hJetTPCPtDists[icone][ipt]->SetXTitle("p_{T #gamma} (GeV/c)");
2360 fListHistos->Add(hJetTPCPtDists[icone][ipt]) ;
2362 hBkgTPCPtDists[icone][ipt] = new TH2F
2363 ("BkgTPCPtDistCone"+fNameCones[icone]+"Pt"+fNamePtThres[ipt],
2364 "x = p_{T i charged}, cone ="+fNameCones[icone]+", pt>" +
2365 fNamePtThres[ipt]+" GeV/c",120,0.,120.,120,0.,120.);
2366 hBkgTPCPtDists[icone][ipt]->SetXTitle("p_{T #gamma} (GeV/c)");
2367 fListHistos->Add(hBkgTPCPtDists[icone][ipt]) ;
2371 }//If we want to study any cone or pt threshold
2375 //____________________________________________________________________________
2376 void AliPHOSGammaJet::MakeJet(TClonesArray * pl,
2377 Double_t ptg, Double_t phig,
2378 Double_t ptl, Double_t phil, Double_t etal,
2379 TString conf, TLorentzVector & jet)
2381 //Fill the jet with the particles around the leading particle with
2382 //R=fCone and pt_th = fPtThres. Calculate the energy of the jet and
2383 //check if we select it. Fill jet histograms
2384 Float_t ptcut = 0. ;
2386 if(ptg > fPtJetSelectionCut) ptcut = 2. ;
2390 TClonesArray * jetList = new TClonesArray("TParticle",1000);
2391 TClonesArray * bkgList = new TClonesArray("TParticle",1000);
2393 TLorentzVector bkg(0,0,0,0);
2394 TLorentzVector lv (0,0,0,0);
2396 Double_t ptjet = 0.0;
2397 Double_t ptbkg = 0.0;
2404 TParticle * particle = 0 ;
2406 while ( (particle = dynamic_cast<TParticle*>(next())) ) {
2412 SetJet(particle, b0, fCone, etal, phil) ;
2415 new((*jetList)[n0++]) TParticle(*particle) ;
2416 particle->Momentum(lv);
2417 if(particle->Pt() > ptcut ){
2419 ptjet+=particle->Pt();
2423 //Background around (phi_gamma-pi, eta_leading)
2424 SetJet(particle, b1, fCone,etal, phig) ;
2427 new((*bkgList)[n1++]) TParticle(*particle) ;
2428 if(particle->Pt() > ptcut ){
2430 ptbkg+=particle->Pt();
2438 if(strstr(fOptionGJ,"deb") || strstr(fOptionGJ,"deb all"))
2439 Info("MakeJet","Gamma pt %f, Jet pt %f, Bkg pt %f",ptg,ptjet,ptbkg);
2444 Double_t rat = ptjet/ptg ;
2445 Double_t ratbkg = ptbkg/ptg ;
2448 (fListHistos->FindObject("Jet"+conf+"Ratio"))->Fill(ptg,rat);
2450 (fListHistos->FindObject("Jet"+conf+"Pt")) ->Fill(ptg,ptjet);
2452 (fListHistos->FindObject("Bkg"+conf+"Ratio"))->Fill(ptg,ratbkg);
2454 (fListHistos->FindObject("Bkg"+conf+"Pt")) ->Fill(ptg,ptbkg);
2457 if(IsJetSelected(ptg,ptjet,conf) || fSelect){
2458 if(strstr(fOptionGJ,"deb") || strstr(fOptionGJ,"deb all"))
2459 Info("MakeJet","JetSelected");
2460 dynamic_cast<TH1F*>(fListHistos->FindObject("N"+conf+"Jet"))->
2462 dynamic_cast<TH2F*>(fListHistos->FindObject("N"+conf+"Leading"))
2464 FillJetHistos(jetList, ptg, conf, "Jet");
2465 FillJetHistos(bkgList, ptg, conf, "Bkg");
2472 //____________________________________________________________________________
2473 void AliPHOSGammaJet::MakeJetAnyConeOrPt(TClonesArray * pl, Double_t ptg,
2474 Double_t phig, Double_t ptl,
2475 Double_t phil, Double_t etal,
2478 //Fill the jet with the particles around the leading particle with
2479 //R=fCone(i) and pt_th = fPtThres(i). Calculate the energy of the jet and
2480 //check if we select it. Fill jet i histograms
2482 TClonesArray * jetList = new TClonesArray("TParticle",1000);
2483 TClonesArray * bkgList = new TClonesArray("TParticle",1000);
2485 Double_t ptjet = 0.0;
2486 Double_t ptbkg = 0.0;
2493 //Create as many jets as cones and pt thresholds are defined
2494 Double_t maxcut = fJetRatioMaxCut;
2495 Double_t mincut = fJetRatioMinCut;
2498 maxcut = fJetTPCRatioMaxCut;
2499 mincut = fJetTPCRatioMinCut;
2502 Double_t ratjet = 0;
2503 Double_t ratbkg = 0;
2505 for(Int_t icone = 0; icone<fNCone; icone++) {
2506 for(Int_t ipt = 0; ipt<fNPt;ipt++) {
2508 TString cone = fNameCones[icone] ;
2509 TString ptcut = fNamePtThres[ipt] ;
2512 TParticle * particle = 0 ;
2517 while ( (particle = dynamic_cast<TParticle*>(next())) ) {
2521 SetJet(particle, b0, fCones[icone],etal, phil) ;
2522 SetJet(particle, b1, fCones[icone],etal, phig) ;
2525 new((*jetList)[n0++]) TParticle(*particle) ;
2526 if(particle->Pt() > fPtThres[ipt] )
2527 ptjet+=particle->Pt();
2530 new((*bkgList)[n1++]) TParticle(*particle) ;
2531 if(particle->Pt() > fPtThres[ipt] )
2532 ptbkg+=particle->Pt();
2540 if(strstr(fOptionGJ,"deb")){
2541 Info("MakeJetAnyPt","cone %f, ptcut %f",fCones[icone],fPtThres[ipt]);
2542 Info("MakeJetAnyPt","pT: Gamma %f, Jet %f, Bkg %f",ptg,ptjet,ptbkg);
2545 ratjet = ptjet /ptg;
2549 (fListHistos->FindObject("Jet"+conf+"RatioCone"+cone+"Pt"+ptcut))
2552 (fListHistos->FindObject("Jet"+conf+"PtCone"+cone+"Pt"+ptcut))
2556 (fListHistos->FindObject("Bkg"+conf+"RatioCone"+cone+"Pt"+ptcut))
2560 (fListHistos->FindObject("Bkg"+conf+"PtCone"+cone+"Pt"+ptcut))
2565 if((ratjet < maxcut) && (ratjet > mincut)){
2568 (fListHistos->FindObject("N"+conf+"JetCone"+cone+"Pt"+ptcut))->
2571 (fListHistos->FindObject("N"+conf+"LeadingCone"+cone+"Pt"+ptcut))
2574 FillJetHistosAnyConeOrPt
2575 (jetList,ptg,conf,"Jet",fNameCones[icone],fNamePtThres[ipt]);
2576 FillJetHistosAnyConeOrPt
2577 (bkgList,ptg,conf,"Bkg",fNameCones[icone],fNamePtThres[ipt]);
2587 //____________________________________________________________________________
2588 void AliPHOSGammaJet::MakePhoton(TLorentzVector & particle)
2590 //Fast reconstruction for photons
2591 Double_t energy = particle.E() ;
2592 Double_t modenergy = MakeEnergy(energy) ;
2593 //Info("MakePhoton","Energy %f, Modif %f",energy,modenergy);
2595 // get the detected direction
2596 TVector3 pos = particle.Vect();
2598 TVector3 modpos = MakePosition(energy, pos) ;
2599 modpos *= modenergy / 460.;
2601 Float_t modtheta = modpos.Theta();
2602 Float_t modphi = modpos.Phi();
2604 // Set the modified 4-momentum of the reconstructed particle
2605 Float_t py = modenergy*TMath::Sin(modphi)*TMath::Sin(modtheta);
2606 Float_t px = modenergy*TMath::Cos(modphi)*TMath::Sin(modtheta);
2607 Float_t pz = modenergy*TMath::Cos(modtheta);
2609 particle.SetPxPyPzE(px,py,pz,modenergy);
2613 //____________________________________________________________________________
2614 TVector3 AliPHOSGammaJet::MakePosition(const Double_t energy, const TVector3 pos)
2616 // Smears the impact position according to the energy dependent position resolution
2617 // A gaussian position distribution is assumed
2621 Double_t sigma = SigmaP(energy) ;
2622 Double_t x = fRan.Gaus( pos.X(), sigma ) ;
2623 Double_t z = fRan.Gaus( pos.Z(), sigma ) ;
2624 Double_t y = pos.Y() ;
2630 // Info("MakePosition","Theta dif %f",pos.Theta()-newpos.Theta());
2631 // Info("MakePosition","Phi dif %f",pos.Phi()-newpos.Phi());
2635 //____________________________________________________________________________
2636 void AliPHOSGammaJet::Pi0Decay(Double_t mPi0, TLorentzVector &p0,
2637 TLorentzVector &p1, TLorentzVector &p2, Double_t &angle) {
2638 // Perform isotropic decay pi0 -> 2 photons
2639 // p0 is pi0 4-momentum (inut)
2640 // p1 and p2 are photon 4-momenta (output)
2641 // cout<<"Boost vector"<<endl;
2642 TVector3 b = p0.BoostVector();
2643 //cout<<"Parameters"<<endl;
2644 //Double_t mPi0 = p0.M();
2645 Double_t phi = TMath::TwoPi() * gRandom->Rndm();
2646 Double_t cosThe = 2 * gRandom->Rndm() - 1;
2647 Double_t cosPhi = TMath::Cos(phi);
2648 Double_t sinPhi = TMath::Sin(phi);
2649 Double_t sinThe = TMath::Sqrt(1-cosThe*cosThe);
2650 Double_t ePi0 = mPi0/2.;
2651 //cout<<"ePi0 "<<ePi0<<endl;
2652 //cout<<"Components"<<endl;
2653 p1.SetPx(+ePi0*cosPhi*sinThe);
2654 p1.SetPy(+ePi0*sinPhi*sinThe);
2655 p1.SetPz(+ePi0*cosThe);
2657 //cout<<"p1: "<<p1.Px()<<" "<<p1.Py()<<" "<<p1.Pz()<<" "<<p1.E()<<endl;
2658 //cout<<"p1 Mass: "<<p1.Px()*p1.Px()+p1.Py()*p1.Py()+p1.Pz()*p1.Pz()-p1.E()*p1.E()<<endl;
2659 p2.SetPx(-ePi0*cosPhi*sinThe);
2660 p2.SetPy(-ePi0*sinPhi*sinThe);
2661 p2.SetPz(-ePi0*cosThe);
2663 //cout<<"p2: "<<p2.Px()<<" "<<p2.Py()<<" "<<p2.Pz()<<" "<<p2.E()<<endl;
2664 //cout<<"p2 Mass: "<<p2.Px()*p2.Px()+p2.Py()*p2.Py()+p2.Pz()*p2.Pz()-p2.E()*p2.E()<<endl;
2665 //cout<<"Boost "<<b.X()<<" "<<b.Y()<<" "<<b.Z()<<endl;
2667 //cout<<"p1: "<<p1.Px()<<" "<<p1.Py()<<" "<<p1.Pz()<<" "<<p1.E()<<endl;
2669 //cout<<"p2: "<<p2.Px()<<" "<<p2.Py()<<" "<<p2.Pz()<<" "<<p2.E()<<endl;
2670 //cout<<"angle"<<endl;
2671 angle = p1.Angle(p2.Vect());
2672 //cout<<angle<<endl;
2675 //____________________________________________________________________________
2676 void AliPHOSGammaJet::Plot(TString what, Option_t * option) const
2678 //Plot some relevant histograms of the analysis
2679 TH2F * h = dynamic_cast<TH2F*>(fOutputFile->Get(what));
2683 else if (what == "all") {
2684 TCanvas * can = new TCanvas("GammaJet", "Gamma-Jet Study",10,40,1600,1200);
2686 can->Range(0,0,22,20);
2687 TPaveLabel *pl1 = new TPaveLabel(1,18,20,19.5,"Titre","br");
2688 pl1->SetFillColor(18);
2689 pl1->SetTextFont(32);
2690 pl1->SetTextColor(49);
2694 Float_t begx = -0.29, begy = 0., endx = 0., endy = 0.30 ;
2695 for (index = 0 ; index < fListHistos->GetEntries() ; index++) {
2696 TString name("pad") ;
2700 if (begx >= 1.0 || endx >= 1.0) {
2706 printf("%f %f %f %f \n", begx, begy, endx, endy) ;
2707 pad = new TPad(name,"This is a pad",begx,begy,endx,endy,33);
2710 fListHistos->At(index)->Draw(option) ;
2718 Info("Draw", "Histogram %s does not exist or unknown option", what.Data()) ;
2723 //____________________________________________________________________________
2724 void AliPHOSGammaJet::Print(const Option_t * opt) const
2727 //Print some relevant parameters set for the analysis
2731 Info("Print", "%s %s", GetName(), GetTitle() ) ;
2733 printf("Eta cut : %f\n", fEtaCut) ;
2734 printf("D phi max cut : %f\n", fPhiMaxCut) ;
2735 printf("D phi min cut : %f\n", fPhiMinCut) ;
2736 printf("Leading Ratio max cut : %f\n", fRatioMaxCut) ;
2737 printf("Leading Ratio min cut : %f\n", fRatioMinCut) ;
2738 printf("Jet Ratio max cut : %f\n", fJetRatioMaxCut) ;
2739 printf("Jet Ratio min cut : %f\n", fJetRatioMinCut) ;
2740 printf("Jet TPC Ratio max cut : %f\n", fJetTPCRatioMaxCut) ;
2741 printf("Jet TPC Ratio min cut : %f\n", fJetTPCRatioMinCut) ;
2742 printf("Fast recons : %d\n", fOptFast);
2743 printf("Inv Mass max cut : %f\n", fInvMassMaxCut) ;
2744 printf("Inv Mass min cut : %f\n", fInvMassMinCut) ;
2748 //__________________________________________________________________________
2749 TChain * AliPHOSGammaJet::ReadESDfromdisk(const UInt_t eventsToRead,
2750 const TString dirName,
2751 const TString esdTreeName,
2752 const char * pattern)
2754 // Reads ESDs from Disk
2757 AliInfo( Form("\nReading files in %s \nESD tree name is %s \nReading %d events",
2758 dirName.Data(), esdTreeName.Data(), eventsToRead) ) ;
2760 // create a TChain of all the files
2761 TChain * cESDTree = new TChain(esdTreeName) ;
2763 // read from the directory file until the require number of events are collected
2764 void * from = gSystem->OpenDirectory(dirName) ;
2766 AliError( Form("Directory %s does not exist") ) ;
2769 else{ // reading file names from directory
2770 const char * subdir ;
2771 // search all subdirectories witch matching pattern
2772 while( (subdir = gSystem->GetDirEntry(from)) &&
2773 (cESDTree->GetEntries() < eventsToRead)) {
2774 if ( strstr(subdir, pattern) != 0 ) {
2776 sprintf(file, "%s%s/AliESDs.root", dirName.Data(), subdir);
2777 AliInfo( Form("Adding %s\n", file) );
2778 cESDTree->Add(file) ;
2782 AliInfo( Form(" %d events read", cESDTree->GetEntriesFast()) ) ;
2785 } // reading file names from directory
2789 //__________________________________________________________________________
2790 TChain * AliPHOSGammaJet::ReadESD(const UInt_t eventsToRead,
2791 const TString dirName,
2792 const TString esdTreeName,
2793 const char * pattern)
2795 // Read AliESDs files and return a Chain of events
2797 if ( dirName == "" ) {
2798 AliError("Give the name of the DIR where to find files") ;
2801 if ( esdTreeName == "" )
2802 return ReadESDfromdisk(eventsToRead, dirName) ;
2803 else if ( strcmp(pattern, "") == 0 )
2804 return ReadESDfromdisk(eventsToRead, dirName, esdTreeName) ;
2806 return ReadESDfromdisk(eventsToRead, dirName, esdTreeName, pattern) ;
2809 //___________________________________________________________________
2810 void AliPHOSGammaJet::SetJet(TParticle * part, Bool_t & b, Float_t cone,
2811 Double_t eta, Double_t phi)
2814 //Check if the particle is inside the cone defined by the leading particle
2817 if(phi > TMath::TwoPi())
2818 phi-=TMath::TwoPi();
2820 phi+=TMath::TwoPi();
2822 Double_t rad = 10000 + cone;
2824 if(TMath::Abs(part->Phi()-phi) <= (TMath::TwoPi() - cone))
2825 rad = TMath::Sqrt(TMath::Power(part->Eta()-eta,2)+
2826 TMath::Power(part->Phi()-phi,2));
2828 if(part->Phi()-phi > TMath::TwoPi() - cone)
2829 rad = TMath::Sqrt(TMath::Power(part->Eta()-eta,2)+
2830 TMath::Power((part->Phi()-TMath::TwoPi())-phi,2));
2831 if(part->Phi()-phi < -(TMath::TwoPi() - cone))
2832 rad = TMath::Sqrt(TMath::Power(part->Eta()-eta,2)+
2833 TMath::Power((part->Phi()+TMath::TwoPi())-phi,2));
2841 //____________________________________________________________________________
2842 Double_t AliPHOSGammaJet::SigmaE(Double_t energy)
2844 // Calculates the energy dependent energy resolution
2848 rv = TMath::Sqrt( TMath::Power(fResPara1/energy, 2)
2849 + TMath::Power(fResPara2/TMath::Sqrt(energy), 2)
2850 + TMath::Power(fResPara3, 2) ) ;
2852 return rv * energy ;
2855 //____________________________________________________________________________
2856 Double_t AliPHOSGammaJet::SigmaP(Double_t energy)
2858 // Calculates the energy dependent position resolution
2860 Double_t sigma = TMath::Sqrt(TMath::Power(fPosParaA,2) +
2861 TMath::Power(fPosParaB,2) / energy) ;
2864 return sigma ; // in cm