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.13 2006/04/26 07:32:37 hristov
21 * Coding conventions, clean-up and related changes
23 * Revision 1.12 2006/03/10 13:23:36 hristov
24 * Using AliESDCaloCluster instead of AliESDtrack
26 * Revision 1.11 2006/01/31 20:30:52 hristov
29 * Revision 1.10 2006/01/23 18:04:08 hristov
30 * Removing meaningless const
32 * Revision 1.9 2006/01/12 16:23:26 schutz
33 * ESD is properly read with methods of macros/AliReadESD.C copied in it
35 * Revision 1.8 2005/12/20 07:08:32 schutz
36 * corrected error in call AliReadESD
38 * Revision 1.6 2005/05/28 14:19:04 schutz
39 * Compilation warnings fixed by T.P.
43 //_________________________________________________________________________
44 // Class for the analysis of gamma-jet correlations
45 // Basically it seaches for a prompt photon in the PHOS acceptance,
46 // if so we construct a jet around the highest pt particle in the opposite
47 // side in azimuth, inside the TPC and EMCAL acceptances. First the leading
48 // particle and then the jet have to fullfill several conditions
49 // (energy, direction ..) to be accepted. Then the fragmentation function
50 // of this jet is constructed
52 //*-- Author: Gustavo Conesa & Yves Schutz (IFIC, CERN)
53 //////////////////////////////////////////////////////////////////////////////
56 // --- ROOT system ---
59 #include <TParticle.h>
61 #include <TPaveLabel.h>
65 #include "AliPHOSGammaJet.h"
66 #include "AliPHOSGetter.h"
67 #include "AliPHOSGeometry.h"
68 #include "AliPHOSFastGlobalReconstruction.h"
70 #include "AliESDtrack.h"
71 #include "AliESDCaloCluster.h"
72 #include "Riostream.h"
75 ClassImp(AliPHOSGammaJet)
77 //____________________________________________________________________________
78 AliPHOSGammaJet::AliPHOSGammaJet() :
79 TTask(), fAnyConeOrPt(0), fOptionGJ(""),
80 fOutputFile(new TFile(gDirectory->GetName())),
81 fOutputFileName(gDirectory->GetName()),
82 fInputFileName(gDirectory->GetName()),
83 fHIJINGFileName(gDirectory->GetName()),
84 fHIJING(0), fESDdata(0), fEtaCut(0.),
85 fOnlyCharged(0), fPhiMaxCut(0.),
86 fPhiMinCut(0.), fPtCut(0.),
87 fNeutralPtCut(0.), fChargedPtCut(0.),
88 fInvMassMaxCut(0.), fInvMassMinCut(0.),
89 fMinDistance(0.), fRatioMaxCut(0.), fRatioMinCut(0.),
90 fTPCCutsLikeEMCAL(0), fDirName(""), fESDTree(""),
91 fPattern(""), fJetTPCRatioMaxCut(0.),
92 fJetTPCRatioMinCut(0.), fJetRatioMaxCut(0.),
93 fJetRatioMinCut(0.), fNEvent(0), fNCone(0),
94 fNPt(0), fCone(0), fPtThreshold(0),
95 fPtJetSelectionCut(0.0),
96 fListHistos(new TObjArray(100)),
97 fFastRec(0), fOptFast(0),
98 fRan(0), fResPara1(0.), fResPara2(0.), fResPara3(0.),
99 fPosParaA(0.), fPosParaB(0.), fAngleMaxParam(), fSelect(0)
102 fAngleMaxParam.Set(4) ;
103 fAngleMaxParam.Reset(0.);
105 for(Int_t i = 0; i<10; i++){
109 fNamePtThres[i] = "" ;
120 fJetSigma1[i] = 0.0 ;
121 fJetSigma2[i] = 0.0 ;
122 fPhiEMCALCut[i] = 0.0 ;
127 TList * list = gDirectory->GetListOfKeys() ;
131 for (index = 0 ; index < list->GetSize()-1 ; index++) {
132 //-1 to avoid GammaJet Task
133 h = dynamic_cast<TH2F*>(gDirectory->Get(list->At(index)->GetName())) ;
134 fListHistos->Add(h) ;
139 //____________________________________________________________________________
140 AliPHOSGammaJet::AliPHOSGammaJet(const TString inputfilename) :
141 TTask("GammaJet","Analysis of gamma-jet correlations"),
142 fAnyConeOrPt(0), fOptionGJ(),
147 fHIJING(0), fESDdata(0), fEtaCut(0.),
148 fOnlyCharged(0), fPhiMaxCut(0.),
149 fPhiMinCut(0.), fPtCut(0.),
150 fNeutralPtCut(0.), fChargedPtCut(0.),
151 fInvMassMaxCut(0.), fInvMassMinCut(0.),
152 fMinDistance(0.), fRatioMaxCut(0.), fRatioMinCut(0.),
153 fTPCCutsLikeEMCAL(0), fDirName(), fESDTree(),
154 fPattern(), fJetTPCRatioMaxCut(0.),
155 fJetTPCRatioMinCut(0.), fJetRatioMaxCut(0.),
156 fJetRatioMinCut(0.), fNEvent(0), fNCone(0),
157 fNPt(0), fCone(0), fPtThreshold(0),
158 fPtJetSelectionCut(0.0),
160 fFastRec(0), fOptFast(0),
161 fRan(0), fResPara1(0.), fResPara2(0.), fResPara3(0.),
162 fPosParaA(0.), fPosParaB(0.), fAngleMaxParam(), fSelect(0)
166 fInputFileName = inputfilename;
167 fFastRec = new AliPHOSFastGlobalReconstruction(fInputFileName);
168 AliPHOSGetter * gime = AliPHOSGetter::Instance(fInputFileName) ;
169 fNEvent = gime->MaxEvent();
173 //____________________________________________________________________________
174 AliPHOSGammaJet::AliPHOSGammaJet(const AliPHOSGammaJet & gj) :
176 fAnyConeOrPt(gj.fAnyConeOrPt), fOptionGJ(gj.fOptionGJ),
177 fOutputFile(gj.fOutputFile),
178 fOutputFileName(gj.fOutputFileName),
179 fInputFileName(gj.fInputFileName),
180 fHIJINGFileName(gj.fHIJINGFileName),
181 fHIJING(gj.fHIJING), fESDdata(gj.fESDdata), fEtaCut(gj.fEtaCut),
182 fOnlyCharged(gj.fOnlyCharged), fPhiMaxCut(gj.fPhiMaxCut),
183 fPhiMinCut(gj.fPhiMinCut), fPtCut(gj.fPtCut),
184 fNeutralPtCut(gj.fNeutralPtCut), fChargedPtCut(gj.fChargedPtCut),
185 fInvMassMaxCut(gj.fInvMassMaxCut), fInvMassMinCut(gj.fInvMassMinCut),
186 fMinDistance(gj.fMinDistance), fRatioMaxCut(gj.fRatioMaxCut),
187 fRatioMinCut(gj.fRatioMinCut), fTPCCutsLikeEMCAL(gj.fTPCCutsLikeEMCAL),
188 fDirName(gj.fDirName), fESDTree(gj.fESDTree),
189 fPattern(gj.fPattern), fJetTPCRatioMaxCut(gj.fJetTPCRatioMaxCut),
190 fJetTPCRatioMinCut(gj.fJetTPCRatioMinCut), fJetRatioMaxCut(gj.fJetRatioMaxCut),
191 fJetRatioMinCut(gj.fJetRatioMinCut), fNEvent(gj.fNEvent), fNCone(gj.fNCone),
192 fNPt(gj.fNPt), fCone(gj.fCone), fPtThreshold(gj.fPtThreshold),
193 fPtJetSelectionCut(gj.fPtJetSelectionCut),
194 fListHistos(0),//?????
195 fFastRec(gj.fFastRec), fOptFast(gj.fOptFast),
197 fResPara1(gj.fResPara1), fResPara2(gj.fResPara2), fResPara3(gj.fResPara3),
198 fPosParaA(gj.fPosParaA), fPosParaB(gj.fPosParaB),
199 fAngleMaxParam(gj.fAngleMaxParam), fSelect(gj.fSelect)
202 SetName (gj.GetName()) ;
203 SetTitle(gj.GetTitle()) ;
205 for(Int_t i = 0; i<10; i++){
206 fCones[i] = gj.fCones[i] ;
207 fNameCones[i] = gj.fNameCones[i] ;
208 fPtThres[i] = gj.fPtThres[i] ;
209 fNamePtThres[i] = gj.fNamePtThres[i] ;
211 fJetXMin1[i] = gj.fJetXMin1[i] ;
212 fJetXMin2[i] = gj.fJetXMin2[i] ;
213 fJetXMax1[i] = gj.fJetXMax1[i] ;
214 fJetXMax2[i] = gj.fJetXMax2[i] ;
215 fBkgMean[i] = gj.fBkgMean[i] ;
216 fBkgRMS[i] = gj.fBkgRMS[i] ;
218 fJetE1[i] = gj.fJetE1[i] ;
219 fJetE2[i] = gj.fJetE2[i] ;
220 fJetSigma1[i] = gj.fJetSigma1[i] ;
221 fJetSigma2[i] = gj.fJetSigma2[i] ;
222 fPhiEMCALCut[i] = gj.fPhiEMCALCut[i] ;
228 //____________________________________________________________________________
229 AliPHOSGammaJet::~AliPHOSGammaJet()
231 fOutputFile->Close() ;
234 //____________________________________________________________________________
235 void AliPHOSGammaJet::AddHIJINGToList(Int_t iEvent, TClonesArray * particleList,
238 TClonesArray * plNePHOS)
241 // List of particles copied to a root file.
243 // Char_t sb[] = "bgrd/";
244 // // cout<<sb<<endl;
246 // Char_t sf[] = "/list.root";
247 // // cout<<sf<<endl;
248 // sprintf(si,"%d",iEvent);
250 // // cout<<si<<endl;
252 // // cout<<si<<endl;
253 // TFile * f = TFile::Open(sb) ;
254 // //cout<<f->GetName()<<endl;
257 sprintf(fi,"bgrd/%d/list.root",iEvent);
258 TFile * f = TFile::Open(fi) ;
259 //cout<<f->GetName()<<endl;
261 TParticle *particle = new TParticle();
262 TTree *t = (TTree*) f->Get("T");
263 TBranch *branch = t->GetBranch("primaries");
264 branch->SetAddress(&particle);
266 Int_t index = particleList->GetEntries() ;
267 Int_t indexNe = plNe->GetEntries() ;
268 Int_t indexCh = plCh->GetEntries() ;
269 Int_t indexNePHOS = plNePHOS->GetEntries() ;
270 Double_t charge = 0.;
271 Int_t iParticle = 0 ;
273 Double_t x = 0., z = 0.;
274 // cout<<"bkg entries "<<t->GetEntries()<<endl;
276 AliPHOSGetter * gime = AliPHOSGetter::Instance() ;
277 const AliPHOSGeometry * geom = gime->PHOSGeometry() ;
280 for (iParticle=0 ; iParticle < t->GetEntries() ; iParticle++) {
281 t->GetEvent(iParticle) ;
287 charge = TDatabasePDG::Instance()
288 ->GetParticle(particle->GetPdgCode())->Charge();
290 if((charge != 0) && (particle->Pt() > fChargedPtCut)){
291 if(TMath::Abs(particle->Eta())<fEtaCut){
292 new((*particleList)[index]) TParticle(*particle) ;
293 (dynamic_cast<TParticle*>(particleList->At(index)))
297 new((*plCh)[indexCh]) TParticle(*particle) ;
298 (dynamic_cast<TParticle*>(plCh->At(indexCh)))->SetStatusCode(0) ;
301 if(strstr(fOptionGJ,"deb all")||strstr(fOptionGJ,"deb")){
302 dynamic_cast<TH1F*>(fListHistos->FindObject("PtSpectra"))->
303 Fill(particle->Pt());
306 else if((charge == 0) && (particle->Pt() > fNeutralPtCut) ){
307 geom->ImpactOnEmc(particle->Theta(),particle->Phi(), m,z,x);
311 if(strstr(fOptionGJ,"deb all")|| strstr(fOptionGJ,"deb"))
312 dynamic_cast<TH1F*>(fListHistos->FindObject("PtSpectra"))->
313 Fill(particle->Pt());
315 new((*plNePHOS)[indexNePHOS]) TParticle(*particle) ;
316 (dynamic_cast<TParticle*>(plNePHOS->At(indexNePHOS)))->SetStatusCode(0) ;
320 if((particle->Phi()>fPhiEMCALCut[0]) &&
321 (particle->Phi()<fPhiEMCALCut[1]) && m == 0)
323 if(strstr(fOptionGJ,"deb all")|| strstr(fOptionGJ,"deb"))
324 dynamic_cast<TH1F*>(fListHistos->FindObject("PtSpectra"))->
325 Fill(particle->Pt());
327 new((*particleList)[index]) TParticle(*particle) ;
328 (dynamic_cast<TParticle*>(particleList->At(index)))
331 new((*plNe)[indexNe]) TParticle(*particle) ;
332 (dynamic_cast<TParticle*>(plNe->At(indexNe)))->SetStatusCode(0) ;
341 Double_t mass = TDatabasePDG::Instance()->GetParticle(111)->Mass();
342 TLorentzVector pPi0, pGamma1, pGamma2 ;
343 Double_t angle = 0, cellDistance = 0.;
346 // fFastRec = new AliPHOSFastGlobalReconstruction(fHIJINGFileName);
347 // fFastRec->FastReconstruction(iiEvent);
349 for (iParticle=0 ; iParticle < t->GetEntries() ; iParticle++) {
350 t->GetEvent(iParticle) ;
354 charge = TDatabasePDG::Instance()
355 ->GetParticle(particle->GetPdgCode())->Charge();
357 if((charge != 0) && (particle->Pt() > fChargedPtCut)
358 && (TMath::Abs(particle->Eta())<fEtaCut)){
360 new((*particleList)[index]) TParticle(*particle) ;
361 (dynamic_cast<TParticle*>(particleList->At(index)))
365 new((*plCh)[indexCh]) TParticle(*particle) ;
366 (dynamic_cast<TParticle*>(plCh->At(indexCh)))->SetStatusCode(0) ;
369 else if(charge == 0){
370 geom->ImpactOnEmc(particle->Theta(),particle->Phi(), m,z,x);
371 if((particle->GetPdgCode() != 111) && (particle->Pt() > fNeutralPtCut) &&
372 (TMath::Abs(particle->Eta())<fEtaCut) ){
373 TLorentzVector part(particle->Px(),particle->Py(),
374 particle->Pz(),particle->Energy());
376 if(part.Pt() > fNeutralPtCut){
378 if(particle->Phi()>fPhiEMCALCut[0] &&
379 particle->Phi()<fPhiEMCALCut[1] && m == 0)
381 new((*particleList)[index]) TParticle(*particle) ;
382 (dynamic_cast<TParticle*>(particleList->At(index)))
384 (dynamic_cast<TParticle*>(particleList->At(index)))
385 ->SetMomentum(part.Px(),part.Py(),part.Pz(),part.E());
389 new((*plNe)[indexNe]) TParticle(*particle) ;
390 (dynamic_cast<TParticle*>(plNe->At(indexNe)))
391 ->SetMomentum(part.Px(),part.Py(),part.Pz(),part.E());
392 (dynamic_cast<TParticle*>(plNe->At(indexNe)))->SetStatusCode(0) ;
397 new((*plNePHOS)[indexNePHOS]) TParticle(*particle) ;
398 (dynamic_cast<TParticle*>(plNePHOS->At(indexNePHOS)))
399 ->SetMomentum(part.Px(),part.Py(),part.Pz(),part.E());
400 (dynamic_cast<TParticle*>(plNePHOS->At(indexNePHOS)))->SetStatusCode(0) ;
405 if((particle->GetPdgCode() == 111) && (particle->Pt() > fNeutralPtCut) &&
406 (TMath::Abs(particle->Eta())<fEtaCut+1))
409 pPi0.SetPxPyPzE(particle->Px(),particle->Py(),particle->Pz(),
414 Pi0Decay(mass,pPi0,pGamma1,pGamma2,angle);
415 //Check if decay photons are too close for PHOS
416 cellDistance = angle*460; //cm
417 if (cellDistance < fMinDistance) {
418 if(strstr(fOptionGJ,"deb all")|| strstr(fOptionGJ,"deb"))
419 dynamic_cast<TH1F*>(fListHistos->FindObject("PtSpectra"))->
420 Fill(particle->Pt());
422 //Pi0 inside phi EMCAL acceptance
424 TLorentzVector part(particle->Px(),particle->Py(),
425 particle->Pz(),particle->Energy());
427 if(part.Pt() > fNeutralPtCut){
428 if(particle->Phi()>fPhiEMCALCut[0] &&
429 particle->Phi()<fPhiEMCALCut[1] && m == 0){
431 new((*particleList)[index]) TParticle(*particle) ;
432 (dynamic_cast<TParticle*>(particleList->At(index)))->SetStatusCode(0) ;
433 (dynamic_cast<TParticle*>(particleList->At(index)))
434 ->SetMomentum(part.Px(),part.Py(),part.Pz(),part.E());
437 new((*plNe)[indexNe]) TParticle(*particle) ;
438 (dynamic_cast<TParticle*>(plNe->At(indexNe))) ->SetStatusCode(0) ;
439 (dynamic_cast<TParticle*>(plNe->At(indexNe)))
440 ->SetMomentum(part.Px(),part.Py(),part.Pz(),part.E());
444 if(strstr(fOptionGJ,"deb all")|| strstr(fOptionGJ,"deb"))
445 dynamic_cast<TH1F*>(fListHistos->FindObject("PtSpectra"))->
446 Fill(particle->Pt());
447 new((*plNePHOS)[indexNePHOS]) TParticle(*particle) ;
448 (dynamic_cast<TParticle*>(plNePHOS->At(indexNePHOS))) ->SetStatusCode(0) ;
449 (dynamic_cast<TParticle*>(plNePHOS->At(indexNePHOS)))
450 ->SetMomentum(part.Px(),part.Py(),part.Pz(),part.E());
457 dynamic_cast<TH2F*>(fListHistos->FindObject("AnglePair"))
458 ->Fill(pPi0.E(),angle);
461 if(pGamma1.Pt() > 0. && TMath::Abs(pGamma1.Eta())<fEtaCut){
465 if(pGamma1.Pt() > fNeutralPtCut ){
467 TParticle * photon1 =
468 new TParticle(22,1,0,0,0,0,pGamma1.Px(),pGamma1.Py(),
469 pGamma1.Pz(),pGamma1.E(),0,0,0,0);
470 geom->ImpactOnEmc(photon1->Theta(),photon1->Phi(), m,z,x);
471 if( photon1->Phi()>fPhiEMCALCut[0] && photon1->Phi()<fPhiEMCALCut[1]
473 if(strstr(fOptionGJ,"deb all") || strstr(fOptionGJ,"deb"))
474 dynamic_cast<TH1F*>(fListHistos->FindObject("PtSpectra"))->
476 new((*particleList)[index]) TParticle(*photon1) ;
477 (dynamic_cast<TParticle*>(particleList->At(index)))->SetStatusCode(0) ;
480 new((*plNe)[indexNe]) TParticle(*photon1) ;
481 (dynamic_cast<TParticle*>(plNe->At(indexNe))) ->SetStatusCode(0) ;
486 if(strstr(fOptionGJ,"deb all") || strstr(fOptionGJ,"deb"))
487 dynamic_cast<TH1F*>(fListHistos->FindObject("PtSpectra"))->
489 new((*plNePHOS)[indexNePHOS]) TParticle(*photon1) ;
490 (dynamic_cast<TParticle*>(plNePHOS->At(indexNePHOS)))->SetStatusCode(0) ;
496 if(pGamma2.Pt() > 0. && TMath::Abs(pGamma2.Eta())<fEtaCut){
499 if(pGamma2.Pt() > fNeutralPtCut){
501 TParticle * photon2 =
502 new TParticle(22,1,0,0,0,0,pGamma2.Px(), pGamma2.Py(),
503 pGamma2.Pz(),pGamma2.E(),0,0,0,0);
504 geom->ImpactOnEmc(photon2->Theta(),photon2->Phi(), m,z,x);
505 if(photon2->Phi()>fPhiEMCALCut[0] &&
506 photon2->Phi()<fPhiEMCALCut[1] && m == 0){
507 if(strstr(fOptionGJ,"deb all") || strstr(fOptionGJ,"deb"))
508 dynamic_cast<TH1F*>(fListHistos->FindObject("PtSpectra"))->
510 new((*particleList)[index]) TParticle(*photon2) ;
511 (dynamic_cast<TParticle*>(particleList->At(index)))->SetStatusCode(0) ;
514 new((*plNe)[indexNe]) TParticle(*photon2) ;
515 (dynamic_cast<TParticle*>(plNe->At(indexNe))) ->SetStatusCode(0) ;
519 if(strstr(fOptionGJ,"deb all") || strstr(fOptionGJ,"deb"))
520 dynamic_cast<TH1F*>(fListHistos->FindObject("PtSpectra"))->
522 new((*plNePHOS)[indexNePHOS]) TParticle(*photon2) ;
523 (dynamic_cast<TParticle*>(plNePHOS->At(indexNePHOS))) ->SetStatusCode(0) ;
527 // e = (pGamma1+pGamma2).E();
528 // if(IsAngleInWindow(angle,e))
530 (fListHistos->FindObject("AnglePairAccepted"))->
531 Fill(pPi0.E(),angle);
534 }//photon2 in acceptance
535 }//if angle > mindist
538 }//for (iParticle<nParticle)
541 //Info("AddHIJINGToList","End HIJING");
544 //____________________________________________________________________________
545 Double_t AliPHOSGammaJet::CalculateJetRatioLimit(const Double_t ptg,
549 //Info("CalculateLimit","x1 %f, x2%f",x[0],x[1]);
550 Double_t epp = par[0] + par[1] * ptg ;
551 Double_t spp = par[2] + par[3] * ptg ;
552 Double_t f = x[0] + x[1] * ptg ;
553 Double_t epb = epp + par[4] ;
554 Double_t spb = TMath::Sqrt(spp*spp+ par[5]*par[5]) ;
555 Double_t rat = (epb - spb * f) / ptg ;
556 //Info("CalculateLimit","epp %f, spp %f, f %f", epp, spp, f);
557 //Info("CalculateLimit","epb %f, spb %f, rat %f", epb, spb, rat);
561 //____________________________________________________________________________
562 void AliPHOSGammaJet::CreateParticleList(Int_t iEvent,
563 TClonesArray * particleList,
566 TClonesArray * plNePHOS)
568 //Info("CreateParticleList","Inside");
569 AliPHOSGetter * gime = AliPHOSGetter::Instance(fInputFileName) ;
570 const AliPHOSGeometry * geom = gime->PHOSGeometry() ;
571 gime->Event(iEvent, "X") ;
574 Int_t index = particleList->GetEntries() ;
575 Int_t indexCh = plCh->GetEntries() ;
576 Int_t indexNe = plNe->GetEntries() ;
577 Int_t indexNePHOS = plNePHOS->GetEntries() ;
578 Int_t iParticle = 0 ;
579 Double_t charge = 0.;
581 Double_t x = 0., z = 0.;
584 for (iParticle=0 ; iParticle < gime->NPrimaries() ; iParticle++) {
585 const TParticle * particle = gime->Primary(iParticle) ;
591 //Keep Stable particles within eta range
592 if((particle->GetStatusCode() == 1) &&
593 (particle->Pt() > 0)){
594 if(TMath::Abs(particle->Eta())<fEtaCut){
597 charge = TDatabasePDG::Instance()
598 ->GetParticle(particle->GetPdgCode())->Charge();
599 if((charge != 0) && (particle->Pt() > fChargedPtCut)){
601 if(strstr(fOptionGJ,"deb all")|| strstr(fOptionGJ,"deb"))
602 dynamic_cast<TH1F*>(fListHistos->FindObject("PtSpectra"))->
603 Fill(particle->Pt());
604 new((*plCh)[indexCh++]) TParticle(*particle) ;
605 new((*particleList)[index++]) TParticle(*particle) ;
607 else if((charge == 0) && (particle->Pt() > fNeutralPtCut)){
608 geom->ImpactOnEmc(particle->Theta(),particle->Phi(), m,z,x);
611 if(strstr(fOptionGJ,"deb all")|| strstr(fOptionGJ,"deb"))
612 dynamic_cast<TH1F*>(fListHistos->FindObject("PtSpectra"))->
613 Fill(particle->Pt());
615 new((*plNePHOS)[indexNePHOS++]) TParticle(*particle) ;
617 if((particle->Phi()>fPhiEMCALCut[0]) &&
618 (particle->Phi()<fPhiEMCALCut[1]) && m == 0)
620 if(strstr(fOptionGJ,"deb all")|| strstr(fOptionGJ,"deb"))
621 dynamic_cast<TH1F*>(fListHistos->FindObject("PtSpectra"))->
622 Fill(particle->Pt());
623 new((*plNe)[indexNe++]) TParticle(*particle) ;
624 new((*particleList)[index++]) TParticle(*particle) ;
628 }//final particle, etacut
629 }//for (iParticle<nParticle)
633 Double_t mass = TDatabasePDG::Instance()->GetParticle(111)->Mass();
634 TLorentzVector pPi0, pGamma1, pGamma2 ;
635 Double_t angle = 0, cellDistance = 0.;
637 fFastRec = new AliPHOSFastGlobalReconstruction(fInputFileName);
638 fFastRec->FastReconstruction(iEvent);
640 for (iParticle=0 ; iParticle < gime->NPrimaries() ; iParticle++) {
641 const TParticle * particle = gime->Primary(iParticle) ;
645 //Keep Stable particles within eta range
646 if((particle->GetStatusCode() == 1) && (particle->Pt() > 0)){
650 charge = TDatabasePDG::Instance()
651 ->GetParticle(particle->GetPdgCode())->Charge();
652 if((charge != 0) && (particle->Pt() > fChargedPtCut) && (TMath::Abs(particle->Eta())<fEtaCut)){
653 new((*plCh)[indexCh++]) TParticle(*particle) ;
654 new((*particleList)[index++]) TParticle(*particle) ;
656 else if(charge == 0) {
657 geom->ImpactOnEmc(particle->Theta(),particle->Phi(), m,z,x);
658 if((particle->GetPdgCode() != 111) && particle->Pt() > 0 &&
659 (TMath::Abs(particle->Eta())<fEtaCut))
662 TLorentzVector part(particle->Px(),particle->Py(),
663 particle->Pz(),particle->Energy());
667 if(part.Pt() > fNeutralPtCut){
668 if(particle->Phi()>fPhiEMCALCut[0] &&
669 particle->Phi()<fPhiEMCALCut[1] && m == 0)
671 new((*particleList)[index]) TParticle(*particle) ;
672 (dynamic_cast<TParticle*>(particleList->At(index)))
673 ->SetMomentum(part.Px(),part.Py(),part.Pz(),part.E());
676 new((*plNe)[indexNe]) TParticle(*particle) ;
677 (dynamic_cast<TParticle*>(plNe->At(indexNe)))
678 ->SetMomentum(part.Px(),part.Py(),part.Pz(),part.E());
683 new((*plNePHOS)[indexNePHOS]) TParticle(*particle) ;
684 (dynamic_cast<TParticle*>(plNePHOS->At(indexNePHOS)))
685 ->SetMomentum(part.Px(),part.Py(),part.Pz(),part.E());
690 if((particle->GetPdgCode() == 111) && (particle->Pt() > fNeutralPtCut) &&
691 (TMath::Abs(particle->Eta())<fEtaCut+1))
694 pPi0.SetPxPyPzE(particle->Px(),particle->Py(),particle->Pz(),
699 Pi0Decay(mass,pPi0,pGamma1,pGamma2,angle);
700 //Check if decay photons are too close for PHOS
701 cellDistance = angle*460; //cm
703 if (cellDistance < fMinDistance) {
705 //Pi0 inside phi EMCAL acceptance
708 TLorentzVector part(particle->Px(),particle->Py(),
709 particle->Pz(),particle->Energy());
712 if(part.Pt() > fNeutralPtCut){
713 if(particle->Phi()>fPhiEMCALCut[0] &&
714 particle->Phi()<fPhiEMCALCut[1] && m == 0){
715 if(strstr(fOptionGJ,"deb all")|| strstr(fOptionGJ,"deb"))
716 dynamic_cast<TH1F*>(fListHistos->FindObject("PtSpectra"))->
717 Fill(particle->Pt());
719 new((*plNe)[indexNe]) TParticle(*particle) ;
720 (dynamic_cast<TParticle*>(plNe->At(indexNe)))
721 ->SetMomentum(part.Px(),part.Py(),part.Pz(),part.E());
722 new((*particleList)[index]) TParticle(*particle) ;
723 (dynamic_cast<TParticle*>(particleList->At(index)))
724 ->SetMomentum(part.Px(),part.Py(),part.Pz(),part.E());
729 if(strstr(fOptionGJ,"deb all")|| strstr(fOptionGJ,"deb"))
730 dynamic_cast<TH1F*>(fListHistos->FindObject("PtSpectra"))->
731 Fill(particle->Pt());
732 new((*plNePHOS)[indexNePHOS]) TParticle(*particle) ;
733 (dynamic_cast<TParticle*>(plNePHOS->At(indexNePHOS)))
734 ->SetMomentum(part.Px(),part.Py(),part.Pz(),part.E());
741 dynamic_cast<TH2F*>(fListHistos->FindObject("AnglePair"))
742 ->Fill(pPi0.E(),angle);
746 if(pGamma1.Pt() > 0 && TMath::Abs(pGamma1.Eta())<fEtaCut){
749 if(pGamma1.Pt() > fNeutralPtCut){
750 TParticle * photon1 =
751 new TParticle(22,1,0,0,0,0,pGamma1.Px(),pGamma1.Py(),
752 pGamma1.Pz(),pGamma1.E(),0,0,0,0);
753 geom->ImpactOnEmc(photon1->Theta(),photon1->Phi(), m,z,x);
754 if( photon1->Phi()>fPhiEMCALCut[0] && photon1->Phi()<fPhiEMCALCut[1]
756 if(strstr(fOptionGJ,"deb all") || strstr(fOptionGJ,"deb"))
757 dynamic_cast<TH1F*>(fListHistos->FindObject("PtSpectra"))->
759 new((*plNe)[indexNe++]) TParticle(*photon1) ;
760 new((*particleList)[index++]) TParticle(*photon1) ;
765 if(strstr(fOptionGJ,"deb all") || strstr(fOptionGJ,"deb"))
766 dynamic_cast<TH1F*>(fListHistos->FindObject("PtSpectra"))->
768 new((*plNePHOS)[indexNePHOS++]) TParticle(*photon1) ;
773 if(pGamma2.Pt() > 0 && TMath::Abs(pGamma2.Eta())<fEtaCut){
776 if(pGamma2.Pt() > fNeutralPtCut ){
777 TParticle * photon2 =
778 new TParticle(22,1,0,0,0,0,pGamma2.Px(), pGamma2.Py(),
779 pGamma2.Pz(),pGamma2.E(),0,0,0,0);
780 geom->ImpactOnEmc(photon2->Theta(),photon2->Phi(), m,z,x);
781 if(photon2->Phi()>fPhiEMCALCut[0] &&
782 photon2->Phi()<fPhiEMCALCut[1] && m == 0){
783 if(strstr(fOptionGJ,"deb all") || strstr(fOptionGJ,"deb"))
784 dynamic_cast<TH1F*>(fListHistos->FindObject("PtSpectra"))->
786 new((*plNe)[indexNe++]) TParticle(*photon2) ;
787 new((*particleList)[index++]) TParticle(*photon2) ;
790 if(strstr(fOptionGJ,"deb all") || strstr(fOptionGJ,"deb"))
791 dynamic_cast<TH1F*>(fListHistos->FindObject("PtSpectra"))->
793 new((*plNePHOS)[indexNePHOS++]) TParticle(*photon2) ;
797 // Float_t e = (pGamma1+pGamma2).E();
798 // if(IsAngleInWindow(angle,e))
800 (fListHistos->FindObject("AnglePairAccepted"))->
801 Fill(pPi0.E(),angle);
804 }//photon2 in acceptance
805 }//if angle > mindist
808 }//final particle, etacut
809 }//for (iParticle<nParticle)
813 //____________________________________________________________________________
814 void AliPHOSGammaJet::CreateParticleListFromESD(TClonesArray * pl,
817 TClonesArray * plNePHOS,
820 //Create a list of particles from the ESD. These particles have been measured
821 //by the Central Tracking system (TPC), PHOS and EMCAL
822 //(EMCAL not available for the moment).
823 //Info("CreateParticleListFromESD","Inside");
825 Int_t index = pl->GetEntries() ;
827 Float_t *pid = new Float_t[AliPID::kSPECIESN];
829 //########### PHOS ##############
830 //Info("CreateParticleListFromESD","Fill ESD PHOS list");
831 Int_t begphos = esd->GetFirstPHOSCluster();
832 Int_t endphos = esd->GetFirstPHOSCluster() +
833 esd->GetNumberOfPHOSClusters() ;
834 Int_t indexNePHOS = plNePHOS->GetEntries() ;
835 if(strstr(fOptionGJ,"deb all"))
836 Info("CreateParticleListFromESD","PHOS: first particle %d, last particle %d",
839 for (npar = begphos; npar < endphos; npar++) {//////////////PHOS track loop
840 AliESDCaloCluster * clus = esd->GetCaloCluster(npar) ; // retrieve track from esd
842 //Create a TParticle to fill the particle list
844 Float_t en = clus->GetClusterEnergy() ;
845 Float_t *p = new Float_t();
846 clus->GetGlobalPosition(p) ;
847 TVector3 pos(p[0],p[1],p[2]) ;
848 Double_t phi = pos.Phi();
849 Double_t theta= pos.Theta();
850 Double_t px = en*TMath::Cos(phi)*TMath::Sin(theta);;
851 Double_t py = en*TMath::Sin(phi)*TMath::Sin(theta);
852 Double_t pz = en*TMath::Cos(theta);
854 TParticle * particle = new TParticle() ;
855 particle->SetMomentum(px,py,pz,en) ;
857 //Select only photons
860 //cout<<"pid "<<pid[AliPID::kPhoton]<<endl ;
861 if( pid[AliPID::kPhoton] > 0.75)
862 new((*plNePHOS)[indexNePHOS++]) TParticle(*particle) ;
865 //########### TPC #####################
866 //Info("CreateParticleListFromESD","Fill ESD TPC list");
868 Int_t endtpc = esd->GetNumberOfTracks() ;
869 Int_t indexCh = plCh->GetEntries() ;
870 if(strstr(fOptionGJ,"deb all"))
871 Info("CreateParticleListFromESD","TPC: first particle %d, last particle %d",
874 for (npar = begtpc; npar < endtpc; npar++) {////////////// track loop
875 AliESDtrack * track = esd->GetTrack(npar) ; // retrieve track from esd
877 Double_t en = track ->GetTPCsignal() ;
879 track->GetPxPyPz(mom) ;
880 Double_t px = mom[0];
881 Double_t py = mom[1];
882 Double_t pz = mom[2]; //Check with TPC people if this is correct.
884 //cout<<"TPC signal "<<en<<endl;
885 //cout<<"px "<<px<<"; py "<<py<<"; pz "<<pz<<endl;
886 TParticle * particle = new TParticle() ;
887 particle->SetMomentum(px,py,pz,en) ;
889 new((*plCh)[indexCh++]) TParticle(*particle) ;
890 new((*pl)[index++]) TParticle(*particle) ;
894 //################ EMCAL ##############
895 Double_t v[3] ; //vertex ;
896 esd->GetVertex()->GetXYZ(v) ;
897 //##########Uncomment when ESD for EMCAL works ##########
898 //Info("CreateParticleListFromESD","Fill ESD EMCAL list");
900 Int_t begem = esd->GetFirstEMCALCluster();
901 Int_t endem = esd->GetFirstEMCALCluster() +
902 esd->GetNumberOfEMCALClusters() ;
903 Int_t indexNe = plNe->GetEntries() ;
904 if(strstr(fOptionGJ,"deb all"))
905 Info("CreateParticleListFromESD","EMCAL: first particle %d, last particle %d",
908 for (npar = begem; npar < endem; npar++) {//////////////EMCAL track loop
909 AliESDCaloCluster * clus = esd->GetCaloCluster(npar) ; // retrieve track from esd
911 Float_t en = clus->GetClusterEnergy() ;
912 Float_t *p = new Float_t();
913 clus->GetGlobalPosition(p) ;
914 TVector3 pos(p[0],p[1],p[2]) ;
915 Double_t phi = pos.Phi();
916 Double_t theta= pos.Theta();
917 Double_t px = en*TMath::Cos(phi)*TMath::Sin(theta);;
918 Double_t py = en*TMath::Sin(phi)*TMath::Sin(theta);
919 Double_t pz = en*TMath::Cos(theta);
920 //cout<<"EMCAL signal "<<en<<endl;
921 //cout<<"px "<<px<<"; py "<<py<<"; pz "<<pz<<endl;
922 //TParticle * particle = new TParticle() ;
923 //particle->SetMomentum(px,py,pz,en) ;
926 // //Uncomment if PID IS WORKING, photon and pi0 idenitification.
927 // //if( pid[AliPID::kPhoton] > 0.75) //This has to be fixen.
929 // //else if( pid[AliPID::kPi0] > 0.75)
931 pdg = 22; //No PID, assume all photons
932 TParticle * particle = new TParticle(pdg, 1, -1, -1, -1, -1,
933 px, py, pz, en, v[0], v[1], v[2], 0);
935 new((*plNe)[indexNe++]) TParticle(*particle) ;
936 new((*pl)[index++]) TParticle(*particle) ;
939 // Info("CreateParticleListFromESD","End Inside");
945 //____________________________________________________________________________
946 void AliPHOSGammaJet::Exec(Option_t *option)
956 // Create chain of esd trees
957 const UInt_t kNevent = static_cast<UInt_t>(GetNEvent()) ;
958 t = ReadESD(kNevent, fDirName, fESDTree, fPattern) ;
960 AliError("Could not create the TChain") ;
965 t->SetBranchAddress("ESD",&esd); // point to the container esd where to put the event from the esdTree
970 // AliGenPythia* pyth = (AliGenPythia*) gAlice->Generator();
973 TClonesArray * particleList = new TClonesArray("TParticle",1000);
974 TClonesArray * plCh = new TClonesArray("TParticle",1000);
975 TClonesArray * plNe = new TClonesArray("TParticle",1000);
976 TClonesArray * plNePHOS = new TClonesArray("TParticle",1000);
978 for (Int_t iEvent = 0 ; iEvent < fNEvent ; iEvent++) {
979 if(strstr(fOptionGJ,"deb")||strstr(fOptionGJ,"deb all"))
980 Info("Exec", "Event %d", iEvent) ;
984 Double_t phig = 0., phil = 0., phich = 0 , phipi = 0;
985 Double_t etag = 0., etal = 0., etach = 0., etapi = 0. ;
986 Double_t ptg = 0., ptl = 0., ptch = 0., ptpi = 0.;
988 TLorentzVector jet (0,0,0,0);
989 TLorentzVector jettpc(0,0,0,0);
993 Int_t iNbytes = t->GetEntry(iEvent); // store event in esd
994 //cout<<"nbytes "<<iNbytes<<endl;
995 if ( iNbytes == 0 ) {
996 AliError("Empty TChain") ;
999 CreateParticleListFromESD(particleList, plCh,plNe,plNePHOS, esd); //,iEvent);
1002 CreateParticleList(iEvent, particleList, plCh,plNe,plNePHOS);
1004 // TLorentzVector pyjet(0,0,0,0);
1007 // Float_t jets[4][10];
1008 // pyth->SetJetReconstructionMode(1);
1009 // pyth->LoadEvent();
1010 // pyth->GetJets(nJ, nJT, jets);
1012 // Float_t pxJ = jets[0][0];
1013 // Float_t pyJ = jets[1][0];
1014 // Float_t pzJ = jets[2][0];
1015 // Float_t eJ = jets[3][0];
1016 // pyjet.SetPxPyPzE(pxJ,pyJ,pzJ,eJ ) ;
1019 // //Info("Exec",">>>>>>>>>>Number of jets !!!! %d",nJT);
1020 // for (Int_t iJ = 1; iJ < nJT; iJ++) {
1021 // Float_t pxJ = jets[0][iJ];
1022 // Float_t pyJ = jets[1][iJ];
1023 // Float_t pzJ = jets[2][iJ];
1024 // Float_t eJ = jets[3][iJ];
1025 // pyjet.SetPxPyPzE(pxJ,pyJ,pzJ,eJ ) ;
1026 // //Info("Exec",">>>>>Pythia Jet: %d, Phi %f, Eta %f, Pt %f",
1027 // // iJ,pyjet.Phi(),pyjet.Eta(),pyjet.Pt());
1033 AddHIJINGToList(iEvent, particleList, plCh,plNe, plNePHOS);
1036 Bool_t iIsInPHOS = kFALSE ;
1037 GetGammaJet(plNePHOS, ptg, phig, etag, iIsInPHOS) ;
1041 //Info("Exec"," In PHOS") ;
1042 dynamic_cast<TH1F*>(fListHistos->FindObject("NGamma"))->Fill(ptg);
1043 dynamic_cast<TH2F*>(fListHistos->FindObject("PhiGamma"))
1045 dynamic_cast<TH2F*>(fListHistos->FindObject("EtaGamma"))
1047 if(strstr(fOptionGJ,"deb")||strstr(fOptionGJ,"deb all"))
1048 Info("Exec", "Gamma: pt %f, phi %f, eta %f", ptg,
1051 // cout<<"n charged "<<plCh->GetEntries()<<endl;
1052 // cout<<"n neutral "<<plNe->GetEntries()<<endl;
1053 // cout<<"n All "<<particleList->GetEntries()<<endl;
1055 GetLeadingCharge(plCh, ptg, phig, ptch, etach, phich) ;
1056 GetLeadingPi0 (plNe, ptg, phig, ptpi, etapi, phipi) ;
1058 // cout<<"n2 charged "<<plCh->GetEntries()<<endl;
1059 // cout<<"n2 neutral "<<plNe->GetEntries()<<endl;
1060 // cout<<"n2 All "<<particleList->GetEntries()<<endl;
1065 //Is the leading cone inside EMCAL?
1066 Bool_t insidech = kFALSE ;
1067 if((phich - fCone) > fPhiEMCALCut[0] &&
1068 (phich + fCone) < fPhiEMCALCut[1]){
1071 Bool_t insidepi = kFALSE ;
1072 if((phipi - fCone) > fPhiEMCALCut[0] &&
1073 (phipi + fCone) < fPhiEMCALCut[1]){
1077 if ((ptch > 0 || ptpi > 0)){
1078 if((ptch > ptpi) && insidech){
1082 dynamic_cast<TH2F*>(fListHistos->FindObject("ChargeRatio"))
1083 ->Fill(ptg,ptch/ptg);
1084 dynamic_cast<TH2F*>(fListHistos->FindObject("DeltaPhiCharge"))
1085 ->Fill(ptg,phig-phich);
1086 dynamic_cast<TH2F*>(fListHistos->FindObject("DeltaEtaCharge"))
1087 ->Fill(ptg,etag-etach);
1088 if(strstr(fOptionGJ,"deb"))
1089 Info("Exec"," Charged Leading") ;
1091 if((ptpi > ptch) && insidepi){
1096 dynamic_cast<TH2F*>(fListHistos->FindObject("Pi0Ratio"))
1097 ->Fill(ptg,ptpi/ptg);
1098 dynamic_cast<TH2F*>(fListHistos->FindObject("DeltaPhiPi0"))
1099 ->Fill(ptg,phig-phipi);
1100 dynamic_cast<TH2F*>(fListHistos->FindObject("DeltaEtaPi0"))
1101 ->Fill(ptg,etag-etapi);
1103 if(ptpi > 0. && strstr(fOptionGJ,"deb"))
1104 Info("Exec"," Pi0 Leading") ;
1107 if(strstr(fOptionGJ,"deb"))
1108 Info("Exec","Leading pt %f, phi %f",ptl,phil);
1109 if(insidech || insidepi){
1112 MakeJet(particleList, ptg, phig, ptl, phil, etal, "", jet);
1114 if(strstr(fOptionGJ,"deb")){
1115 // Info("Exec","Pythia Jet: Phi %f, Eta %f, Pt %f",
1116 // pyjet.Phi(),pyjet.Eta(),pyjet.Pt());
1117 Info("Exec","TPC+EMCAL Jet: Phi %f, Eta %f, Pt %f",
1118 jet.Phi(),jet.Eta(),jet.Pt());
1120 // dynamic_cast<TH2F*>(fListHistos->FindObject("DeltaPhiJet"))
1121 // ->Fill(ptg,pyjet.Phi()-jet.Phi());
1122 // dynamic_cast<TH2F*>(fListHistos->FindObject("DeltaEtaJet"))
1123 // ->Fill(ptg,pyjet.Eta()-jet.Eta());
1124 // dynamic_cast<TH2F*>(fListHistos->FindObject("DeltaPtJet"))
1125 // ->Fill(ptg,pyjet.Pt()-jet.Pt());
1128 MakeJetAnyConeOrPt(particleList, ptg, phig, ptl, phil, etal, "");
1132 if(fOnlyCharged && ptch > 0.)
1134 if(strstr(fOptionGJ,"deb"))
1135 Info("Exec","Leading TPC pt %f, phi %f",ptch,phich);
1137 dynamic_cast<TH2F*>(fListHistos->FindObject("TPCRatio"))
1138 ->Fill(ptg,ptch/ptg);
1139 dynamic_cast<TH2F*>(fListHistos->FindObject("DeltaPhiTPC"))
1140 ->Fill(ptg,phig-phich);
1141 dynamic_cast<TH2F*>(fListHistos->FindObject("DeltaEtaTPC"))
1142 ->Fill(ptg,etag-etach);
1146 MakeJet(plCh, ptg, phig, ptch, phich, etach, "TPC",jettpc);
1148 if(strstr(fOptionGJ,"deb")){
1149 // Info("Exec","Pythia Jet: Phi %f, Eta %f, Pt %f",
1150 // pyjet.Phi(),pyjet.Eta(),pyjet.Pt());
1151 Info("Exec","TPC Jet: Phi %f, Eta %f, Pt %f",
1152 jettpc.Phi(),jettpc.Eta(),jettpc.Pt());
1154 // dynamic_cast<TH2F*>(fListHistos->FindObject("DeltaPhiTPCJet"))
1155 // ->Fill(ptg,pyjet.Phi()-jettpc.Phi());
1156 // dynamic_cast<TH2F*>(fListHistos->FindObject("DeltaEtaTPCJet"))
1157 // ->Fill(ptg,pyjet.Eta()-jettpc.Eta());
1158 // dynamic_cast<TH2F*>(fListHistos->FindObject("DeltaPtTPCJet"))
1159 // ->Fill(ptg,pyjet.Pt()-jettpc.Pt());
1162 MakeJetAnyConeOrPt(plCh, ptg, phig, ptch, phich, etach, "TPC");
1168 particleList->Delete() ;
1171 plNePHOS->Delete() ;
1176 delete particleList ;
1178 fOutputFile->Write() ;
1183 //____________________________________________________________________________
1184 void AliPHOSGammaJet::FillJetHistos(TClonesArray * pl, Double_t ptg,
1185 TString conf, TString type)
1187 //Fill jet fragmentation histograms if !fAnyCone,
1188 //only for fCone and fPtThres
1189 TParticle * particle = 0 ;
1194 while ( (particle = dynamic_cast<TParticle*>(next())) ) {
1196 Double_t pt = particle->Pt();
1198 charge = TDatabasePDG::Instance()
1199 ->GetParticle(particle->GetPdgCode())->Charge();
1200 if(charge != 0){//Only jet Charged particles
1202 (fListHistos->FindObject(type+conf+"Fragment"))
1205 (fListHistos->FindObject(type+conf+"PtDist"))
1211 (fListHistos->FindObject("NBkg"+conf))->Fill(ipr);
1214 //____________________________________________________________________________
1215 void AliPHOSGammaJet::FillJetHistosAnyConeOrPt(TClonesArray * pl, Double_t ptg,
1216 TString conf, TString type,
1217 TString cone, TString ptcut)
1219 //Fill jet fragmentation histograms if fAnyCone,
1220 //for several cones and pt thresholds
1221 TParticle *particle = 0;
1226 while ( (particle = dynamic_cast<TParticle*>(next())) ) {
1228 Double_t pt = particle->Pt();
1229 charge = TDatabasePDG::Instance()
1230 ->GetParticle(particle->GetPdgCode())->Charge();
1231 if(charge != 0){//Only jet Charged particles
1233 (fListHistos->FindObject(type+conf+"FragmentCone"+cone+"Pt"+ptcut))
1236 (fListHistos->FindObject(type+conf+"PtDistCone"+cone+"Pt"+ptcut))
1243 (fListHistos->FindObject("NBkg"+conf+"Cone"+cone+"Pt"+ptcut))
1248 //____________________________________________________________________________
1249 void AliPHOSGammaJet::GetGammaJet(TClonesArray * pl, Double_t &pt,
1250 Double_t &phi, Double_t &eta, Bool_t &Is) const
1252 //Search for the prompt photon in PHOS with pt > fPtCut
1257 for(Int_t ipr = 0;ipr < pl->GetEntries() ; ipr ++ ){
1258 TParticle * particle = dynamic_cast<TParticle *>(pl->At(ipr)) ;
1260 if((particle->Pt() > fPtCut) && (particle->Pt() > pt)){
1262 pt = particle->Pt();
1263 phi = particle->Phi() ;
1264 eta = particle->Eta() ;
1270 //____________________________________________________________________________
1271 void AliPHOSGammaJet::GetLeadingCharge(TClonesArray * pl,
1272 Double_t ptg, Double_t phig,
1273 Double_t &pt, Double_t &eta, Double_t &phi) const
1275 //Search for the charged particle with highest with
1276 //Phi=Phi_gamma-Pi and pT=0.1E_gamma
1281 for(Int_t ipr = 0;ipr < pl->GetEntries() ; ipr ++ ){
1283 TParticle * particle = dynamic_cast<TParticle *>(pl->At(ipr)) ;
1285 Double_t ptl = particle->Pt();
1286 Double_t rat = ptl/ptg ;
1287 Double_t phil = particle->Phi() ;
1289 if(((phig-phil)> fPhiMinCut) && ((phig-phil)<fPhiMaxCut) &&
1290 (rat > fRatioMinCut) && (rat < fRatioMaxCut) && (ptl > pt)) {
1291 eta = particle->Eta() ;
1294 //printf("GetLeadingCharge: %f %f %f %f \n", pt, eta, phi,rat) ;
1297 //printf("GetLeadingCharge: %f %f %f \n", pt, eta, phi) ;
1302 //____________________________________________________________________________
1303 void AliPHOSGammaJet::GetLeadingPi0(TClonesArray * pl,
1304 Double_t ptg, Double_t phig,
1305 Double_t &pt, Double_t &eta, Double_t &phi)
1308 //Search for the neutral pion with highest with
1309 //Phi=Phi_gamma-Pi and pT=0.1E_gamma
1313 Double_t ptl = -100.;
1314 Double_t rat = -100.;
1315 Double_t phil = -100. ;
1318 TParticle * particle = 0;
1322 while ( (particle = (TParticle*)next()) ) {
1323 if( particle->GetPdgCode() == 111){
1324 ptl = particle->Pt();
1326 phil = particle->Phi() ;
1327 e = particle->Energy();
1329 (fListHistos->FindObject("AnglePairNoCut"))->
1332 (fListHistos->FindObject("InvMassPairNoCut"))->
1335 if(((phig-phil)> fPhiMinCut) && ((phig-phil)<fPhiMaxCut) &&
1336 (rat > fRatioMinCut) && (rat < fRatioMaxCut)) {
1339 (fListHistos->FindObject("AnglePairLeadingCut"))->
1342 (fListHistos->FindObject("InvMassPairLeadingCut"))->
1346 (fListHistos->FindObject("AnglePairAngleCut"))->
1349 (fListHistos->FindObject("InvMassPairAngleCut"))->
1353 (fListHistos->FindObject("InvMassPairAllCut"))->
1356 (fListHistos->FindObject("AnglePairAllCut"))->
1361 eta = particle->Eta() ;
1365 //printf("GetLeadingPi0: %f %f %f %f %f \n", pt, eta, phi, rat, ptg) ;
1371 (fListHistos->FindObject("InvMassPairLeading"))->
1374 (fListHistos->FindObject("AnglePairLeading"))->
1379 Int_t iPrimary = -1;
1380 TLorentzVector gammai,gammaj;
1381 Double_t angle = 0., e = 0., invmass = 0.;
1382 Double_t anglef = 0., ef = 0., invmassf = 0.;
1386 while ( (particle = (TParticle*)next()) ) {
1389 ksPdg = particle->GetPdgCode();
1390 ptl = particle->Pt();
1391 if(ksPdg == 111){ //2 gamma
1393 phil = particle->Phi() ;
1394 if((ptl> pt)&& (rat > fRatioMinCut) && (rat < fRatioMaxCut) &&
1395 ((phig-phil)>fPhiMinCut)&&((phig-phil)<fPhiMaxCut)){
1396 eta = particle->Eta() ;
1401 if(ksPdg == 22){//1 gamma
1402 particle->Momentum(gammai);
1405 while ( (particle = (TParticle*)next2()) ) {
1407 if(jPrimary>iPrimary){
1408 ksPdg = particle->GetPdgCode();
1410 particle->Momentum(gammaj);
1411 //Info("GetLeadingPi0","Egammai %f, Egammaj %f",
1412 //gammai.Pt(),gammaj.Pt());
1414 ptl = (gammai+gammaj).Pt();
1415 phil = (gammai+gammaj).Phi();
1417 phil+=TMath::TwoPi();
1419 invmass = (gammai+gammaj).M();
1420 angle = gammaj.Angle(gammai.Vect());
1421 //Info("GetLeadingPi0","Angle %f", angle);
1422 e = (gammai+gammaj).E();
1425 (fListHistos->FindObject("AnglePairNoCut"))->
1428 (fListHistos->FindObject("InvMassPairNoCut"))->
1431 if((rat > fRatioMinCut) && (rat < fRatioMaxCut) &&
1432 ((phig-phil)>fPhiMinCut)&&((phig-phil)<fPhiMaxCut)){
1435 (fListHistos->FindObject("AnglePairLeadingCut"))->
1438 (fListHistos->FindObject("InvMassPairLeadingCut"))->
1441 if(IsAngleInWindow(angle,e)){
1443 (fListHistos->FindObject("AnglePairAngleCut"))->
1446 (fListHistos->FindObject("InvMassPairAngleCut"))->
1449 //Info("GetLeadingPi0","InvMass %f", invmass);
1450 if((invmass>fInvMassMinCut) && (invmass<fInvMassMaxCut)){
1452 (fListHistos->FindObject("InvMassPairAllCut"))->
1455 (fListHistos->FindObject("AnglePairAllCut"))->
1459 eta = particle->Eta() ;
1463 invmassf = invmass ;
1467 }//(invmass>0.125) && (invmass<0.145)
1468 }//gammaj.Angle(gammai.Vect())<0.04
1470 }//iprimary<jprimary
1478 (fListHistos->FindObject("InvMassPairLeading"))->
1481 (fListHistos->FindObject("AnglePairLeading"))->
1485 // printf("GetLeadingPi0: %f %f %f \n", pt, eta, phi) ;
1489 //____________________________________________________________________________
1490 void AliPHOSGammaJet::InitParameters()
1493 //Initialize the parameters of the analysis.
1495 fAngleMaxParam.Set(4) ;
1496 fAngleMaxParam.AddAt(0.4,0);//={0.4,-0.25,0.025,-2e-4};
1497 fAngleMaxParam.AddAt(-0.25,1) ;
1498 fAngleMaxParam.AddAt(0.025,2) ;
1499 fAngleMaxParam.AddAt(-2e-4,3) ;
1500 fAnyConeOrPt = kFALSE ;
1501 fOutputFileName = "GammaJet.root" ;
1503 fHIJINGFileName = "galice.root" ;
1505 fMinDistance = 3.6 ;
1507 fInvMassMaxCut = 0.15 ;
1508 fInvMassMinCut = 0.12 ;
1509 fOnlyCharged = kFALSE ;
1512 fPhiEMCALCut[0] = 60. *TMath::Pi()/180.;
1513 fPhiEMCALCut[1] = 180.*TMath::Pi()/180.;
1517 fNeutralPtCut = 0.5 ;
1518 fChargedPtCut = 0.5 ;
1519 fTPCCutsLikeEMCAL = kFALSE ;
1520 //Jet selection parameters
1522 fRatioMaxCut = 1.0 ;
1523 fRatioMinCut = 0.1 ;
1524 fJetRatioMaxCut = 1.2 ;
1525 fJetRatioMinCut = 0.8 ;
1526 fJetTPCRatioMaxCut = 1.2 ;
1527 fJetTPCRatioMinCut = 0.3 ;
1531 fESDTree = "esdTree" ;
1534 //Cut depending on gamma energy
1536 fPtJetSelectionCut = 20.; //For Low pt jets+BKG, another limits applyed
1537 //Reconstructed jet energy dependence parameters
1538 //e_jet = a1+e_gamma b2.
1539 //Index 0-> Pt>2 GeV r = 0.3; Index 1-> Pt>0.5 GeV r = 0.3
1540 fJetE1[0] = -5.75; fJetE1[1] = -4.1;
1541 fJetE2[0] = 1.005; fJetE2[1] = 1.05;
1543 //Reconstructed sigma of jet energy dependence parameters
1544 //s_jet = a1+e_gamma b2.
1545 //Index 0-> Pt>2 GeV r = 0.3; Index 1-> Pt>0.5 GeV r = 0.3
1546 fJetSigma1[0] = 2.65; fJetSigma1[1] = 2.75;
1547 fJetSigma2[0] = 0.0018; fJetSigma2[1] = 0.033;
1549 //Background mean energy and RMS
1550 //Index 0-> No BKG; Index 1-> BKG > 2 GeV;
1551 //Index 2-> (low pt jets)BKG > 0.5 GeV;
1552 //Index > 2, same for TPC conf
1553 fBkgMean[0] = 0.; fBkgMean[1] = 8.8 ; fBkgMean[2] = 69.5;
1554 fBkgMean[3] = 0.; fBkgMean[4] = 6.4; fBkgMean[5] = 48.6;
1555 fBkgRMS[0] = 0.; fBkgRMS[1] = 7.5; fBkgRMS[2] = 22.0;
1556 fBkgRMS[3] = 0.; fBkgRMS[4] = 5.4; fBkgRMS[5] = 13.2;
1558 //Factor x of min/max = E -+ x * sigma. Obtained after selecting the
1559 //limits for monoenergetic jets.
1560 //Index 0-> No BKG; Index 1-> BKG > 2 GeV;
1561 //Index 2-> (low pt jets) BKG > 0.5 GeV;
1562 //Index > 2, same for TPC conf
1564 fJetXMin1[0] =-0.69 ; fJetXMin1[1] = 0.39 ; fJetXMin1[2] =-0.88 ;
1565 fJetXMin1[3] =-2.0 ; fJetXMin1[4] =-0.442 ; fJetXMin1[5] =-1.1 ;
1566 fJetXMin2[0] = 0.066; fJetXMin2[1] = 0.038; fJetXMin2[2] = 0.034;
1567 fJetXMin2[3] = 0.25 ; fJetXMin2[4] = 0.113; fJetXMin2[5] = 0.077 ;
1568 fJetXMax1[0] =-3.8 ; fJetXMax1[1] =-0.76 ; fJetXMax1[2] =-3.6 ;
1569 fJetXMax1[3] =-2.7 ; fJetXMax1[4] =-1.21 ; fJetXMax1[5] =-3.7 ;
1570 fJetXMax2[0] =-0.012; fJetXMax2[1] =-0.022; fJetXMax2[2] = 0.016;
1571 fJetXMax2[3] =-0.024; fJetXMax2[4] =-0.008; fJetXMax2[5] = 0.027;
1574 //Photon fast reconstruction
1575 fResPara1 = 0.0255 ; // GeV
1576 fResPara2 = 0.0272 ;
1577 fResPara3 = 0.0129 ;
1579 fPosParaA = 0.096 ; // cm
1582 //Different cones and pt thresholds to construct the jet
1588 fCones[0] = 0.3 ; fNameCones[0] = "03" ;
1589 fPtThres[0] = 0.5 ; fNamePtThres[0] = "05" ;
1593 //__________________________________________________________________________-
1594 Bool_t AliPHOSGammaJet::IsAngleInWindow(const Float_t angle,const Float_t e) {
1595 //Check if the opening angle of the candidate pairs is inside
1596 //our selection windowd
1597 Bool_t result = kFALSE;
1598 Double_t mpi0 = 0.1349766;
1599 Double_t max = fAngleMaxParam.At(0)*TMath::Exp(fAngleMaxParam.At(1)*e)
1600 +fAngleMaxParam.At(2)+fAngleMaxParam.At(3)*e;
1601 Double_t arg = (e*e-2*mpi0*mpi0)/(e*e);
1602 Double_t min = 100. ;
1604 min = TMath::ACos(arg);
1606 if((angle<max)&&(angle>=min))
1612 //__________________________________________________________________________-
1613 Bool_t AliPHOSGammaJet::IsJetSelected(const Double_t ptg, const Double_t ptj,
1614 const TString type ){
1615 //Check if the energy of the reconstructed jet is within an energy window
1623 if(type == "TPC" && !fTPCCutsLikeEMCAL){
1624 iTPC = 3 ;//If(fTPCCutsLikeEMCAL) take jet energy cuts like EMCAL
1629 //Phythia alone, jets with pt_th > 0.2, r = 0.3
1630 par[0] = fJetE1[0]; par[1] = fJetE2[0];
1631 //Energy of the jet peak
1632 //e_jet = fJetE1[0]+fJetE2[0]*e_gamma, simulation fit
1633 par[2] = fJetSigma1[0]; par[3] = fJetSigma2[0];
1634 //Sigma of the jet peak
1635 //sigma_jet = fJetSigma1[0]+fJetSigma2[0]*e_gamma, simulation fit
1636 par[4] = fBkgMean[0 + iTPC]; par[5] = fBkgRMS[0 + iTPC];
1637 //Parameters reserved for HIJING bkg.
1638 xmax[0] = fJetXMax1[0 + iTPC]; xmax[1] = fJetXMax2[0 + iTPC];
1639 xmin[0] = fJetXMin1[0 + iTPC]; xmin[1] = fJetXMin2[0 + iTPC];
1640 //Factor that multiplies sigma to obtain the best limits,
1641 //by observation, of mono jet ratios (ptjet/ptg)
1642 //X_jet = fJetX1[0]+fJetX2[0]*e_gamma
1646 if(ptg > fPtJetSelectionCut){
1647 //Phythia +HIJING with pt_th > 2 GeV/c, r = 0.3
1648 par[0] = fJetE1[0]; par[1] = fJetE2[0];
1649 //Energy of the jet peak, same as in pp
1650 //e_jet = fJetE1[0]+fJetE2[0]*e_gamma, simulation fit
1651 par[2] = fJetSigma1[0]; par[3] = fJetSigma2[0];
1652 //Sigma of the jet peak, same as in pp
1653 //sigma_jet = fJetSigma1[0]+fJetSigma2[0]*e_gamma, simulation fit
1654 par[4] = fBkgMean[1 + iTPC]; par[5] = fBkgRMS[1 + iTPC];
1655 //Mean value and RMS of HIJING Bkg
1656 xmax[0] = fJetXMax1[1 + iTPC]; xmax[1] = fJetXMax2[1 + iTPC];
1657 xmin[0] = fJetXMin1[1 + iTPC]; xmin[1] = fJetXMin2[1 + iTPC];
1658 //Factor that multiplies sigma to obtain the best limits,
1659 //by observation, of mono jet ratios (ptjet/ptg) mixed with HIJING Bkg,
1660 //pt_th > 2 GeV, r = 0.3
1661 //X_jet = fJetX1[0]+fJetX2[0]*e_gamma
1665 //Phythia + HIJING with pt_th > 0.5 GeV/c, r = 0.3
1666 par[0] = fJetE1[1]; par[1] = fJetE2[1];
1667 //Energy of the jet peak, pt_th > 2 GeV/c, r = 0.3
1668 //e_jet = fJetE1[0]+fJetE2[0]*e_gamma, simulation fit
1669 par[2] = fJetSigma1[1]; par[3] = fJetSigma2[1];
1670 //Sigma of the jet peak, pt_th > 2 GeV/c, r = 0.3
1671 //sigma_jet = fJetSigma1[0]+fJetSigma2[0]*e_gamma, simulation fit
1672 par[4] = fBkgMean[2 + iTPC]; par[5] = fBkgRMS[2 + iTPC];
1673 //Mean value and RMS of HIJING Bkg in a 0.3 cone, pt > 2 GeV.
1674 xmax[0] = fJetXMax1[2 + iTPC]; xmax[1] = fJetXMax2[2 + iTPC];
1675 xmin[0] = fJetXMin1[2 + iTPC]; xmin[1] = fJetXMin2[2 + iTPC];
1676 //Factor that multiplies sigma to obtain the best limits,
1677 //by observation, of mono jet ratios (ptjet/ptg) mixed with HIJING Bkg,
1678 //pt_th > 2 GeV, r = 0.3
1679 //X_jet = fJetX1[0]+fJetX2[0]*e_gamma
1681 }//If low pt jet in bkg
1684 //Calculate minimum and maximum limits of the jet ratio.
1685 Double_t min = CalculateJetRatioLimit(ptg, par, xmin);
1686 Double_t max = CalculateJetRatioLimit(ptg, par, xmax);
1688 //Info("IsJetSeleted","%s : Limits min %f, max %f, ptg / ptj %f",
1689 // type.Data(),min,max,ptj/ptg);
1690 if(( min < ptj/ptg ) && ( max > ptj/ptg))
1697 //____________________________________________________________________________
1698 void AliPHOSGammaJet::List() const
1702 Info("List", "%d histograms found", fListHistos->GetEntries() ) ;
1703 TIter next(fListHistos) ;
1705 while ( (h = dynamic_cast<TH2F*>(next())) )
1706 Info("List", "%s", h->GetName()) ;
1709 //____________________________________________________________________________
1710 Double_t AliPHOSGammaJet::MakeEnergy(const Double_t energy)
1712 // Smears the energy according to the energy dependent energy resolution.
1713 // A gaussian distribution is assumed
1715 Double_t sigma = SigmaE(energy) ;
1716 return fRan.Gaus(energy, sigma) ;
1720 //____________________________________________________________________________
1721 void AliPHOSGammaJet::MakeHistos()
1723 // Create histograms to be saved in output file and
1724 // stores them in a TObjectArray
1726 fOutputFile = new TFile(fOutputFileName, "recreate") ;
1728 fListHistos = new TObjArray(10000) ;
1730 // Histos gamma pt vs leading pt
1732 TH1F * hPtSpectra = new TH1F
1733 ("PtSpectra","p_{T i} vs p_{T #gamma}",200,0,200);
1734 hPtSpectra->SetXTitle("p_{T} (GeV/c)");
1735 fListHistos->Add(hPtSpectra) ;
1737 //Histos ratio charged leading pt / gamma pt vs pt
1739 TH2F * hChargeRatio = new TH2F
1740 ("ChargeRatio","p_{T leading charge} /p_{T #gamma} vs p_{T #gamma}",
1742 hChargeRatio->SetYTitle("p_{T lead charge} /p_{T #gamma}");
1743 hChargeRatio->SetXTitle("p_{T #gamma} (GeV/c)");
1744 fListHistos->Add(hChargeRatio) ;
1746 TH2F * hTPCRatio = new TH2F
1747 ("TPCRatio","p_{T leading charge} /p_{T #gamma} vs p_{T #gamma}",
1749 hTPCRatio->SetYTitle("p_{T lead charge} /p_{T #gamma}");
1750 hTPCRatio->SetXTitle("p_{T #gamma} (GeV/c)");
1751 fListHistos->Add(hTPCRatio) ;
1754 TH2F * hPi0Ratio = new TH2F
1755 ("Pi0Ratio","p_{T leading #pi^{0}} /p_{T #gamma} vs p_{T #gamma}",
1757 hPi0Ratio->SetYTitle("p_{T lead #pi^{0}} /p_{T #gamma}");
1758 hPi0Ratio->SetXTitle("p_{T #gamma} (GeV/c)");
1759 fListHistos->Add(hPi0Ratio) ;
1761 TH2F * hPhiGamma = new TH2F
1762 ("PhiGamma","#phi_{#gamma}",200,0,120,200,0,7);
1763 hPhiGamma->SetYTitle("#phi");
1764 hPhiGamma->SetXTitle("p_{T #gamma} (GeV/c)");
1765 fListHistos->Add(hPhiGamma) ;
1767 TH2F * hEtaGamma = new TH2F
1768 ("EtaGamma","#phi_{#gamma}",200,0,120,200,-0.8,0.8);
1769 hEtaGamma->SetYTitle("#eta");
1770 hEtaGamma->SetXTitle("p_{T #gamma} (GeV/c)");
1771 fListHistos->Add(hEtaGamma) ;
1773 // //Jet reconstruction check
1774 // TH2F * hDeltaPhiJet = new TH2F
1775 // ("DeltaPhiJet","#phi_{jet} - #phi_{pyth jet} vs p_{T #gamma}",
1776 // 200,0,120,200,-1.5,1.5);
1777 // hDeltaPhiJet->SetYTitle("#Delta #phi");
1778 // hDeltaPhiJet->SetXTitle("p_{T #gamma} (GeV/c)");
1779 // fListHistos->Add(hDeltaPhiJet) ;
1781 // TH2F * hDeltaPhiTPCJet = new TH2F
1782 // ("DeltaPhiTPCJet","#phi_{jet TPC} - #phi_{pyth jet} vs p_{T #gamma}",
1783 // 200,0,120,200,-1.5,1.5);
1784 // hDeltaPhiTPCJet->SetYTitle("#Delta #phi");
1785 // hDeltaPhiTPCJet->SetXTitle("p_{T #gamma} (GeV/c)");
1786 // fListHistos->Add(hDeltaPhiTPCJet) ;
1788 // TH2F * hDeltaEtaJet = new TH2F
1789 // ("DeltaEtaJet","#phi_{jet} - #phi_{pyth jet} vs p_{T #gamma}",
1790 // 200,0,120,200,-1.5,1.5);
1791 // hDeltaEtaJet->SetYTitle("#Delta #phi");
1792 // hDeltaEtaJet->SetXTitle("p_{T #gamma} (GeV/c)");
1793 // fListHistos->Add(hDeltaEtaJet) ;
1795 // TH2F * hDeltaEtaTPCJet = new TH2F
1796 // ("DeltaEtaTPCJet","#phi_{jet TPC} - #phi_{pyth jet} vs p_{T #gamma}",
1797 // 200,0,120,200,-1.5,1.5);
1798 // hDeltaEtaTPCJet->SetYTitle("#Delta #phi");
1799 // hDeltaEtaTPCJet->SetXTitle("p_{T #gamma} (GeV/c)");
1800 // fListHistos->Add(hDeltaEtaTPCJet) ;
1802 // TH2F * hDeltaPtJet = new TH2F
1803 // ("DeltaPtJet","#phi_{jet} - #phi_{pyth jet} vs p_{T #gamma}",
1804 // 200,0,120,200,0.,100.);
1805 // hDeltaPtJet->SetYTitle("#Delta #phi");
1806 // hDeltaPtJet->SetXTitle("p_{T #gamma} (GeV/c)");
1807 // fListHistos->Add(hDeltaPtJet) ;
1809 // TH2F * hDeltaPtTPCJet = new TH2F
1810 // ("DeltaPtTPCJet","#phi_{jet TPC} - #phi_{pyth jet} vs p_{T #gamma}",
1811 // 200,0,120,200,0.,100.);
1812 // hDeltaPtTPCJet->SetYTitle("#Delta #phi");
1813 // hDeltaPtTPCJet->SetXTitle("p_{T #gamma} (GeV/c)");
1814 // fListHistos->Add(hDeltaPtTPCJet) ;
1817 TH2F * hDeltaPhiCharge = new TH2F
1818 ("DeltaPhiCharge","#phi_{#gamma} - #phi_{charge} vs p_{T #gamma}",
1819 200,0,120,200,0,6.4);
1820 hDeltaPhiCharge->SetYTitle("#Delta #phi");
1821 hDeltaPhiCharge->SetXTitle("p_{T #gamma} (GeV/c)");
1822 fListHistos->Add(hDeltaPhiCharge) ;
1824 TH2F * hDeltaPhiTPC = new TH2F
1825 ("DeltaPhiTPC","#phi_{#gamma} - #phi_{charge} vs p_{T #gamma}",
1826 200,0,120,200,0,6.4);
1827 hDeltaPhiTPC->SetYTitle("#Delta #phi");
1828 hDeltaPhiTPC->SetXTitle("p_{T #gamma} (GeV/c)");
1829 fListHistos->Add(hDeltaPhiTPC) ;
1830 TH2F * hDeltaPhiPi0 = new TH2F
1831 ("DeltaPhiPi0","#phi_{#gamma} - #phi_{ #pi^{0}} vs p_{T #gamma}",
1832 200,0,120,200,0,6.4);
1833 hDeltaPhiPi0->SetYTitle("#Delta #phi");
1834 hDeltaPhiPi0->SetXTitle("p_{T #gamma} (GeV/c)");
1835 fListHistos->Add(hDeltaPhiPi0) ;
1837 TH2F * hDeltaEtaCharge = new TH2F
1838 ("DeltaEtaCharge","#eta_{#gamma} - #eta_{charge} vs p_{T #gamma}",
1839 200,0,120,200,-2,2);
1840 hDeltaEtaCharge->SetYTitle("#Delta #eta");
1841 hDeltaEtaCharge->SetXTitle("p_{T #gamma} (GeV/c)");
1842 fListHistos->Add(hDeltaEtaCharge) ;
1844 TH2F * hDeltaEtaTPC = new TH2F
1845 ("DeltaEtaTPC","#eta_{#gamma} - #eta_{charge} vs p_{T #gamma}",
1846 200,0,120,200,-2,2);
1847 hDeltaEtaTPC->SetYTitle("#Delta #eta");
1848 hDeltaEtaTPC->SetXTitle("p_{T #gamma} (GeV/c)");
1849 fListHistos->Add(hDeltaEtaTPC) ;
1851 TH2F * hDeltaEtaPi0 = new TH2F
1852 ("DeltaEtaPi0","#eta_{#gamma} - #eta_{ #pi^{0}} vs p_{T #gamma}",
1853 200,0,120,200,-2,2);
1854 hDeltaEtaPi0->SetYTitle("#Delta #eta");
1855 hDeltaEtaPi0->SetXTitle("p_{T #gamma} (GeV/c)");
1856 fListHistos->Add(hDeltaEtaPi0) ;
1860 TH2F * hAnglePair = new TH2F
1862 "Angle between #pi^{0} #gamma pair vs p_{T #pi^{0}}",
1863 200,0,50,200,0,0.2);
1864 hAnglePair->SetYTitle("Angle (rad)");
1865 hAnglePair->SetXTitle("E_{ #pi^{0}} (GeV/c)");
1866 fListHistos->Add(hAnglePair) ;
1868 TH2F * hAnglePairAccepted = new TH2F
1869 ("AnglePairAccepted",
1870 "Angle between #pi^{0} #gamma pair vs p_{T #pi^{0}}, both #gamma in eta<0.7, inside window",
1871 200,0,50,200,0,0.2);
1872 hAnglePairAccepted->SetYTitle("Angle (rad)");
1873 hAnglePairAccepted->SetXTitle("E_{ #pi^{0}} (GeV/c)");
1874 fListHistos->Add(hAnglePairAccepted) ;
1876 TH2F * hAnglePairNoCut = new TH2F
1878 "Angle between all #gamma pair vs p_{T #pi^{0}}",200,0,50,200,0,0.2);
1879 hAnglePairNoCut->SetYTitle("Angle (rad)");
1880 hAnglePairNoCut->SetXTitle("E_{ #pi^{0}} (GeV/c)");
1881 fListHistos->Add(hAnglePairNoCut) ;
1883 TH2F * hAnglePairLeadingCut = new TH2F
1884 ("AnglePairLeadingCut",
1885 "Angle between all #gamma pair that have a good phi and pt vs p_{T #pi^{0}}",
1886 200,0,50,200,0,0.2);
1887 hAnglePairLeadingCut->SetYTitle("Angle (rad)");
1888 hAnglePairLeadingCut->SetXTitle("E_{ #pi^{0}} (GeV/c)");
1889 fListHistos->Add(hAnglePairLeadingCut) ;
1891 TH2F * hAnglePairAngleCut = new TH2F
1892 ("AnglePairAngleCut",
1893 "Angle between all #gamma pair (angle + leading cut) vs p_{T #pi^{0}}"
1894 ,200,0,50,200,0,0.2);
1895 hAnglePairAngleCut->SetYTitle("Angle (rad)");
1896 hAnglePairAngleCut->SetXTitle("E_{ #pi^{0}} (GeV/c)");
1897 fListHistos->Add(hAnglePairAngleCut) ;
1899 TH2F * hAnglePairAllCut = new TH2F
1901 "Angle between all #gamma pair (angle + inv mass cut+leading) vs p_{T #pi^{0}}"
1902 ,200,0,50,200,0,0.2);
1903 hAnglePairAllCut->SetYTitle("Angle (rad)");
1904 hAnglePairAllCut->SetXTitle("E_{ #pi^{0}} (GeV/c)");
1905 fListHistos->Add(hAnglePairAllCut) ;
1907 TH2F * hAnglePairLeading = new TH2F
1908 ("AnglePairLeading",
1909 "Angle between all #gamma pair finally selected vs p_{T #pi^{0}}",
1910 200,0,50,200,0,0.2);
1911 hAnglePairLeading->SetYTitle("Angle (rad)");
1912 hAnglePairLeading->SetXTitle("E_{ #pi^{0}} (GeV/c)");
1913 fListHistos->Add(hAnglePairLeading) ;
1916 TH2F * hInvMassPairNoCut = new TH2F
1917 ("InvMassPairNoCut","Invariant Mass of all #gamma pair vs p_{T #gamma}",
1918 120,0,120,360,0,0.5);
1919 hInvMassPairNoCut->SetYTitle("Invariant Mass (GeV/c^{2})");
1920 hInvMassPairNoCut->SetXTitle("p_{T #gamma} (GeV/c)");
1921 fListHistos->Add(hInvMassPairNoCut) ;
1923 TH2F * hInvMassPairLeadingCut = new TH2F
1924 ("InvMassPairLeadingCut",
1925 "Invariant Mass of #gamma pair (leading cuts) vs p_{T #gamma}",
1926 120,0,120,360,0,0.5);
1927 hInvMassPairLeadingCut->SetYTitle("Invariant Mass (GeV/c^{2})");
1928 hInvMassPairLeadingCut->SetXTitle("p_{T #gamma} (GeV/c)");
1929 fListHistos->Add(hInvMassPairLeadingCut) ;
1931 TH2F * hInvMassPairAngleCut = new TH2F
1932 ("InvMassPairAngleCut",
1933 "Invariant Mass of #gamma pair (angle cut) vs p_{T #gamma}",
1934 120,0,120,360,0,0.5);
1935 hInvMassPairAngleCut->SetYTitle("Invariant Mass (GeV/c^{2})");
1936 hInvMassPairAngleCut->SetXTitle("p_{T #gamma} (GeV/c)");
1937 fListHistos->Add(hInvMassPairAngleCut) ;
1940 TH2F * hInvMassPairAllCut = new TH2F
1941 ("InvMassPairAllCut",
1942 "Invariant Mass of #gamma pair (angle+invmass cut+leading) vs p_{T #gamma}",
1943 120,0,120,360,0,0.5);
1944 hInvMassPairAllCut->SetYTitle("Invariant Mass (GeV/c^{2})");
1945 hInvMassPairAllCut->SetXTitle("p_{T #gamma} (GeV/c)");
1946 fListHistos->Add(hInvMassPairAllCut) ;
1948 TH2F * hInvMassPairLeading = new TH2F
1949 ("InvMassPairLeading",
1950 "Invariant Mass of #gamma pair selected vs p_{T #gamma}",
1951 120,0,120,360,0,0.5);
1952 hInvMassPairLeading->SetYTitle("Invariant Mass (GeV/c^{2})");
1953 hInvMassPairLeading->SetXTitle("p_{T #gamma} (GeV/c)");
1954 fListHistos->Add(hInvMassPairLeading) ;
1959 TH1F * hNGamma = new TH1F("NGamma","Number of #gamma over PHOS",240,0,120);
1960 hNGamma->SetYTitle("N");
1961 hNGamma->SetXTitle("p_{T #gamma}(GeV/c)");
1962 fListHistos->Add(hNGamma) ;
1964 TH1F * hNBkg = new TH1F("NBkg","bkg multiplicity",9000,0,9000);
1965 hNBkg->SetYTitle("counts");
1966 hNBkg->SetXTitle("N");
1967 fListHistos->Add(hNBkg) ;
1969 TH2F * hNLeading = new TH2F
1970 ("NLeading","Accepted Jet Leading", 240,0,120,240,0,120);
1971 hNLeading->SetYTitle("p_{T charge} (GeV/c)");
1972 hNLeading->SetXTitle("p_{T #gamma}(GeV/c)");
1973 fListHistos->Add(hNLeading) ;
1976 TH1F * hN = new TH1F("NJet","Accepted jets",240,0,120);
1978 hN->SetXTitle("p_{T #gamma}(GeV/c)");
1979 fListHistos->Add(hN) ;
1982 //Ratios and Pt dist of reconstructed (not selected) jets
1984 TH2F * hJetRatio = new TH2F
1985 ("JetRatio","p_{T jet lead}/p_{T #gamma} vs p_{T #gamma}",
1986 240,0,120,200,0,10);
1987 hJetRatio->SetYTitle("p_{T jet lead #pi^{0}}/p_{T #gamma}");
1988 hJetRatio->SetXTitle("p_{T #gamma} (GeV/c)");
1989 fListHistos->Add(hJetRatio) ;
1992 TH2F * hJetPt = new TH2F
1993 ("JetPt", "p_{T jet lead} vs p_{T #gamma}",240,0,120,400,0,200);
1994 hJetPt->SetYTitle("p_{T jet lead #pi^{0}}/p_{T #gamma}");
1995 hJetPt->SetXTitle("p_{T #gamma} (GeV/c)");
1996 fListHistos->Add(hJetPt) ;
2001 TH2F * hBkgRatio = new TH2F
2002 ("BkgRatio","p_{T bkg lead}/p_{T #gamma} vs p_{T #gamma}",
2003 240,0,120,200,0,10);
2004 hBkgRatio->SetYTitle("p_{T bkg lead charge}/p_{T #gamma}");
2005 hBkgRatio->SetXTitle("p_{T #gamma} (GeV/c)");
2006 fListHistos->Add(hBkgRatio) ;
2009 TH2F * hBkgPt = new TH2F
2010 ("BkgPt","p_{T jet lead} vs p_{T #gamma}",240,0,120,400,0,200);
2011 hBkgPt->SetYTitle("p_{T jet lead charge}/p_{T #gamma}");
2012 hBkgPt->SetXTitle("p_{T #gamma} (GeV/c)");
2013 fListHistos->Add(hBkgPt) ;
2018 TH2F * hJetFragment =
2019 new TH2F("JetFragment","x = p_{T i charged}/p_{T #gamma}",
2020 240,0.,120.,1000,0.,1.2);
2021 hJetFragment->SetYTitle("x_{T}");
2022 hJetFragment->SetXTitle("p_{T #gamma}");
2023 fListHistos->Add(hJetFragment) ;
2025 TH2F * hBkgFragment = new TH2F
2026 ("BkgFragment","x = p_{T i charged}/p_{T #gamma}",
2027 240,0.,120.,1000,0.,1.2);
2028 hBkgFragment->SetYTitle("x_{T}");
2029 hBkgFragment->SetXTitle("p_{T #gamma}");
2030 fListHistos->Add(hBkgFragment) ;
2033 new TH2F("JetPtDist","x = p_{T i charged}",240,0.,120.,400,0.,200.);
2034 hJetPtDist->SetXTitle("p_{T #gamma} (GeV/c)");
2035 fListHistos->Add(hJetPtDist) ;
2037 TH2F * hBkgPtDist = new TH2F
2038 ("BkgPtDist","x = p_{T i charged}",240,0.,120.,400,0.,200.);
2039 hBkgPtDist->SetXTitle("p_{T #gamma} (GeV/c)");
2040 fListHistos->Add(hBkgPtDist) ;
2045 TH1F * hNBkgTPC = new TH1F
2046 ("NBkgTPC","TPC bkg multiplicity ",9000,0,9000);
2047 hNBkgTPC->SetYTitle("counts");
2048 hNBkgTPC->SetXTitle("N");
2049 fListHistos->Add(hNBkgTPC) ;
2051 TH2F * hNTPCLeading = new TH2F
2052 ("NTPCLeading","Accepted TPC jet leading",240,0,120,240,0,120);
2053 hNTPCLeading->SetYTitle("p_{T charge} (GeV/c)");
2054 hNTPCLeading->SetXTitle("p_{T #gamma}(GeV/c)");
2055 fListHistos->Add(hNTPCLeading) ;
2057 TH1F * hNTPC = new TH1F("NTPCJet","Number of TPC jets",240,0,120);
2058 hNTPC->SetYTitle("N");
2059 hNTPC->SetXTitle("p_{T #gamma}(GeV/c)");
2060 fListHistos->Add(hNTPC) ;
2062 TH2F * hJetTPCRatio = new TH2F
2063 ("JetTPCRatio", "p_{T jet lead TPC}/p_{T #gamma} vs p_{T #gamma}",
2064 240,0,120,200,0,10);
2065 hJetTPCRatio->SetYTitle("p_{T jet lead TPC}/p_{T #gamma}");
2066 hJetTPCRatio->SetXTitle("p_{T #gamma} (GeV/c)");
2067 fListHistos->Add(hJetTPCRatio) ;
2069 TH2F * hBkgTPCRatio = new TH2F
2070 ("BkgTPCRatio","p_{T bkg lead TPC}/p_{T #gamma} vs p_{T #gamma}",
2071 240,0,120,200,0,10);
2072 hBkgTPCRatio->SetYTitle("p_{T bkg lead TPC}/p_{T #gamma}");
2073 hBkgTPCRatio->SetXTitle("p_{T #gamma} (GeV/c)");
2074 fListHistos->Add(hBkgTPCRatio) ;
2076 TH2F * hJetTPCPt = new TH2F
2077 ("JetTPCPt", "p_{T jet lead TPC} vs p_{T #gamma}",240,0,120,400,0,200);
2078 hJetTPCPt->SetYTitle("p_{T jet lead TPC}/p_{T #gamma}");
2079 hJetTPCPt->SetXTitle("p_{T #gamma} (GeV/c)");
2080 fListHistos->Add(hJetTPCPt) ;
2082 TH2F * hBkgTPCPt = new TH2F
2083 ("BkgTPCPt", "p_{T bkg lead TPC} vs p_{T #gamma}",240,0,120,400,0,200);
2084 hBkgTPCPt->SetYTitle("p_{T bkg lead TPC}/p_{T #gamma}");
2085 hBkgTPCPt->SetXTitle("p_{T #gamma} (GeV/c)");
2086 fListHistos->Add(hBkgTPCPt) ;
2090 TH2F * hJetTPCFragment =
2091 new TH2F("JetTPCFragment","x = p_{T i charged}/p_{T #gamma}",
2092 240,0.,120.,1000,0.,1.2);
2093 hJetTPCFragment->SetYTitle("x_{T}");
2094 hJetTPCFragment->SetXTitle("p_{T #gamma}");
2095 fListHistos->Add(hJetTPCFragment) ;
2097 TH2F * hBkgTPCFragment = new TH2F
2098 ("BkgTPCFragment","x = p_{T i charged}/p_{T #gamma}",
2099 240,0.,120.,1000,0.,1.2);
2100 hBkgTPCFragment->SetYTitle("x_{T}");
2101 hBkgTPCFragment->SetXTitle("p_{T #gamma}");
2102 fListHistos->Add(hBkgTPCFragment) ;
2105 TH2F * hJetTPCPtDist = new TH2F("JetTPCPtDist",
2106 "x = p_{T i charged}",240,0.,120.,400,0.,200.);
2107 hJetTPCPtDist->SetXTitle("p_{T #gamma} (GeV/c)");
2108 fListHistos->Add(hJetTPCPtDist) ;
2110 TH2F * hBkgTPCPtDist = new TH2F
2111 ("BkgTPCPtDist","x = p_{T i charged}",240,0.,120.,400,0.,200.);
2112 hBkgTPCPtDist->SetXTitle("p_{T #gamma} (GeV/c)");
2113 fListHistos->Add(hBkgTPCPtDist) ;
2119 //If we want to study the jet for different cones and pt. Old version
2121 TH2F * hJetRatios[5][5];
2122 TH2F * hJetTPCRatios[5][5];
2124 TH2F * hJetPts[5][5];
2125 TH2F * hJetTPCPts[5][5];
2127 TH2F * hBkgRatios[5][5];
2128 TH2F * hBkgTPCRatios[5][5];
2130 TH2F * hBkgPts[5][5];
2131 TH2F * hBkgTPCPts[5][5];
2133 TH2F * hNLeadings[5][5];
2134 TH2F * hNTPCLeadings[5][5];
2137 TH1F * hNTPCs[5][5];
2139 TH1F * hNBkgs[5][5];
2140 TH1F * hNBkgTPCs[5][5];
2142 TH2F * hJetFragments[5][5];
2143 TH2F * hBkgFragments[5][5];
2144 TH2F * hJetPtDists[5][5];
2145 TH2F * hBkgPtDists[5][5];
2147 TH2F * hJetTPCFragments[5][5];
2148 TH2F * hBkgTPCFragments[5][5];
2149 TH2F * hJetTPCPtDists[5][5];
2150 TH2F * hBkgTPCPtDists[5][5];
2153 for(Int_t icone = 0; icone<fNCone; icone++){
2154 for(Int_t ipt = 0; ipt<fNPt;ipt++){
2159 hJetRatios[icone][ipt] = new TH2F
2160 ("JetRatioCone"+fNameCones[icone]+"Pt"+fNamePtThres[ipt],
2161 "p_{T jet lead #pi^{0}}/p_{T #gamma} vs p_{T #gamma}, cone ="
2162 +fNameCones[icone]+", pt>" +fNamePtThres[ipt]+" GeV/c",
2163 240,0,120,200,0,10);
2164 hJetRatios[icone][ipt]->
2165 SetYTitle("p_{T jet lead #pi^{0}}/p_{T #gamma}");
2166 hJetRatios[icone][ipt]->SetXTitle("p_{T #gamma} (GeV/c)");
2167 fListHistos->Add(hJetRatios[icone][ipt]) ;
2170 hJetPts[icone][ipt] = new TH2F
2171 ("JetPtCone"+fNameCones[icone]+"Pt"+fNamePtThres[ipt],
2172 "p_{T jet lead #pi^{0}}/p_{T #gamma} vs p_{T #gamma}, cone ="
2173 +fNameCones[icone]+", pt>" +fNamePtThres[ipt]+" GeV/c",
2174 240,0,120,400,0,200);
2175 hJetPts[icone][ipt]->
2176 SetYTitle("p_{T jet lead #pi^{0}}/p_{T #gamma}");
2177 hJetPts[icone][ipt]->SetXTitle("p_{T #gamma} (GeV/c)");
2178 fListHistos->Add(hJetPts[icone][ipt]) ;
2182 hBkgRatios[icone][ipt] = new TH2F
2183 ("BkgRatioCone"+fNameCones[icone]+"Pt"+fNamePtThres[ipt],
2184 "p_{T bkg lead #pi^{0}}/p_{T #gamma} vs p_{T #gamma}, cone ="
2185 +fNameCones[icone]+", pt>" +fNamePtThres[ipt]+" GeV/c",
2186 240,0,120,200,0,10);
2187 hBkgRatios[icone][ipt]->
2188 SetYTitle("p_{T bkg lead #pi^{0}}/p_{T #gamma}");
2189 hBkgRatios[icone][ipt]->SetXTitle("p_{T #gamma} (GeV/c)");
2190 fListHistos->Add(hBkgRatios[icone][ipt]) ;
2194 hBkgPts[icone][ipt] = new TH2F
2195 ("BkgPtCone"+fNameCones[icone]+"Pt"+fNamePtThres[ipt],
2196 "p_{T jet lead #pi^{0}}/p_{T #gamma} vs p_{T #gamma}, cone ="
2197 +fNameCones[icone]+", pt>" +fNamePtThres[ipt]+" GeV/c",
2198 240,0,120,400,0,200);
2199 hBkgPts[icone][ipt]->
2200 SetYTitle("p_{T jet lead #pi^{0}}/p_{T #gamma}");
2201 hBkgPts[icone][ipt]->SetXTitle("p_{T #gamma} (GeV/c)");
2202 fListHistos->Add(hBkgPts[icone][ipt]) ;
2207 hJetTPCRatios[icone][ipt] = new TH2F
2208 ("JetTPCRatioCone"+fNameCones[icone]+"Pt"+fNamePtThres[ipt],
2209 "p_{T jet lead TPC}/p_{T #gamma} vs p_{T #gamma}, cone ="
2210 +fNameCones[icone]+", pt>" +fNamePtThres[ipt]+" GeV/c",
2211 240,0,120,200,0,10);
2212 hJetTPCRatios[icone][ipt]->SetYTitle("p_{T jet lead TPC}/p_{T #gamma}");
2213 hJetTPCRatios[icone][ipt]->SetXTitle("p_{T #gamma} (GeV/c)");
2214 fListHistos->Add(hJetTPCRatios[icone][ipt]) ;
2216 hBkgTPCRatios[icone][ipt] = new TH2F
2217 ("BkgTPCRatioCone"+fNameCones[icone]+"Pt"+fNamePtThres[ipt],
2218 "p_{T bkg lead TPC}/p_{T #gamma} vs p_{T #gamma}, cone ="
2219 +fNameCones[icone]+", pt>" +fNamePtThres[ipt]+" GeV/c",
2220 240,0,120,200,0,10);
2221 hBkgTPCRatios[icone][ipt]->SetYTitle("p_{T bkg lead TPC}/p_{T #gamma}");
2222 hBkgTPCRatios[icone][ipt]->SetXTitle("p_{T #gamma} (GeV/c)");
2223 fListHistos->Add(hBkgTPCRatios[icone][ipt]) ;
2225 hJetTPCPts[icone][ipt] = new TH2F
2226 ("JetTPCPtCone"+fNameCones[icone]+"Pt"+fNamePtThres[ipt],
2227 "p_{T jet lead TPC}/p_{T #gamma} vs p_{T #gamma}, cone ="
2228 +fNameCones[icone]+", pt>" +fNamePtThres[ipt]+" GeV/c",
2229 240,0,120,400,0,200);
2230 hJetTPCPts[icone][ipt]->SetYTitle("p_{T jet lead TPC}/p_{T #gamma}");
2231 hJetTPCPts[icone][ipt]->SetXTitle("p_{T #gamma} (GeV/c)");
2232 fListHistos->Add(hJetTPCPts[icone][ipt]) ;
2234 hBkgTPCPts[icone][ipt] = new TH2F
2235 ("BkgTPCPtCone"+fNameCones[icone]+"Pt"+fNamePtThres[ipt],
2236 "p_{T bkg lead TPC}/p_{T #gamma} vs p_{T #gamma}, cone ="
2237 +fNameCones[icone]+", pt>" +fNamePtThres[ipt]+" GeV/c",
2238 240,0,120,400,0,200);
2239 hBkgTPCPts[icone][ipt]->SetYTitle("p_{T bkg lead TPC}/p_{T #gamma}");
2240 hBkgTPCPts[icone][ipt]->SetXTitle("p_{T #gamma} (GeV/c)");
2241 fListHistos->Add(hBkgTPCPts[icone][ipt]) ;
2246 hNBkgs[icone][ipt] = new TH1F
2247 ("NBkgCone"+fNameCones[icone]+"Pt"+fNamePtThres[ipt],
2248 "bkg multiplicity cone ="+fNameCones[icone]+", pt>"
2249 +fNamePtThres[ipt]+" GeV/c",9000,0,9000);
2250 hNBkgs[icone][ipt]->SetYTitle("counts");
2251 hNBkgs[icone][ipt]->SetXTitle("N");
2252 fListHistos->Add(hNBkgs[icone][ipt]) ;
2254 hNBkgTPCs[icone][ipt] = new TH1F
2255 ("NBkgTPCCone"+fNameCones[icone]+"Pt"+fNamePtThres[ipt],
2256 "bkg multiplicity cone ="+fNameCones[icone]+", pt>"
2257 +fNamePtThres[ipt]+" GeV/c",9000,0,9000);
2258 hNBkgTPCs[icone][ipt]->SetYTitle("counts");
2259 hNBkgTPCs[icone][ipt]->SetXTitle("N");
2260 fListHistos->Add(hNBkgTPCs[icone][ipt]) ;
2263 hNLeadings[icone][ipt] = new TH2F
2264 ("NLeadingCone"+fNameCones[icone]+"Pt"+fNamePtThres[ipt],
2265 "p_{T #gamma} vs p_{T #pi^{0}} cone ="+fNameCones[icone]+", pt>"
2266 +fNamePtThres[ipt]+" GeV/c",120,0,120,120,0,120);
2267 hNLeadings[icone][ipt]->SetYTitle("p_{T #pi^{0}}(GeV/c)");
2268 hNLeadings[icone][ipt]->SetXTitle("p_{T #gamma}(GeV/c)");
2269 fListHistos->Add(hNLeadings[icone][ipt]) ;
2271 hNs[icone][ipt] = new TH1F
2272 ("NJetCone"+fNameCones[icone]+"Pt"+fNamePtThres[ipt],
2273 "Number of neutral jets, cone ="+fNameCones[icone]+", pt>"
2274 +fNamePtThres[ipt]+" GeV/c",120,0,120);
2275 hNs[icone][ipt]->SetYTitle("N");
2276 hNs[icone][ipt]->SetXTitle("p_{T #gamma}(GeV/c)");
2277 fListHistos->Add(hNs[icone][ipt]) ;
2279 //Fragmentation Function
2280 hJetFragments[icone][ipt] = new TH2F
2281 ("JetFragmentCone"+fNameCones[icone]+"Pt"+fNamePtThres[ipt],
2282 "x_{T} = p_{T i}/p_{T #gamma}, cone ="+fNameCones[icone]+", pt>"
2283 +fNamePtThres[ipt]+" GeV/c",120,0.,120.,240,0.,1.2);
2284 hJetFragments[icone][ipt]->SetYTitle("x_{T}");
2285 hJetFragments[icone][ipt]->SetXTitle("p_{T #gamma}");
2286 fListHistos->Add(hJetFragments[icone][ipt]) ;
2288 hBkgFragments[icone][ipt] = new TH2F
2289 ("BkgFragmentCone"+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 hBkgFragments[icone][ipt]->SetYTitle("x_{T}");
2293 hBkgFragments[icone][ipt]->SetXTitle("p_{T #gamma}");
2294 fListHistos->Add(hBkgFragments[icone][ipt]) ;
2297 //Jet particle distribution
2299 hJetPtDists[icone][ipt] = new TH2F
2300 ("JetPtDistCone"+fNameCones[icone]+"Pt"+fNamePtThres[ipt],
2301 "p_{T i}, cone ="+fNameCones[icone]+", pt>" +fNamePtThres[ipt]+
2302 " GeV/c",120,0.,120.,120,0.,120.);
2303 hJetPtDists[icone][ipt]->SetXTitle("p_{T #gamma} (GeV/c)");
2304 fListHistos->Add(hJetPtDists[icone][ipt]) ;
2306 hBkgPtDists[icone][ipt] = new TH2F
2307 ("BkgPtDistCone"+fNameCones[icone]+"Pt"+fNamePtThres[ipt],
2308 "p_{T i}, cone ="+fNameCones[icone]+", pt>" +fNamePtThres[ipt]+
2309 " GeV/c",120,0.,120.,120,0.,120.);
2310 hBkgPtDists[icone][ipt]->SetXTitle("p_{T #gamma} (GeV/c)");
2311 fListHistos->Add(hBkgPtDists[icone][ipt]) ;
2315 hNTPCLeadings[icone][ipt] = new TH2F
2316 ("NTPCLeadingCone"+fNameCones[icone]+"Pt"+fNamePtThres[ipt],
2317 "p_{T #gamma} vs p_{T charge} cone ="+fNameCones[icone]+", pt>"
2318 +fNamePtThres[ipt]+" GeV/c",120,0,120,120,0,120);
2319 hNTPCLeadings[icone][ipt]->SetYTitle("p_{T charge} (GeV/c)");
2320 hNTPCLeadings[icone][ipt]->SetXTitle("p_{T #gamma}(GeV/c)");
2321 fListHistos->Add(hNTPCLeadings[icone][ipt]) ;
2323 hNTPCs[icone][ipt] = new TH1F
2324 ("NTPCJetCone"+fNameCones[icone]+"Pt"+fNamePtThres[ipt],
2325 "Number of charged jets, cone ="+fNameCones[icone]+", pt>"
2326 +fNamePtThres[ipt]+" GeV/c",120,0,120);
2327 hNTPCs[icone][ipt]->SetYTitle("N");
2328 hNTPCs[icone][ipt]->SetXTitle("p_{T #gamma}(GeV/c)");
2329 fListHistos->Add(hNTPCs[icone][ipt]) ;
2331 hJetTPCFragments[icone][ipt] = new TH2F
2332 ("JetTPCFragmentCone"+fNameCones[icone]+"Pt"+fNamePtThres[ipt],
2333 "x = p_{T i charged}/p_{T #gamma}, cone ="+fNameCones[icone]+", pt>"
2334 +fNamePtThres[ipt]+" GeV/c",120,0.,120.,240,0.,1.2);
2335 hJetTPCFragments[icone][ipt]->SetYTitle("x_{T}");
2336 hJetTPCFragments[icone][ipt]->SetXTitle("p_{T #gamma}");
2337 fListHistos->Add(hJetTPCFragments[icone][ipt]) ;
2339 hBkgTPCFragments[icone][ipt] = new TH2F
2340 ("BkgTPCFragmentCone"+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 hBkgTPCFragments[icone][ipt]->SetYTitle("x_{T}");
2344 hBkgTPCFragments[icone][ipt]->SetXTitle("p_{T #gamma}");
2345 fListHistos->Add(hBkgTPCFragments[icone][ipt]) ;
2347 hJetTPCPtDists[icone][ipt] = new TH2F
2348 ("JetTPCPtDistCone"+fNameCones[icone]+"Pt"+fNamePtThres[ipt],
2349 "x = p_{T i charged}, cone ="+fNameCones[icone]+", pt>"
2350 +fNamePtThres[ipt]+" GeV/c",120,0.,120.,120,0.,120.);
2351 hJetTPCPtDists[icone][ipt]->SetXTitle("p_{T #gamma} (GeV/c)");
2352 fListHistos->Add(hJetTPCPtDists[icone][ipt]) ;
2354 hBkgTPCPtDists[icone][ipt] = new TH2F
2355 ("BkgTPCPtDistCone"+fNameCones[icone]+"Pt"+fNamePtThres[ipt],
2356 "x = p_{T i charged}, cone ="+fNameCones[icone]+", pt>" +
2357 fNamePtThres[ipt]+" GeV/c",120,0.,120.,120,0.,120.);
2358 hBkgTPCPtDists[icone][ipt]->SetXTitle("p_{T #gamma} (GeV/c)");
2359 fListHistos->Add(hBkgTPCPtDists[icone][ipt]) ;
2363 }//If we want to study any cone or pt threshold
2367 //____________________________________________________________________________
2368 void AliPHOSGammaJet::MakeJet(TClonesArray * pl,
2369 Double_t ptg, Double_t phig,
2370 Double_t ptl, Double_t phil, Double_t etal,
2371 TString conf, TLorentzVector & jet)
2373 //Fill the jet with the particles around the leading particle with
2374 //R=fCone and pt_th = fPtThres. Calculate the energy of the jet and
2375 //check if we select it. Fill jet histograms
2376 Float_t ptcut = 0. ;
2378 if(ptg > fPtJetSelectionCut) ptcut = 2. ;
2382 TClonesArray * jetList = new TClonesArray("TParticle",1000);
2383 TClonesArray * bkgList = new TClonesArray("TParticle",1000);
2385 TLorentzVector bkg(0,0,0,0);
2386 TLorentzVector lv (0,0,0,0);
2388 Double_t ptjet = 0.0;
2389 Double_t ptbkg = 0.0;
2396 TParticle * particle = 0 ;
2398 while ( (particle = dynamic_cast<TParticle*>(next())) ) {
2404 SetJet(particle, b0, fCone, etal, phil) ;
2407 new((*jetList)[n0++]) TParticle(*particle) ;
2408 particle->Momentum(lv);
2409 if(particle->Pt() > ptcut ){
2411 ptjet+=particle->Pt();
2415 //Background around (phi_gamma-pi, eta_leading)
2416 SetJet(particle, b1, fCone,etal, phig) ;
2419 new((*bkgList)[n1++]) TParticle(*particle) ;
2420 if(particle->Pt() > ptcut ){
2422 ptbkg+=particle->Pt();
2430 if(strstr(fOptionGJ,"deb") || strstr(fOptionGJ,"deb all"))
2431 Info("MakeJet","Gamma pt %f, Jet pt %f, Bkg pt %f",ptg,ptjet,ptbkg);
2436 Double_t rat = ptjet/ptg ;
2437 Double_t ratbkg = ptbkg/ptg ;
2440 (fListHistos->FindObject("Jet"+conf+"Ratio"))->Fill(ptg,rat);
2442 (fListHistos->FindObject("Jet"+conf+"Pt")) ->Fill(ptg,ptjet);
2444 (fListHistos->FindObject("Bkg"+conf+"Ratio"))->Fill(ptg,ratbkg);
2446 (fListHistos->FindObject("Bkg"+conf+"Pt")) ->Fill(ptg,ptbkg);
2449 if(IsJetSelected(ptg,ptjet,conf) || fSelect){
2450 if(strstr(fOptionGJ,"deb") || strstr(fOptionGJ,"deb all"))
2451 Info("MakeJet","JetSelected");
2452 dynamic_cast<TH1F*>(fListHistos->FindObject("N"+conf+"Jet"))->
2454 dynamic_cast<TH2F*>(fListHistos->FindObject("N"+conf+"Leading"))
2456 FillJetHistos(jetList, ptg, conf, "Jet");
2457 FillJetHistos(bkgList, ptg, conf, "Bkg");
2464 //____________________________________________________________________________
2465 void AliPHOSGammaJet::MakeJetAnyConeOrPt(TClonesArray * pl, Double_t ptg,
2466 Double_t phig, Double_t ptl,
2467 Double_t phil, Double_t etal,
2470 //Fill the jet with the particles around the leading particle with
2471 //R=fCone(i) and pt_th = fPtThres(i). Calculate the energy of the jet and
2472 //check if we select it. Fill jet i histograms
2474 TClonesArray * jetList = new TClonesArray("TParticle",1000);
2475 TClonesArray * bkgList = new TClonesArray("TParticle",1000);
2477 Double_t ptjet = 0.0;
2478 Double_t ptbkg = 0.0;
2485 //Create as many jets as cones and pt thresholds are defined
2486 Double_t maxcut = fJetRatioMaxCut;
2487 Double_t mincut = fJetRatioMinCut;
2490 maxcut = fJetTPCRatioMaxCut;
2491 mincut = fJetTPCRatioMinCut;
2494 Double_t ratjet = 0;
2495 Double_t ratbkg = 0;
2497 for(Int_t icone = 0; icone<fNCone; icone++) {
2498 for(Int_t ipt = 0; ipt<fNPt;ipt++) {
2500 TString cone = fNameCones[icone] ;
2501 TString ptcut = fNamePtThres[ipt] ;
2504 TParticle * particle = 0 ;
2509 while ( (particle = dynamic_cast<TParticle*>(next())) ) {
2513 SetJet(particle, b0, fCones[icone],etal, phil) ;
2514 SetJet(particle, b1, fCones[icone],etal, phig) ;
2517 new((*jetList)[n0++]) TParticle(*particle) ;
2518 if(particle->Pt() > fPtThres[ipt] )
2519 ptjet+=particle->Pt();
2522 new((*bkgList)[n1++]) TParticle(*particle) ;
2523 if(particle->Pt() > fPtThres[ipt] )
2524 ptbkg+=particle->Pt();
2532 if(strstr(fOptionGJ,"deb")){
2533 Info("MakeJetAnyPt","cone %f, ptcut %f",fCones[icone],fPtThres[ipt]);
2534 Info("MakeJetAnyPt","pT: Gamma %f, Jet %f, Bkg %f",ptg,ptjet,ptbkg);
2537 ratjet = ptjet /ptg;
2541 (fListHistos->FindObject("Jet"+conf+"RatioCone"+cone+"Pt"+ptcut))
2544 (fListHistos->FindObject("Jet"+conf+"PtCone"+cone+"Pt"+ptcut))
2548 (fListHistos->FindObject("Bkg"+conf+"RatioCone"+cone+"Pt"+ptcut))
2552 (fListHistos->FindObject("Bkg"+conf+"PtCone"+cone+"Pt"+ptcut))
2557 if((ratjet < maxcut) && (ratjet > mincut)){
2560 (fListHistos->FindObject("N"+conf+"JetCone"+cone+"Pt"+ptcut))->
2563 (fListHistos->FindObject("N"+conf+"LeadingCone"+cone+"Pt"+ptcut))
2566 FillJetHistosAnyConeOrPt
2567 (jetList,ptg,conf,"Jet",fNameCones[icone],fNamePtThres[ipt]);
2568 FillJetHistosAnyConeOrPt
2569 (bkgList,ptg,conf,"Bkg",fNameCones[icone],fNamePtThres[ipt]);
2579 //____________________________________________________________________________
2580 void AliPHOSGammaJet::MakePhoton(TLorentzVector & particle)
2582 //Fast reconstruction for photons
2583 Double_t energy = particle.E() ;
2584 Double_t modenergy = MakeEnergy(energy) ;
2585 //Info("MakePhoton","Energy %f, Modif %f",energy,modenergy);
2587 // get the detected direction
2588 TVector3 pos = particle.Vect();
2590 TVector3 modpos = MakePosition(energy, pos) ;
2591 modpos *= modenergy / 460.;
2593 Float_t modtheta = modpos.Theta();
2594 Float_t modphi = modpos.Phi();
2596 // Set the modified 4-momentum of the reconstructed particle
2597 Float_t py = modenergy*TMath::Sin(modphi)*TMath::Sin(modtheta);
2598 Float_t px = modenergy*TMath::Cos(modphi)*TMath::Sin(modtheta);
2599 Float_t pz = modenergy*TMath::Cos(modtheta);
2601 particle.SetPxPyPzE(px,py,pz,modenergy);
2605 //____________________________________________________________________________
2606 TVector3 AliPHOSGammaJet::MakePosition(const Double_t energy, const TVector3 pos)
2608 // Smears the impact position according to the energy dependent position resolution
2609 // A gaussian position distribution is assumed
2613 Double_t sigma = SigmaP(energy) ;
2614 Double_t x = fRan.Gaus( pos.X(), sigma ) ;
2615 Double_t z = fRan.Gaus( pos.Z(), sigma ) ;
2616 Double_t y = pos.Y() ;
2622 // Info("MakePosition","Theta dif %f",pos.Theta()-newpos.Theta());
2623 // Info("MakePosition","Phi dif %f",pos.Phi()-newpos.Phi());
2627 //____________________________________________________________________________
2628 void AliPHOSGammaJet::Pi0Decay(Double_t mPi0, TLorentzVector &p0,
2629 TLorentzVector &p1, TLorentzVector &p2, Double_t &angle) {
2630 // Perform isotropic decay pi0 -> 2 photons
2631 // p0 is pi0 4-momentum (inut)
2632 // p1 and p2 are photon 4-momenta (output)
2633 // cout<<"Boost vector"<<endl;
2634 TVector3 b = p0.BoostVector();
2635 //cout<<"Parameters"<<endl;
2636 //Double_t mPi0 = p0.M();
2637 Double_t phi = TMath::TwoPi() * gRandom->Rndm();
2638 Double_t cosThe = 2 * gRandom->Rndm() - 1;
2639 Double_t cosPhi = TMath::Cos(phi);
2640 Double_t sinPhi = TMath::Sin(phi);
2641 Double_t sinThe = TMath::Sqrt(1-cosThe*cosThe);
2642 Double_t ePi0 = mPi0/2.;
2643 //cout<<"ePi0 "<<ePi0<<endl;
2644 //cout<<"Components"<<endl;
2645 p1.SetPx(+ePi0*cosPhi*sinThe);
2646 p1.SetPy(+ePi0*sinPhi*sinThe);
2647 p1.SetPz(+ePi0*cosThe);
2649 //cout<<"p1: "<<p1.Px()<<" "<<p1.Py()<<" "<<p1.Pz()<<" "<<p1.E()<<endl;
2650 //cout<<"p1 Mass: "<<p1.Px()*p1.Px()+p1.Py()*p1.Py()+p1.Pz()*p1.Pz()-p1.E()*p1.E()<<endl;
2651 p2.SetPx(-ePi0*cosPhi*sinThe);
2652 p2.SetPy(-ePi0*sinPhi*sinThe);
2653 p2.SetPz(-ePi0*cosThe);
2655 //cout<<"p2: "<<p2.Px()<<" "<<p2.Py()<<" "<<p2.Pz()<<" "<<p2.E()<<endl;
2656 //cout<<"p2 Mass: "<<p2.Px()*p2.Px()+p2.Py()*p2.Py()+p2.Pz()*p2.Pz()-p2.E()*p2.E()<<endl;
2657 //cout<<"Boost "<<b.X()<<" "<<b.Y()<<" "<<b.Z()<<endl;
2659 //cout<<"p1: "<<p1.Px()<<" "<<p1.Py()<<" "<<p1.Pz()<<" "<<p1.E()<<endl;
2661 //cout<<"p2: "<<p2.Px()<<" "<<p2.Py()<<" "<<p2.Pz()<<" "<<p2.E()<<endl;
2662 //cout<<"angle"<<endl;
2663 angle = p1.Angle(p2.Vect());
2664 //cout<<angle<<endl;
2667 //____________________________________________________________________________
2668 void AliPHOSGammaJet::Plot(TString what, Option_t * option) const
2670 //Plot some relevant histograms of the analysis
2671 TH2F * h = dynamic_cast<TH2F*>(fOutputFile->Get(what));
2675 else if (what == "all") {
2676 TCanvas * can = new TCanvas("GammaJet", "Gamma-Jet Study",10,40,1600,1200);
2678 can->Range(0,0,22,20);
2679 TPaveLabel *pl1 = new TPaveLabel(1,18,20,19.5,"Titre","br");
2680 pl1->SetFillColor(18);
2681 pl1->SetTextFont(32);
2682 pl1->SetTextColor(49);
2686 Float_t begx = -0.29, begy = 0., endx = 0., endy = 0.30 ;
2687 for (index = 0 ; index < fListHistos->GetEntries() ; index++) {
2688 TString name("pad") ;
2692 if (begx >= 1.0 || endx >= 1.0) {
2698 printf("%f %f %f %f \n", begx, begy, endx, endy) ;
2699 pad = new TPad(name,"This is a pad",begx,begy,endx,endy,33);
2702 fListHistos->At(index)->Draw(option) ;
2710 Info("Draw", "Histogram %s does not exist or unknown option", what.Data()) ;
2715 //____________________________________________________________________________
2716 void AliPHOSGammaJet::Print(const Option_t * opt) const
2719 //Print some relevant parameters set for the analysis
2723 Info("Print", "%s %s", GetName(), GetTitle() ) ;
2725 printf("Eta cut : %f\n", fEtaCut) ;
2726 printf("D phi max cut : %f\n", fPhiMaxCut) ;
2727 printf("D phi min cut : %f\n", fPhiMinCut) ;
2728 printf("Leading Ratio max cut : %f\n", fRatioMaxCut) ;
2729 printf("Leading Ratio min cut : %f\n", fRatioMinCut) ;
2730 printf("Jet Ratio max cut : %f\n", fJetRatioMaxCut) ;
2731 printf("Jet Ratio min cut : %f\n", fJetRatioMinCut) ;
2732 printf("Jet TPC Ratio max cut : %f\n", fJetTPCRatioMaxCut) ;
2733 printf("Jet TPC Ratio min cut : %f\n", fJetTPCRatioMinCut) ;
2734 printf("Fast recons : %d\n", fOptFast);
2735 printf("Inv Mass max cut : %f\n", fInvMassMaxCut) ;
2736 printf("Inv Mass min cut : %f\n", fInvMassMinCut) ;
2740 //__________________________________________________________________________
2741 TChain * AliPHOSGammaJet::ReadESDfromdisk(const UInt_t eventsToRead,
2742 const TString dirName,
2743 const TString esdTreeName,
2744 const char * pattern)
2746 // Reads ESDs from Disk
2749 AliInfo( Form("\nReading files in %s \nESD tree name is %s \nReading %d events",
2750 dirName.Data(), esdTreeName.Data(), eventsToRead) ) ;
2752 // create a TChain of all the files
2753 TChain * cESDTree = new TChain(esdTreeName) ;
2755 // read from the directory file until the require number of events are collected
2756 void * from = gSystem->OpenDirectory(dirName) ;
2758 AliError( Form("Directory %s does not exist") ) ;
2761 else{ // reading file names from directory
2762 const char * subdir ;
2763 // search all subdirectories witch matching pattern
2764 while( (subdir = gSystem->GetDirEntry(from)) &&
2765 (cESDTree->GetEntries() < eventsToRead)) {
2766 if ( strstr(subdir, pattern) != 0 ) {
2768 sprintf(file, "%s%s/AliESDs.root", dirName.Data(), subdir);
2769 AliInfo( Form("Adding %s\n", file) );
2770 cESDTree->Add(file) ;
2774 AliInfo( Form(" %d events read", cESDTree->GetEntriesFast()) ) ;
2777 } // reading file names from directory
2781 //__________________________________________________________________________
2782 TChain * AliPHOSGammaJet::ReadESD(const UInt_t eventsToRead,
2783 const TString dirName,
2784 const TString esdTreeName,
2785 const char * pattern)
2787 // Read AliESDs files and return a Chain of events
2789 if ( dirName == "" ) {
2790 AliError("Give the name of the DIR where to find files") ;
2793 if ( esdTreeName == "" )
2794 return ReadESDfromdisk(eventsToRead, dirName) ;
2795 else if ( strcmp(pattern, "") == 0 )
2796 return ReadESDfromdisk(eventsToRead, dirName, esdTreeName) ;
2798 return ReadESDfromdisk(eventsToRead, dirName, esdTreeName, pattern) ;
2801 //___________________________________________________________________
2802 void AliPHOSGammaJet::SetJet(TParticle * part, Bool_t & b, Float_t cone,
2803 Double_t eta, Double_t phi)
2806 //Check if the particle is inside the cone defined by the leading particle
2809 if(phi > TMath::TwoPi())
2810 phi-=TMath::TwoPi();
2812 phi+=TMath::TwoPi();
2814 Double_t rad = 10000 + cone;
2816 if(TMath::Abs(part->Phi()-phi) <= (TMath::TwoPi() - cone))
2817 rad = TMath::Sqrt(TMath::Power(part->Eta()-eta,2)+
2818 TMath::Power(part->Phi()-phi,2));
2820 if(part->Phi()-phi > TMath::TwoPi() - cone)
2821 rad = TMath::Sqrt(TMath::Power(part->Eta()-eta,2)+
2822 TMath::Power((part->Phi()-TMath::TwoPi())-phi,2));
2823 if(part->Phi()-phi < -(TMath::TwoPi() - cone))
2824 rad = TMath::Sqrt(TMath::Power(part->Eta()-eta,2)+
2825 TMath::Power((part->Phi()+TMath::TwoPi())-phi,2));
2833 //____________________________________________________________________________
2834 Double_t AliPHOSGammaJet::SigmaE(Double_t energy)
2836 // Calculates the energy dependent energy resolution
2840 rv = TMath::Sqrt( TMath::Power(fResPara1/energy, 2)
2841 + TMath::Power(fResPara2/TMath::Sqrt(energy), 2)
2842 + TMath::Power(fResPara3, 2) ) ;
2844 return rv * energy ;
2847 //____________________________________________________________________________
2848 Double_t AliPHOSGammaJet::SigmaP(Double_t energy)
2850 // Calculates the energy dependent position resolution
2852 Double_t sigma = TMath::Sqrt(TMath::Power(fPosParaA,2) +
2853 TMath::Power(fPosParaB,2) / energy) ;
2856 return sigma ; // in cm