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.11 2006/01/31 20:30:52 hristov
23 * Revision 1.10 2006/01/23 18:04:08 hristov
24 * Removing meaningless const
26 * Revision 1.9 2006/01/12 16:23:26 schutz
27 * ESD is properly read with methods of macros/AliReadESD.C copied in it
29 * Revision 1.8 2005/12/20 07:08:32 schutz
30 * corrected error in call AliReadESD
32 * Revision 1.6 2005/05/28 14:19:04 schutz
33 * Compilation warnings fixed by T.P.
37 //_________________________________________________________________________
38 // Class for the analysis of gamma-jet correlations
39 // Basically it seaches for a prompt photon in the PHOS acceptance,
40 // if so we construct a jet around the highest pt particle in the opposite
41 // side in azimuth, inside the TPC and EMCAL acceptances. First the leading
42 // particle and then the jet have to fullfill several conditions
43 // (energy, direction ..) to be accepted. Then the fragmentation function
44 // of this jet is constructed
46 //*-- Author: Gustavo Conesa & Yves Schutz (IFIC, CERN)
47 //////////////////////////////////////////////////////////////////////////////
50 // --- ROOT system ---
53 #include <TParticle.h>
55 #include <TPaveLabel.h>
59 #include "AliPHOSGammaJet.h"
60 #include "AliPHOSGetter.h"
61 #include "AliPHOSGeometry.h"
62 #include "AliPHOSFastGlobalReconstruction.h"
64 #include "AliESDtrack.h"
65 #include "AliESDCaloCluster.h"
66 #include "Riostream.h"
69 ClassImp(AliPHOSGammaJet)
71 //____________________________________________________________________________
72 AliPHOSGammaJet::AliPHOSGammaJet() : TTask()
75 fAngleMaxParam.Set(4) ;
76 fAngleMaxParam.Reset(0.);
83 fJetRatioMaxCut = 0. ;
84 fJetRatioMinCut = 0. ;
85 fJetTPCRatioMaxCut = 0. ;
86 fJetTPCRatioMinCut = 0. ;
102 fTPCCutsLikeEMCAL = 0 ;
109 for(Int_t i = 0; i<10; i++){
113 fPtJetSelectionCut = 0.0 ;
114 fNamePtThres[i] = "" ;
125 fJetSigma1[i] = 0.0 ;
126 fJetSigma2[i] = 0.0 ;
127 fPhiEMCALCut[i] = 0.0 ;
133 fOutputFile = new TFile(gDirectory->GetName()) ;
134 fInputFileName = gDirectory->GetName() ;
135 fOutputFileName = gDirectory->GetName() ;
136 fHIJINGFileName = gDirectory->GetName() ;
145 fListHistos = new TObjArray(100) ;
146 TList * list = gDirectory->GetListOfKeys() ;
150 for (index = 0 ; index < list->GetSize()-1 ; index++) {
151 //-1 to avoid GammaJet Task
152 h = dynamic_cast<TH2F*>(gDirectory->Get(list->At(index)->GetName())) ;
153 fListHistos->Add(h) ;
158 //____________________________________________________________________________
159 AliPHOSGammaJet::AliPHOSGammaJet(const TString inputfilename) :
160 TTask("GammaJet","Analysis of gamma-jet correlations")
163 fInputFileName = inputfilename;
164 fFastRec = new AliPHOSFastGlobalReconstruction(fInputFileName);
165 AliPHOSGetter * gime = AliPHOSGetter::Instance(fInputFileName) ;
166 fNEvent = gime->MaxEvent();
171 //____________________________________________________________________________
172 AliPHOSGammaJet::AliPHOSGammaJet(const AliPHOSGammaJet & gj) : TTask(gj)
175 fAngleMaxParam = gj.fAngleMaxParam;
176 fAnyConeOrPt = gj.fAnyConeOrPt;
177 fESDdata = gj.fESDdata;
178 fEtaCut = gj.fEtaCut ;
179 fInvMassMaxCut = gj.fInvMassMaxCut ;
180 fInvMassMinCut = gj.fInvMassMinCut ;
181 fFastRec = gj.fFastRec ;
182 fOptionGJ = gj.fOptionGJ ;
183 fMinDistance = gj.fMinDistance ;
184 fOptFast = gj.fOptFast ;
185 fOnlyCharged = gj.fOnlyCharged ;
186 fOutputFile = gj.fOutputFile ;
187 fInputFileName = gj.fInputFileName ;
188 fOutputFileName = gj.fOutputFileName ;
189 fHIJINGFileName = gj.fHIJINGFileName ;
190 fHIJING = gj.fHIJING ;
191 fRatioMaxCut = gj.fRatioMaxCut ;
192 fRatioMinCut = gj.fRatioMinCut ;
193 fJetRatioMaxCut = gj.fJetRatioMaxCut ;
194 fJetRatioMinCut = gj.fJetRatioMinCut ;
195 fJetTPCRatioMaxCut = gj.fJetRatioMaxCut ;
196 fJetTPCRatioMinCut = gj.fJetRatioMinCut ;
197 fNEvent = gj.fNEvent ;
200 fResPara1 = gj.fResPara1 ;
201 fResPara2 = gj.fResPara2 ;
202 fResPara3 = gj.fResPara3 ;
204 fNeutralPtCut = gj.fNeutralPtCut ;
205 fChargedPtCut = gj.fChargedPtCut ;
206 fPtJetSelectionCut = gj.fPtJetSelectionCut ;
207 fPhiMaxCut = gj.fPhiMaxCut ;
208 fPhiMinCut = gj.fPhiMinCut ;
209 fPosParaA = gj.fPosParaA ;
210 fPosParaB = gj.fPosParaB ;
211 fSelect = gj.fSelect ;
212 fTPCCutsLikeEMCAL = gj.fTPCCutsLikeEMCAL ;
214 fPtThreshold = gj.fPtThreshold ;
216 fDirName = gj.fDirName ;
217 fESDTree = gj.fESDTree ;
218 fPattern = gj.fPattern ;
220 SetName (gj.GetName()) ;
221 SetTitle(gj.GetTitle()) ;
223 for(Int_t i = 0; i<10; i++){
224 fCones[i] = gj.fCones[i] ;
225 fNameCones[i] = gj.fNameCones[i] ;
226 fPtThres[i] = gj.fPtThres[i] ;
227 fNamePtThres[i] = gj.fNamePtThres[i] ;
229 fJetXMin1[i] = gj.fJetXMin1[i] ;
230 fJetXMin2[i] = gj.fJetXMin2[i] ;
231 fJetXMax1[i] = gj.fJetXMax1[i] ;
232 fJetXMax2[i] = gj.fJetXMax2[i] ;
233 fBkgMean[i] = gj.fBkgMean[i] ;
234 fBkgRMS[i] = gj.fBkgRMS[i] ;
236 fJetE1[i] = gj.fJetE1[i] ;
237 fJetE2[i] = gj.fJetE2[i] ;
238 fJetSigma1[i] = gj.fJetSigma1[i] ;
239 fJetSigma2[i] = gj.fJetSigma2[i] ;
240 fPhiEMCALCut[i] = gj.fPhiEMCALCut[i] ;
246 //____________________________________________________________________________
247 AliPHOSGammaJet::~AliPHOSGammaJet()
249 fOutputFile->Close() ;
253 //____________________________________________________________________________
254 void AliPHOSGammaJet::AddHIJINGToList(Int_t iEvent, TClonesArray * particleList,
257 TClonesArray * plNePHOS)
260 // List of particles copied to a root file.
262 // Char_t sb[] = "bgrd/";
263 // // cout<<sb<<endl;
265 // Char_t sf[] = "/list.root";
266 // // cout<<sf<<endl;
267 // sprintf(si,"%d",iEvent);
269 // // cout<<si<<endl;
271 // // cout<<si<<endl;
272 // TFile * f = TFile::Open(sb) ;
273 // //cout<<f->GetName()<<endl;
276 sprintf(fi,"bgrd/%d/list.root",iEvent);
277 TFile * f = TFile::Open(fi) ;
278 //cout<<f->GetName()<<endl;
280 TParticle *particle = new TParticle();
281 TTree *t = (TTree*) f->Get("T");
282 TBranch *branch = t->GetBranch("primaries");
283 branch->SetAddress(&particle);
285 Int_t index = particleList->GetEntries() ;
286 Int_t indexNe = plNe->GetEntries() ;
287 Int_t indexCh = plCh->GetEntries() ;
288 Int_t indexNePHOS = plNePHOS->GetEntries() ;
289 Double_t charge = 0.;
290 Int_t iParticle = 0 ;
292 Double_t x = 0., z = 0.;
293 // cout<<"bkg entries "<<t->GetEntries()<<endl;
295 AliPHOSGetter * gime = AliPHOSGetter::Instance() ;
296 const AliPHOSGeometry * geom = gime->PHOSGeometry() ;
299 for (iParticle=0 ; iParticle < t->GetEntries() ; iParticle++) {
300 t->GetEvent(iParticle) ;
306 charge = TDatabasePDG::Instance()
307 ->GetParticle(particle->GetPdgCode())->Charge();
309 if((charge != 0) && (particle->Pt() > fChargedPtCut)){
310 if(TMath::Abs(particle->Eta())<fEtaCut){
311 new((*particleList)[index]) TParticle(*particle) ;
312 (dynamic_cast<TParticle*>(particleList->At(index)))
316 new((*plCh)[indexCh]) TParticle(*particle) ;
317 (dynamic_cast<TParticle*>(plCh->At(indexCh)))->SetStatusCode(0) ;
320 if(strstr(fOptionGJ,"deb all")||strstr(fOptionGJ,"deb")){
321 dynamic_cast<TH1F*>(fListHistos->FindObject("PtSpectra"))->
322 Fill(particle->Pt());
325 else if((charge == 0) && (particle->Pt() > fNeutralPtCut) ){
326 geom->ImpactOnEmc(particle->Theta(),particle->Phi(), m,z,x);
330 if(strstr(fOptionGJ,"deb all")|| strstr(fOptionGJ,"deb"))
331 dynamic_cast<TH1F*>(fListHistos->FindObject("PtSpectra"))->
332 Fill(particle->Pt());
334 new((*plNePHOS)[indexNePHOS]) TParticle(*particle) ;
335 (dynamic_cast<TParticle*>(plNePHOS->At(indexNePHOS)))->SetStatusCode(0) ;
339 if((particle->Phi()>fPhiEMCALCut[0]) &&
340 (particle->Phi()<fPhiEMCALCut[1]) && m == 0)
342 if(strstr(fOptionGJ,"deb all")|| strstr(fOptionGJ,"deb"))
343 dynamic_cast<TH1F*>(fListHistos->FindObject("PtSpectra"))->
344 Fill(particle->Pt());
346 new((*particleList)[index]) TParticle(*particle) ;
347 (dynamic_cast<TParticle*>(particleList->At(index)))
350 new((*plNe)[indexNe]) TParticle(*particle) ;
351 (dynamic_cast<TParticle*>(plNe->At(indexNe)))->SetStatusCode(0) ;
360 Double_t mass = TDatabasePDG::Instance()->GetParticle(111)->Mass();
361 TLorentzVector pPi0, pGamma1, pGamma2 ;
362 Double_t angle = 0, cellDistance = 0.;
365 // fFastRec = new AliPHOSFastGlobalReconstruction(fHIJINGFileName);
366 // fFastRec->FastReconstruction(iiEvent);
368 for (iParticle=0 ; iParticle < t->GetEntries() ; iParticle++) {
369 t->GetEvent(iParticle) ;
373 charge = TDatabasePDG::Instance()
374 ->GetParticle(particle->GetPdgCode())->Charge();
376 if((charge != 0) && (particle->Pt() > fChargedPtCut)
377 && (TMath::Abs(particle->Eta())<fEtaCut)){
379 new((*particleList)[index]) TParticle(*particle) ;
380 (dynamic_cast<TParticle*>(particleList->At(index)))
384 new((*plCh)[indexCh]) TParticle(*particle) ;
385 (dynamic_cast<TParticle*>(plCh->At(indexCh)))->SetStatusCode(0) ;
388 else if(charge == 0){
389 geom->ImpactOnEmc(particle->Theta(),particle->Phi(), m,z,x);
390 if((particle->GetPdgCode() != 111) && (particle->Pt() > fNeutralPtCut) &&
391 (TMath::Abs(particle->Eta())<fEtaCut) ){
392 TLorentzVector part(particle->Px(),particle->Py(),
393 particle->Pz(),particle->Energy());
395 if(part.Pt() > fNeutralPtCut){
397 if(particle->Phi()>fPhiEMCALCut[0] &&
398 particle->Phi()<fPhiEMCALCut[1] && m == 0)
400 new((*particleList)[index]) TParticle(*particle) ;
401 (dynamic_cast<TParticle*>(particleList->At(index)))
403 (dynamic_cast<TParticle*>(particleList->At(index)))
404 ->SetMomentum(part.Px(),part.Py(),part.Pz(),part.E());
408 new((*plNe)[indexNe]) TParticle(*particle) ;
409 (dynamic_cast<TParticle*>(plNe->At(indexNe)))
410 ->SetMomentum(part.Px(),part.Py(),part.Pz(),part.E());
411 (dynamic_cast<TParticle*>(plNe->At(indexNe)))->SetStatusCode(0) ;
416 new((*plNePHOS)[indexNePHOS]) TParticle(*particle) ;
417 (dynamic_cast<TParticle*>(plNePHOS->At(indexNePHOS)))
418 ->SetMomentum(part.Px(),part.Py(),part.Pz(),part.E());
419 (dynamic_cast<TParticle*>(plNePHOS->At(indexNePHOS)))->SetStatusCode(0) ;
424 if((particle->GetPdgCode() == 111) && (particle->Pt() > fNeutralPtCut) &&
425 (TMath::Abs(particle->Eta())<fEtaCut+1))
428 pPi0.SetPxPyPzE(particle->Px(),particle->Py(),particle->Pz(),
433 Pi0Decay(mass,pPi0,pGamma1,pGamma2,angle);
434 //Check if decay photons are too close for PHOS
435 cellDistance = angle*460; //cm
436 if (cellDistance < fMinDistance) {
437 if(strstr(fOptionGJ,"deb all")|| strstr(fOptionGJ,"deb"))
438 dynamic_cast<TH1F*>(fListHistos->FindObject("PtSpectra"))->
439 Fill(particle->Pt());
441 //Pi0 inside phi EMCAL acceptance
443 TLorentzVector part(particle->Px(),particle->Py(),
444 particle->Pz(),particle->Energy());
446 if(part.Pt() > fNeutralPtCut){
447 if(particle->Phi()>fPhiEMCALCut[0] &&
448 particle->Phi()<fPhiEMCALCut[1] && m == 0){
450 new((*particleList)[index]) TParticle(*particle) ;
451 (dynamic_cast<TParticle*>(particleList->At(index)))->SetStatusCode(0) ;
452 (dynamic_cast<TParticle*>(particleList->At(index)))
453 ->SetMomentum(part.Px(),part.Py(),part.Pz(),part.E());
456 new((*plNe)[indexNe]) TParticle(*particle) ;
457 (dynamic_cast<TParticle*>(plNe->At(indexNe))) ->SetStatusCode(0) ;
458 (dynamic_cast<TParticle*>(plNe->At(indexNe)))
459 ->SetMomentum(part.Px(),part.Py(),part.Pz(),part.E());
463 if(strstr(fOptionGJ,"deb all")|| strstr(fOptionGJ,"deb"))
464 dynamic_cast<TH1F*>(fListHistos->FindObject("PtSpectra"))->
465 Fill(particle->Pt());
466 new((*plNePHOS)[indexNePHOS]) TParticle(*particle) ;
467 (dynamic_cast<TParticle*>(plNePHOS->At(indexNePHOS))) ->SetStatusCode(0) ;
468 (dynamic_cast<TParticle*>(plNePHOS->At(indexNePHOS)))
469 ->SetMomentum(part.Px(),part.Py(),part.Pz(),part.E());
476 dynamic_cast<TH2F*>(fListHistos->FindObject("AnglePair"))
477 ->Fill(pPi0.E(),angle);
480 if(pGamma1.Pt() > 0. && TMath::Abs(pGamma1.Eta())<fEtaCut){
484 if(pGamma1.Pt() > fNeutralPtCut ){
486 TParticle * photon1 =
487 new TParticle(22,1,0,0,0,0,pGamma1.Px(),pGamma1.Py(),
488 pGamma1.Pz(),pGamma1.E(),0,0,0,0);
489 geom->ImpactOnEmc(photon1->Theta(),photon1->Phi(), m,z,x);
490 if( photon1->Phi()>fPhiEMCALCut[0] && photon1->Phi()<fPhiEMCALCut[1]
492 if(strstr(fOptionGJ,"deb all") || strstr(fOptionGJ,"deb"))
493 dynamic_cast<TH1F*>(fListHistos->FindObject("PtSpectra"))->
495 new((*particleList)[index]) TParticle(*photon1) ;
496 (dynamic_cast<TParticle*>(particleList->At(index)))->SetStatusCode(0) ;
499 new((*plNe)[indexNe]) TParticle(*photon1) ;
500 (dynamic_cast<TParticle*>(plNe->At(indexNe))) ->SetStatusCode(0) ;
505 if(strstr(fOptionGJ,"deb all") || strstr(fOptionGJ,"deb"))
506 dynamic_cast<TH1F*>(fListHistos->FindObject("PtSpectra"))->
508 new((*plNePHOS)[indexNePHOS]) TParticle(*photon1) ;
509 (dynamic_cast<TParticle*>(plNePHOS->At(indexNePHOS)))->SetStatusCode(0) ;
515 if(pGamma2.Pt() > 0. && TMath::Abs(pGamma2.Eta())<fEtaCut){
518 if(pGamma2.Pt() > fNeutralPtCut){
520 TParticle * photon2 =
521 new TParticle(22,1,0,0,0,0,pGamma2.Px(), pGamma2.Py(),
522 pGamma2.Pz(),pGamma2.E(),0,0,0,0);
523 geom->ImpactOnEmc(photon2->Theta(),photon2->Phi(), m,z,x);
524 if(photon2->Phi()>fPhiEMCALCut[0] &&
525 photon2->Phi()<fPhiEMCALCut[1] && m == 0){
526 if(strstr(fOptionGJ,"deb all") || strstr(fOptionGJ,"deb"))
527 dynamic_cast<TH1F*>(fListHistos->FindObject("PtSpectra"))->
529 new((*particleList)[index]) TParticle(*photon2) ;
530 (dynamic_cast<TParticle*>(particleList->At(index)))->SetStatusCode(0) ;
533 new((*plNe)[indexNe]) TParticle(*photon2) ;
534 (dynamic_cast<TParticle*>(plNe->At(indexNe))) ->SetStatusCode(0) ;
538 if(strstr(fOptionGJ,"deb all") || strstr(fOptionGJ,"deb"))
539 dynamic_cast<TH1F*>(fListHistos->FindObject("PtSpectra"))->
541 new((*plNePHOS)[indexNePHOS]) TParticle(*photon2) ;
542 (dynamic_cast<TParticle*>(plNePHOS->At(indexNePHOS))) ->SetStatusCode(0) ;
546 // e = (pGamma1+pGamma2).E();
547 // if(IsAngleInWindow(angle,e))
549 (fListHistos->FindObject("AnglePairAccepted"))->
550 Fill(pPi0.E(),angle);
553 }//photon2 in acceptance
554 }//if angle > mindist
557 }//for (iParticle<nParticle)
560 //Info("AddHIJINGToList","End HIJING");
563 //____________________________________________________________________________
564 Double_t AliPHOSGammaJet::CalculateJetRatioLimit(const Double_t ptg,
568 //Info("CalculateLimit","x1 %f, x2%f",x[0],x[1]);
569 Double_t epp = par[0] + par[1] * ptg ;
570 Double_t spp = par[2] + par[3] * ptg ;
571 Double_t f = x[0] + x[1] * ptg ;
572 Double_t epb = epp + par[4] ;
573 Double_t spb = TMath::Sqrt(spp*spp+ par[5]*par[5]) ;
574 Double_t rat = (epb - spb * f) / ptg ;
575 //Info("CalculateLimit","epp %f, spp %f, f %f", epp, spp, f);
576 //Info("CalculateLimit","epb %f, spb %f, rat %f", epb, spb, rat);
580 //____________________________________________________________________________
581 void AliPHOSGammaJet::CreateParticleList(Int_t iEvent,
582 TClonesArray * particleList,
585 TClonesArray * plNePHOS)
587 //Info("CreateParticleList","Inside");
588 AliPHOSGetter * gime = AliPHOSGetter::Instance(fInputFileName) ;
589 const AliPHOSGeometry * geom = gime->PHOSGeometry() ;
590 gime->Event(iEvent, "X") ;
593 Int_t index = particleList->GetEntries() ;
594 Int_t indexCh = plCh->GetEntries() ;
595 Int_t indexNe = plNe->GetEntries() ;
596 Int_t indexNePHOS = plNePHOS->GetEntries() ;
597 Int_t iParticle = 0 ;
598 Double_t charge = 0.;
600 Double_t x = 0., z = 0.;
603 for (iParticle=0 ; iParticle < gime->NPrimaries() ; iParticle++) {
604 const TParticle * particle = gime->Primary(iParticle) ;
610 //Keep Stable particles within eta range
611 if((particle->GetStatusCode() == 1) &&
612 (particle->Pt() > 0)){
613 if(TMath::Abs(particle->Eta())<fEtaCut){
616 charge = TDatabasePDG::Instance()
617 ->GetParticle(particle->GetPdgCode())->Charge();
618 if((charge != 0) && (particle->Pt() > fChargedPtCut)){
620 if(strstr(fOptionGJ,"deb all")|| strstr(fOptionGJ,"deb"))
621 dynamic_cast<TH1F*>(fListHistos->FindObject("PtSpectra"))->
622 Fill(particle->Pt());
623 new((*plCh)[indexCh++]) TParticle(*particle) ;
624 new((*particleList)[index++]) TParticle(*particle) ;
626 else if((charge == 0) && (particle->Pt() > fNeutralPtCut)){
627 geom->ImpactOnEmc(particle->Theta(),particle->Phi(), m,z,x);
630 if(strstr(fOptionGJ,"deb all")|| strstr(fOptionGJ,"deb"))
631 dynamic_cast<TH1F*>(fListHistos->FindObject("PtSpectra"))->
632 Fill(particle->Pt());
634 new((*plNePHOS)[indexNePHOS++]) TParticle(*particle) ;
636 if((particle->Phi()>fPhiEMCALCut[0]) &&
637 (particle->Phi()<fPhiEMCALCut[1]) && m == 0)
639 if(strstr(fOptionGJ,"deb all")|| strstr(fOptionGJ,"deb"))
640 dynamic_cast<TH1F*>(fListHistos->FindObject("PtSpectra"))->
641 Fill(particle->Pt());
642 new((*plNe)[indexNe++]) TParticle(*particle) ;
643 new((*particleList)[index++]) TParticle(*particle) ;
647 }//final particle, etacut
648 }//for (iParticle<nParticle)
652 Double_t mass = TDatabasePDG::Instance()->GetParticle(111)->Mass();
653 TLorentzVector pPi0, pGamma1, pGamma2 ;
654 Double_t angle = 0, cellDistance = 0.;
656 fFastRec = new AliPHOSFastGlobalReconstruction(fInputFileName);
657 fFastRec->FastReconstruction(iEvent);
659 for (iParticle=0 ; iParticle < gime->NPrimaries() ; iParticle++) {
660 const TParticle * particle = gime->Primary(iParticle) ;
664 //Keep Stable particles within eta range
665 if((particle->GetStatusCode() == 1) && (particle->Pt() > 0)){
669 charge = TDatabasePDG::Instance()
670 ->GetParticle(particle->GetPdgCode())->Charge();
671 if((charge != 0) && (particle->Pt() > fChargedPtCut) && (TMath::Abs(particle->Eta())<fEtaCut)){
672 new((*plCh)[indexCh++]) TParticle(*particle) ;
673 new((*particleList)[index++]) TParticle(*particle) ;
675 else if(charge == 0) {
676 geom->ImpactOnEmc(particle->Theta(),particle->Phi(), m,z,x);
677 if((particle->GetPdgCode() != 111) && particle->Pt() > 0 &&
678 (TMath::Abs(particle->Eta())<fEtaCut))
681 TLorentzVector part(particle->Px(),particle->Py(),
682 particle->Pz(),particle->Energy());
686 if(part.Pt() > fNeutralPtCut){
687 if(particle->Phi()>fPhiEMCALCut[0] &&
688 particle->Phi()<fPhiEMCALCut[1] && m == 0)
690 new((*particleList)[index]) TParticle(*particle) ;
691 (dynamic_cast<TParticle*>(particleList->At(index)))
692 ->SetMomentum(part.Px(),part.Py(),part.Pz(),part.E());
695 new((*plNe)[indexNe]) TParticle(*particle) ;
696 (dynamic_cast<TParticle*>(plNe->At(indexNe)))
697 ->SetMomentum(part.Px(),part.Py(),part.Pz(),part.E());
702 new((*plNePHOS)[indexNePHOS]) TParticle(*particle) ;
703 (dynamic_cast<TParticle*>(plNePHOS->At(indexNePHOS)))
704 ->SetMomentum(part.Px(),part.Py(),part.Pz(),part.E());
709 if((particle->GetPdgCode() == 111) && (particle->Pt() > fNeutralPtCut) &&
710 (TMath::Abs(particle->Eta())<fEtaCut+1))
713 pPi0.SetPxPyPzE(particle->Px(),particle->Py(),particle->Pz(),
718 Pi0Decay(mass,pPi0,pGamma1,pGamma2,angle);
719 //Check if decay photons are too close for PHOS
720 cellDistance = angle*460; //cm
722 if (cellDistance < fMinDistance) {
724 //Pi0 inside phi EMCAL acceptance
727 TLorentzVector part(particle->Px(),particle->Py(),
728 particle->Pz(),particle->Energy());
731 if(part.Pt() > fNeutralPtCut){
732 if(particle->Phi()>fPhiEMCALCut[0] &&
733 particle->Phi()<fPhiEMCALCut[1] && m == 0){
734 if(strstr(fOptionGJ,"deb all")|| strstr(fOptionGJ,"deb"))
735 dynamic_cast<TH1F*>(fListHistos->FindObject("PtSpectra"))->
736 Fill(particle->Pt());
738 new((*plNe)[indexNe]) TParticle(*particle) ;
739 (dynamic_cast<TParticle*>(plNe->At(indexNe)))
740 ->SetMomentum(part.Px(),part.Py(),part.Pz(),part.E());
741 new((*particleList)[index]) TParticle(*particle) ;
742 (dynamic_cast<TParticle*>(particleList->At(index)))
743 ->SetMomentum(part.Px(),part.Py(),part.Pz(),part.E());
748 if(strstr(fOptionGJ,"deb all")|| strstr(fOptionGJ,"deb"))
749 dynamic_cast<TH1F*>(fListHistos->FindObject("PtSpectra"))->
750 Fill(particle->Pt());
751 new((*plNePHOS)[indexNePHOS]) TParticle(*particle) ;
752 (dynamic_cast<TParticle*>(plNePHOS->At(indexNePHOS)))
753 ->SetMomentum(part.Px(),part.Py(),part.Pz(),part.E());
760 dynamic_cast<TH2F*>(fListHistos->FindObject("AnglePair"))
761 ->Fill(pPi0.E(),angle);
765 if(pGamma1.Pt() > 0 && TMath::Abs(pGamma1.Eta())<fEtaCut){
768 if(pGamma1.Pt() > fNeutralPtCut){
769 TParticle * photon1 =
770 new TParticle(22,1,0,0,0,0,pGamma1.Px(),pGamma1.Py(),
771 pGamma1.Pz(),pGamma1.E(),0,0,0,0);
772 geom->ImpactOnEmc(photon1->Theta(),photon1->Phi(), m,z,x);
773 if( photon1->Phi()>fPhiEMCALCut[0] && photon1->Phi()<fPhiEMCALCut[1]
775 if(strstr(fOptionGJ,"deb all") || strstr(fOptionGJ,"deb"))
776 dynamic_cast<TH1F*>(fListHistos->FindObject("PtSpectra"))->
778 new((*plNe)[indexNe++]) TParticle(*photon1) ;
779 new((*particleList)[index++]) TParticle(*photon1) ;
784 if(strstr(fOptionGJ,"deb all") || strstr(fOptionGJ,"deb"))
785 dynamic_cast<TH1F*>(fListHistos->FindObject("PtSpectra"))->
787 new((*plNePHOS)[indexNePHOS++]) TParticle(*photon1) ;
792 if(pGamma2.Pt() > 0 && TMath::Abs(pGamma2.Eta())<fEtaCut){
795 if(pGamma2.Pt() > fNeutralPtCut ){
796 TParticle * photon2 =
797 new TParticle(22,1,0,0,0,0,pGamma2.Px(), pGamma2.Py(),
798 pGamma2.Pz(),pGamma2.E(),0,0,0,0);
799 geom->ImpactOnEmc(photon2->Theta(),photon2->Phi(), m,z,x);
800 if(photon2->Phi()>fPhiEMCALCut[0] &&
801 photon2->Phi()<fPhiEMCALCut[1] && m == 0){
802 if(strstr(fOptionGJ,"deb all") || strstr(fOptionGJ,"deb"))
803 dynamic_cast<TH1F*>(fListHistos->FindObject("PtSpectra"))->
805 new((*plNe)[indexNe++]) TParticle(*photon2) ;
806 new((*particleList)[index++]) TParticle(*photon2) ;
809 if(strstr(fOptionGJ,"deb all") || strstr(fOptionGJ,"deb"))
810 dynamic_cast<TH1F*>(fListHistos->FindObject("PtSpectra"))->
812 new((*plNePHOS)[indexNePHOS++]) TParticle(*photon2) ;
816 // Float_t e = (pGamma1+pGamma2).E();
817 // if(IsAngleInWindow(angle,e))
819 (fListHistos->FindObject("AnglePairAccepted"))->
820 Fill(pPi0.E(),angle);
823 }//photon2 in acceptance
824 }//if angle > mindist
827 }//final particle, etacut
828 }//for (iParticle<nParticle)
832 //____________________________________________________________________________
833 void AliPHOSGammaJet::CreateParticleListFromESD(TClonesArray * pl,
836 TClonesArray * plNePHOS,
839 //Create a list of particles from the ESD. These particles have been measured
840 //by the Central Tracking system (TPC), PHOS and EMCAL
841 //(EMCAL not available for the moment).
842 //Info("CreateParticleListFromESD","Inside");
844 Int_t index = pl->GetEntries() ;
846 Float_t *pid = new Float_t[AliPID::kSPECIESN];
848 //########### PHOS ##############
849 //Info("CreateParticleListFromESD","Fill ESD PHOS list");
850 Int_t begphos = esd->GetFirstPHOSCluster();
851 Int_t endphos = esd->GetFirstPHOSCluster() +
852 esd->GetNumberOfPHOSClusters() ;
853 Int_t indexNePHOS = plNePHOS->GetEntries() ;
854 if(strstr(fOptionGJ,"deb all"))
855 Info("CreateParticleListFromESD","PHOS: first particle %d, last particle %d",
858 for (npar = begphos; npar < endphos; npar++) {//////////////PHOS track loop
859 AliESDCaloCluster * clus = esd->GetCaloCluster(npar) ; // retrieve track from esd
861 //Create a TParticle to fill the particle list
863 Float_t en = clus->GetClusterEnergy() ;
864 Float_t *p = new Float_t();
865 clus->GetGlobalPosition(p) ;
866 TVector3 pos(p[0],p[1],p[2]) ;
867 Double_t phi = pos.Phi();
868 Double_t theta= pos.Theta();
869 Double_t px = en*TMath::Cos(phi)*TMath::Sin(theta);;
870 Double_t py = en*TMath::Sin(phi)*TMath::Sin(theta);
871 Double_t pz = en*TMath::Cos(theta);
873 TParticle * particle = new TParticle() ;
874 particle->SetMomentum(px,py,pz,en) ;
876 //Select only photons
879 //cout<<"pid "<<pid[AliPID::kPhoton]<<endl ;
880 if( pid[AliPID::kPhoton] > 0.75)
881 new((*plNePHOS)[indexNePHOS++]) TParticle(*particle) ;
884 //########### TPC #####################
885 //Info("CreateParticleListFromESD","Fill ESD TPC list");
887 Int_t endtpc = esd->GetNumberOfTracks() ;
888 Int_t indexCh = plCh->GetEntries() ;
889 if(strstr(fOptionGJ,"deb all"))
890 Info("CreateParticleListFromESD","TPC: first particle %d, last particle %d",
893 for (npar = begtpc; npar < endtpc; npar++) {////////////// track loop
894 AliESDtrack * track = esd->GetTrack(npar) ; // retrieve track from esd
896 Double_t en = track ->GetTPCsignal() ;
897 TVector3 mom = track->P3() ;
898 Double_t px = mom.Px();
899 Double_t py = mom.Py();
900 Double_t pz = mom.Pz(); //Check with TPC people if this is correct.
902 //cout<<"TPC signal "<<en<<endl;
903 //cout<<"px "<<px<<"; py "<<py<<"; pz "<<pz<<endl;
904 TParticle * particle = new TParticle() ;
905 particle->SetMomentum(px,py,pz,en) ;
907 new((*plCh)[indexCh++]) TParticle(*particle) ;
908 new((*pl)[index++]) TParticle(*particle) ;
912 //################ EMCAL ##############
913 Double_t v[3] ; //vertex ;
914 esd->GetVertex()->GetXYZ(v) ;
915 //##########Uncomment when ESD for EMCAL works ##########
916 //Info("CreateParticleListFromESD","Fill ESD EMCAL list");
918 Int_t begem = esd->GetFirstEMCALCluster();
919 Int_t endem = esd->GetFirstEMCALCluster() +
920 esd->GetNumberOfEMCALClusters() ;
921 Int_t indexNe = plNe->GetEntries() ;
922 if(strstr(fOptionGJ,"deb all"))
923 Info("CreateParticleListFromESD","EMCAL: first particle %d, last particle %d",
926 for (npar = begem; npar < endem; npar++) {//////////////EMCAL track loop
927 AliESDCaloCluster * clus = esd->GetCaloCluster(npar) ; // retrieve track from esd
929 Float_t en = clus->GetClusterEnergy() ;
930 Float_t *p = new Float_t();
931 clus->GetGlobalPosition(p) ;
932 TVector3 pos(p[0],p[1],p[2]) ;
933 Double_t phi = pos.Phi();
934 Double_t theta= pos.Theta();
935 Double_t px = en*TMath::Cos(phi)*TMath::Sin(theta);;
936 Double_t py = en*TMath::Sin(phi)*TMath::Sin(theta);
937 Double_t pz = en*TMath::Cos(theta);
938 //cout<<"EMCAL signal "<<en<<endl;
939 //cout<<"px "<<px<<"; py "<<py<<"; pz "<<pz<<endl;
940 //TParticle * particle = new TParticle() ;
941 //particle->SetMomentum(px,py,pz,en) ;
944 // //Uncomment if PID IS WORKING, photon and pi0 idenitification.
945 // //if( pid[AliPID::kPhoton] > 0.75) //This has to be fixen.
947 // //else if( pid[AliPID::kPi0] > 0.75)
949 pdg = 22; //No PID, assume all photons
950 TParticle * particle = new TParticle(pdg, 1, -1, -1, -1, -1,
951 px, py, pz, en, v[0], v[1], v[2], 0);
953 new((*plNe)[indexNe++]) TParticle(*particle) ;
954 new((*pl)[index++]) TParticle(*particle) ;
957 // Info("CreateParticleListFromESD","End Inside");
963 //____________________________________________________________________________
964 void AliPHOSGammaJet::Exec(Option_t *option)
974 // Create chain of esd trees
975 const UInt_t kNevent = static_cast<UInt_t>(GetNEvent()) ;
976 t = ReadESD(kNevent, fDirName, fESDTree, fPattern) ;
978 AliError("Could not create the TChain") ;
983 t->SetBranchAddress("ESD",&esd); // point to the container esd where to put the event from the esdTree
988 // AliGenPythia* pyth = (AliGenPythia*) gAlice->Generator();
991 TClonesArray * particleList = new TClonesArray("TParticle",1000);
992 TClonesArray * plCh = new TClonesArray("TParticle",1000);
993 TClonesArray * plNe = new TClonesArray("TParticle",1000);
994 TClonesArray * plNePHOS = new TClonesArray("TParticle",1000);
996 for (Int_t iEvent = 0 ; iEvent < fNEvent ; iEvent++) {
997 if(strstr(fOptionGJ,"deb")||strstr(fOptionGJ,"deb all"))
998 Info("Exec", "Event %d", iEvent) ;
1002 Double_t phig = 0., phil = 0., phich = 0 , phipi = 0;
1003 Double_t etag = 0., etal = 0., etach = 0., etapi = 0. ;
1004 Double_t ptg = 0., ptl = 0., ptch = 0., ptpi = 0.;
1006 TLorentzVector jet (0,0,0,0);
1007 TLorentzVector jettpc(0,0,0,0);
1011 Int_t iNbytes = t->GetEntry(iEvent); // store event in esd
1012 //cout<<"nbytes "<<iNbytes<<endl;
1013 if ( iNbytes == 0 ) {
1014 AliError("Empty TChain") ;
1017 CreateParticleListFromESD(particleList, plCh,plNe,plNePHOS, esd); //,iEvent);
1020 CreateParticleList(iEvent, particleList, plCh,plNe,plNePHOS);
1022 // TLorentzVector pyjet(0,0,0,0);
1025 // Float_t jets[4][10];
1026 // pyth->SetJetReconstructionMode(1);
1027 // pyth->LoadEvent();
1028 // pyth->GetJets(nJ, nJT, jets);
1030 // Float_t pxJ = jets[0][0];
1031 // Float_t pyJ = jets[1][0];
1032 // Float_t pzJ = jets[2][0];
1033 // Float_t eJ = jets[3][0];
1034 // pyjet.SetPxPyPzE(pxJ,pyJ,pzJ,eJ ) ;
1037 // //Info("Exec",">>>>>>>>>>Number of jets !!!! %d",nJT);
1038 // for (Int_t iJ = 1; iJ < nJT; iJ++) {
1039 // Float_t pxJ = jets[0][iJ];
1040 // Float_t pyJ = jets[1][iJ];
1041 // Float_t pzJ = jets[2][iJ];
1042 // Float_t eJ = jets[3][iJ];
1043 // pyjet.SetPxPyPzE(pxJ,pyJ,pzJ,eJ ) ;
1044 // //Info("Exec",">>>>>Pythia Jet: %d, Phi %f, Eta %f, Pt %f",
1045 // // iJ,pyjet.Phi(),pyjet.Eta(),pyjet.Pt());
1051 AddHIJINGToList(iEvent, particleList, plCh,plNe, plNePHOS);
1054 Bool_t iIsInPHOS = kFALSE ;
1055 GetGammaJet(plNePHOS, ptg, phig, etag, iIsInPHOS) ;
1059 //Info("Exec"," In PHOS") ;
1060 dynamic_cast<TH1F*>(fListHistos->FindObject("NGamma"))->Fill(ptg);
1061 dynamic_cast<TH2F*>(fListHistos->FindObject("PhiGamma"))
1063 dynamic_cast<TH2F*>(fListHistos->FindObject("EtaGamma"))
1065 if(strstr(fOptionGJ,"deb")||strstr(fOptionGJ,"deb all"))
1066 Info("Exec", "Gamma: pt %f, phi %f, eta %f", ptg,
1069 // cout<<"n charged "<<plCh->GetEntries()<<endl;
1070 // cout<<"n neutral "<<plNe->GetEntries()<<endl;
1071 // cout<<"n All "<<particleList->GetEntries()<<endl;
1073 GetLeadingCharge(plCh, ptg, phig, ptch, etach, phich) ;
1074 GetLeadingPi0 (plNe, ptg, phig, ptpi, etapi, phipi) ;
1076 // cout<<"n2 charged "<<plCh->GetEntries()<<endl;
1077 // cout<<"n2 neutral "<<plNe->GetEntries()<<endl;
1078 // cout<<"n2 All "<<particleList->GetEntries()<<endl;
1083 //Is the leading cone inside EMCAL?
1084 Bool_t insidech = kFALSE ;
1085 if((phich - fCone) > fPhiEMCALCut[0] &&
1086 (phich + fCone) < fPhiEMCALCut[1]){
1089 Bool_t insidepi = kFALSE ;
1090 if((phipi - fCone) > fPhiEMCALCut[0] &&
1091 (phipi + fCone) < fPhiEMCALCut[1]){
1095 if ((ptch > 0 || ptpi > 0)){
1096 if((ptch > ptpi) && insidech){
1100 dynamic_cast<TH2F*>(fListHistos->FindObject("ChargeRatio"))
1101 ->Fill(ptg,ptch/ptg);
1102 dynamic_cast<TH2F*>(fListHistos->FindObject("DeltaPhiCharge"))
1103 ->Fill(ptg,phig-phich);
1104 dynamic_cast<TH2F*>(fListHistos->FindObject("DeltaEtaCharge"))
1105 ->Fill(ptg,etag-etach);
1106 if(strstr(fOptionGJ,"deb"))
1107 Info("Exec"," Charged Leading") ;
1109 if((ptpi > ptch) && insidepi){
1114 dynamic_cast<TH2F*>(fListHistos->FindObject("Pi0Ratio"))
1115 ->Fill(ptg,ptpi/ptg);
1116 dynamic_cast<TH2F*>(fListHistos->FindObject("DeltaPhiPi0"))
1117 ->Fill(ptg,phig-phipi);
1118 dynamic_cast<TH2F*>(fListHistos->FindObject("DeltaEtaPi0"))
1119 ->Fill(ptg,etag-etapi);
1121 if(ptpi > 0. && strstr(fOptionGJ,"deb"))
1122 Info("Exec"," Pi0 Leading") ;
1125 if(strstr(fOptionGJ,"deb"))
1126 Info("Exec","Leading pt %f, phi %f",ptl,phil);
1127 if(insidech || insidepi){
1130 MakeJet(particleList, ptg, phig, ptl, phil, etal, "", jet);
1132 if(strstr(fOptionGJ,"deb")){
1133 // Info("Exec","Pythia Jet: Phi %f, Eta %f, Pt %f",
1134 // pyjet.Phi(),pyjet.Eta(),pyjet.Pt());
1135 Info("Exec","TPC+EMCAL Jet: Phi %f, Eta %f, Pt %f",
1136 jet.Phi(),jet.Eta(),jet.Pt());
1138 // dynamic_cast<TH2F*>(fListHistos->FindObject("DeltaPhiJet"))
1139 // ->Fill(ptg,pyjet.Phi()-jet.Phi());
1140 // dynamic_cast<TH2F*>(fListHistos->FindObject("DeltaEtaJet"))
1141 // ->Fill(ptg,pyjet.Eta()-jet.Eta());
1142 // dynamic_cast<TH2F*>(fListHistos->FindObject("DeltaPtJet"))
1143 // ->Fill(ptg,pyjet.Pt()-jet.Pt());
1146 MakeJetAnyConeOrPt(particleList, ptg, phig, ptl, phil, etal, "");
1150 if(fOnlyCharged && ptch > 0.)
1152 if(strstr(fOptionGJ,"deb"))
1153 Info("Exec","Leading TPC pt %f, phi %f",ptch,phich);
1155 dynamic_cast<TH2F*>(fListHistos->FindObject("TPCRatio"))
1156 ->Fill(ptg,ptch/ptg);
1157 dynamic_cast<TH2F*>(fListHistos->FindObject("DeltaPhiTPC"))
1158 ->Fill(ptg,phig-phich);
1159 dynamic_cast<TH2F*>(fListHistos->FindObject("DeltaEtaTPC"))
1160 ->Fill(ptg,etag-etach);
1164 MakeJet(plCh, ptg, phig, ptch, phich, etach, "TPC",jettpc);
1166 if(strstr(fOptionGJ,"deb")){
1167 // Info("Exec","Pythia Jet: Phi %f, Eta %f, Pt %f",
1168 // pyjet.Phi(),pyjet.Eta(),pyjet.Pt());
1169 Info("Exec","TPC Jet: Phi %f, Eta %f, Pt %f",
1170 jettpc.Phi(),jettpc.Eta(),jettpc.Pt());
1172 // dynamic_cast<TH2F*>(fListHistos->FindObject("DeltaPhiTPCJet"))
1173 // ->Fill(ptg,pyjet.Phi()-jettpc.Phi());
1174 // dynamic_cast<TH2F*>(fListHistos->FindObject("DeltaEtaTPCJet"))
1175 // ->Fill(ptg,pyjet.Eta()-jettpc.Eta());
1176 // dynamic_cast<TH2F*>(fListHistos->FindObject("DeltaPtTPCJet"))
1177 // ->Fill(ptg,pyjet.Pt()-jettpc.Pt());
1180 MakeJetAnyConeOrPt(plCh, ptg, phig, ptch, phich, etach, "TPC");
1186 particleList->Delete() ;
1189 plNePHOS->Delete() ;
1194 delete particleList ;
1196 fOutputFile->Write() ;
1201 //____________________________________________________________________________
1202 void AliPHOSGammaJet::FillJetHistos(TClonesArray * pl, Double_t ptg,
1203 TString conf, TString type)
1205 //Fill jet fragmentation histograms if !fAnyCone,
1206 //only for fCone and fPtThres
1207 TParticle * particle = 0 ;
1212 while ( (particle = dynamic_cast<TParticle*>(next())) ) {
1214 Double_t pt = particle->Pt();
1216 charge = TDatabasePDG::Instance()
1217 ->GetParticle(particle->GetPdgCode())->Charge();
1218 if(charge != 0){//Only jet Charged particles
1220 (fListHistos->FindObject(type+conf+"Fragment"))
1223 (fListHistos->FindObject(type+conf+"PtDist"))
1229 (fListHistos->FindObject("NBkg"+conf))->Fill(ipr);
1232 //____________________________________________________________________________
1233 void AliPHOSGammaJet::FillJetHistosAnyConeOrPt(TClonesArray * pl, Double_t ptg,
1234 TString conf, TString type,
1235 TString cone, TString ptcut)
1237 //Fill jet fragmentation histograms if fAnyCone,
1238 //for several cones and pt thresholds
1239 TParticle *particle = 0;
1244 while ( (particle = dynamic_cast<TParticle*>(next())) ) {
1246 Double_t pt = particle->Pt();
1247 charge = TDatabasePDG::Instance()
1248 ->GetParticle(particle->GetPdgCode())->Charge();
1249 if(charge != 0){//Only jet Charged particles
1251 (fListHistos->FindObject(type+conf+"FragmentCone"+cone+"Pt"+ptcut))
1254 (fListHistos->FindObject(type+conf+"PtDistCone"+cone+"Pt"+ptcut))
1261 (fListHistos->FindObject("NBkg"+conf+"Cone"+cone+"Pt"+ptcut))
1266 //____________________________________________________________________________
1267 void AliPHOSGammaJet::GetGammaJet(TClonesArray * pl, Double_t &pt,
1268 Double_t &phi, Double_t &eta, Bool_t &Is) const
1270 //Search for the prompt photon in PHOS with pt > fPtCut
1275 for(Int_t ipr = 0;ipr < pl->GetEntries() ; ipr ++ ){
1276 TParticle * particle = dynamic_cast<TParticle *>(pl->At(ipr)) ;
1278 if((particle->Pt() > fPtCut) && (particle->Pt() > pt)){
1280 pt = particle->Pt();
1281 phi = particle->Phi() ;
1282 eta = particle->Eta() ;
1288 //____________________________________________________________________________
1289 void AliPHOSGammaJet::GetLeadingCharge(TClonesArray * pl,
1290 Double_t ptg, Double_t phig,
1291 Double_t &pt, Double_t &eta, Double_t &phi) const
1293 //Search for the charged particle with highest with
1294 //Phi=Phi_gamma-Pi and pT=0.1E_gamma
1299 for(Int_t ipr = 0;ipr < pl->GetEntries() ; ipr ++ ){
1301 TParticle * particle = dynamic_cast<TParticle *>(pl->At(ipr)) ;
1303 Double_t ptl = particle->Pt();
1304 Double_t rat = ptl/ptg ;
1305 Double_t phil = particle->Phi() ;
1307 if(((phig-phil)> fPhiMinCut) && ((phig-phil)<fPhiMaxCut) &&
1308 (rat > fRatioMinCut) && (rat < fRatioMaxCut) && (ptl > pt)) {
1309 eta = particle->Eta() ;
1312 //printf("GetLeadingCharge: %f %f %f %f \n", pt, eta, phi,rat) ;
1315 //printf("GetLeadingCharge: %f %f %f \n", pt, eta, phi) ;
1320 //____________________________________________________________________________
1321 void AliPHOSGammaJet::GetLeadingPi0(TClonesArray * pl,
1322 Double_t ptg, Double_t phig,
1323 Double_t &pt, Double_t &eta, Double_t &phi)
1326 //Search for the neutral pion with highest with
1327 //Phi=Phi_gamma-Pi and pT=0.1E_gamma
1331 Double_t ptl = -100.;
1332 Double_t rat = -100.;
1333 Double_t phil = -100. ;
1336 TParticle * particle = 0;
1340 while ( (particle = (TParticle*)next()) ) {
1341 if( particle->GetPdgCode() == 111){
1342 ptl = particle->Pt();
1344 phil = particle->Phi() ;
1345 e = particle->Energy();
1347 (fListHistos->FindObject("AnglePairNoCut"))->
1350 (fListHistos->FindObject("InvMassPairNoCut"))->
1353 if(((phig-phil)> fPhiMinCut) && ((phig-phil)<fPhiMaxCut) &&
1354 (rat > fRatioMinCut) && (rat < fRatioMaxCut)) {
1357 (fListHistos->FindObject("AnglePairLeadingCut"))->
1360 (fListHistos->FindObject("InvMassPairLeadingCut"))->
1364 (fListHistos->FindObject("AnglePairAngleCut"))->
1367 (fListHistos->FindObject("InvMassPairAngleCut"))->
1371 (fListHistos->FindObject("InvMassPairAllCut"))->
1374 (fListHistos->FindObject("AnglePairAllCut"))->
1379 eta = particle->Eta() ;
1383 //printf("GetLeadingPi0: %f %f %f %f %f \n", pt, eta, phi, rat, ptg) ;
1389 (fListHistos->FindObject("InvMassPairLeading"))->
1392 (fListHistos->FindObject("AnglePairLeading"))->
1397 Int_t iPrimary = -1;
1398 TLorentzVector gammai,gammaj;
1399 Double_t angle = 0., e = 0., invmass = 0.;
1400 Double_t anglef = 0., ef = 0., invmassf = 0.;
1404 while ( (particle = (TParticle*)next()) ) {
1407 ksPdg = particle->GetPdgCode();
1408 ptl = particle->Pt();
1409 if(ksPdg == 111){ //2 gamma
1411 phil = particle->Phi() ;
1412 if((ptl> pt)&& (rat > fRatioMinCut) && (rat < fRatioMaxCut) &&
1413 ((phig-phil)>fPhiMinCut)&&((phig-phil)<fPhiMaxCut)){
1414 eta = particle->Eta() ;
1419 if(ksPdg == 22){//1 gamma
1420 particle->Momentum(gammai);
1423 while ( (particle = (TParticle*)next2()) ) {
1425 if(jPrimary>iPrimary){
1426 ksPdg = particle->GetPdgCode();
1428 particle->Momentum(gammaj);
1429 //Info("GetLeadingPi0","Egammai %f, Egammaj %f",
1430 //gammai.Pt(),gammaj.Pt());
1432 ptl = (gammai+gammaj).Pt();
1433 phil = (gammai+gammaj).Phi();
1435 phil+=TMath::TwoPi();
1437 invmass = (gammai+gammaj).M();
1438 angle = gammaj.Angle(gammai.Vect());
1439 //Info("GetLeadingPi0","Angle %f", angle);
1440 e = (gammai+gammaj).E();
1443 (fListHistos->FindObject("AnglePairNoCut"))->
1446 (fListHistos->FindObject("InvMassPairNoCut"))->
1449 if((rat > fRatioMinCut) && (rat < fRatioMaxCut) &&
1450 ((phig-phil)>fPhiMinCut)&&((phig-phil)<fPhiMaxCut)){
1453 (fListHistos->FindObject("AnglePairLeadingCut"))->
1456 (fListHistos->FindObject("InvMassPairLeadingCut"))->
1459 if(IsAngleInWindow(angle,e)){
1461 (fListHistos->FindObject("AnglePairAngleCut"))->
1464 (fListHistos->FindObject("InvMassPairAngleCut"))->
1467 //Info("GetLeadingPi0","InvMass %f", invmass);
1468 if((invmass>fInvMassMinCut) && (invmass<fInvMassMaxCut)){
1470 (fListHistos->FindObject("InvMassPairAllCut"))->
1473 (fListHistos->FindObject("AnglePairAllCut"))->
1477 eta = particle->Eta() ;
1481 invmassf = invmass ;
1485 }//(invmass>0.125) && (invmass<0.145)
1486 }//gammaj.Angle(gammai.Vect())<0.04
1488 }//iprimary<jprimary
1496 (fListHistos->FindObject("InvMassPairLeading"))->
1499 (fListHistos->FindObject("AnglePairLeading"))->
1503 // printf("GetLeadingPi0: %f %f %f \n", pt, eta, phi) ;
1507 //____________________________________________________________________________
1508 void AliPHOSGammaJet::InitParameters()
1511 //Initialize the parameters of the analysis.
1513 fAngleMaxParam.Set(4) ;
1514 fAngleMaxParam.AddAt(0.4,0);//={0.4,-0.25,0.025,-2e-4};
1515 fAngleMaxParam.AddAt(-0.25,1) ;
1516 fAngleMaxParam.AddAt(0.025,2) ;
1517 fAngleMaxParam.AddAt(-2e-4,3) ;
1518 fAnyConeOrPt = kFALSE ;
1519 fOutputFileName = "GammaJet.root" ;
1521 fHIJINGFileName = "galice.root" ;
1523 fMinDistance = 3.6 ;
1525 fInvMassMaxCut = 0.15 ;
1526 fInvMassMinCut = 0.12 ;
1527 fOnlyCharged = kFALSE ;
1530 fPhiEMCALCut[0] = 60. *TMath::Pi()/180.;
1531 fPhiEMCALCut[1] = 180.*TMath::Pi()/180.;
1535 fNeutralPtCut = 0.5 ;
1536 fChargedPtCut = 0.5 ;
1537 fTPCCutsLikeEMCAL = kFALSE ;
1538 //Jet selection parameters
1540 fRatioMaxCut = 1.0 ;
1541 fRatioMinCut = 0.1 ;
1542 fJetRatioMaxCut = 1.2 ;
1543 fJetRatioMinCut = 0.8 ;
1544 fJetTPCRatioMaxCut = 1.2 ;
1545 fJetTPCRatioMinCut = 0.3 ;
1549 fESDTree = "esdTree" ;
1552 //Cut depending on gamma energy
1554 fPtJetSelectionCut = 20.; //For Low pt jets+BKG, another limits applyed
1555 //Reconstructed jet energy dependence parameters
1556 //e_jet = a1+e_gamma b2.
1557 //Index 0-> Pt>2 GeV r = 0.3; Index 1-> Pt>0.5 GeV r = 0.3
1558 fJetE1[0] = -5.75; fJetE1[1] = -4.1;
1559 fJetE2[0] = 1.005; fJetE2[1] = 1.05;
1561 //Reconstructed sigma of jet energy dependence parameters
1562 //s_jet = a1+e_gamma b2.
1563 //Index 0-> Pt>2 GeV r = 0.3; Index 1-> Pt>0.5 GeV r = 0.3
1564 fJetSigma1[0] = 2.65; fJetSigma1[1] = 2.75;
1565 fJetSigma2[0] = 0.0018; fJetSigma2[1] = 0.033;
1567 //Background mean energy and RMS
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
1571 fBkgMean[0] = 0.; fBkgMean[1] = 8.8 ; fBkgMean[2] = 69.5;
1572 fBkgMean[3] = 0.; fBkgMean[4] = 6.4; fBkgMean[5] = 48.6;
1573 fBkgRMS[0] = 0.; fBkgRMS[1] = 7.5; fBkgRMS[2] = 22.0;
1574 fBkgRMS[3] = 0.; fBkgRMS[4] = 5.4; fBkgRMS[5] = 13.2;
1576 //Factor x of min/max = E -+ x * sigma. Obtained after selecting the
1577 //limits for monoenergetic jets.
1578 //Index 0-> No BKG; Index 1-> BKG > 2 GeV;
1579 //Index 2-> (low pt jets) BKG > 0.5 GeV;
1580 //Index > 2, same for TPC conf
1582 fJetXMin1[0] =-0.69 ; fJetXMin1[1] = 0.39 ; fJetXMin1[2] =-0.88 ;
1583 fJetXMin1[3] =-2.0 ; fJetXMin1[4] =-0.442 ; fJetXMin1[5] =-1.1 ;
1584 fJetXMin2[0] = 0.066; fJetXMin2[1] = 0.038; fJetXMin2[2] = 0.034;
1585 fJetXMin2[3] = 0.25 ; fJetXMin2[4] = 0.113; fJetXMin2[5] = 0.077 ;
1586 fJetXMax1[0] =-3.8 ; fJetXMax1[1] =-0.76 ; fJetXMax1[2] =-3.6 ;
1587 fJetXMax1[3] =-2.7 ; fJetXMax1[4] =-1.21 ; fJetXMax1[5] =-3.7 ;
1588 fJetXMax2[0] =-0.012; fJetXMax2[1] =-0.022; fJetXMax2[2] = 0.016;
1589 fJetXMax2[3] =-0.024; fJetXMax2[4] =-0.008; fJetXMax2[5] = 0.027;
1592 //Photon fast reconstruction
1593 fResPara1 = 0.0255 ; // GeV
1594 fResPara2 = 0.0272 ;
1595 fResPara3 = 0.0129 ;
1597 fPosParaA = 0.096 ; // cm
1600 //Different cones and pt thresholds to construct the jet
1606 fCones[0] = 0.3 ; fNameCones[0] = "03" ;
1607 fPtThres[0] = 0.5 ; fNamePtThres[0] = "05" ;
1611 //__________________________________________________________________________-
1612 Bool_t AliPHOSGammaJet::IsAngleInWindow(const Float_t angle,const Float_t e) {
1613 //Check if the opening angle of the candidate pairs is inside
1614 //our selection windowd
1615 Bool_t result = kFALSE;
1616 Double_t mpi0 = 0.1349766;
1617 Double_t max = fAngleMaxParam.At(0)*TMath::Exp(fAngleMaxParam.At(1)*e)
1618 +fAngleMaxParam.At(2)+fAngleMaxParam.At(3)*e;
1619 Double_t arg = (e*e-2*mpi0*mpi0)/(e*e);
1620 Double_t min = 100. ;
1622 min = TMath::ACos(arg);
1624 if((angle<max)&&(angle>=min))
1630 //__________________________________________________________________________-
1631 Bool_t AliPHOSGammaJet::IsJetSelected(const Double_t ptg, const Double_t ptj,
1632 const TString type ){
1633 //Check if the energy of the reconstructed jet is within an energy window
1641 if(type == "TPC" && !fTPCCutsLikeEMCAL){
1642 iTPC = 3 ;//If(fTPCCutsLikeEMCAL) take jet energy cuts like EMCAL
1647 //Phythia alone, jets with pt_th > 0.2, r = 0.3
1648 par[0] = fJetE1[0]; par[1] = fJetE2[0];
1649 //Energy of the jet peak
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
1653 //sigma_jet = fJetSigma1[0]+fJetSigma2[0]*e_gamma, simulation fit
1654 par[4] = fBkgMean[0 + iTPC]; par[5] = fBkgRMS[0 + iTPC];
1655 //Parameters reserved for HIJING bkg.
1656 xmax[0] = fJetXMax1[0 + iTPC]; xmax[1] = fJetXMax2[0 + iTPC];
1657 xmin[0] = fJetXMin1[0 + iTPC]; xmin[1] = fJetXMin2[0 + iTPC];
1658 //Factor that multiplies sigma to obtain the best limits,
1659 //by observation, of mono jet ratios (ptjet/ptg)
1660 //X_jet = fJetX1[0]+fJetX2[0]*e_gamma
1664 if(ptg > fPtJetSelectionCut){
1665 //Phythia +HIJING with pt_th > 2 GeV/c, r = 0.3
1666 par[0] = fJetE1[0]; par[1] = fJetE2[0];
1667 //Energy of the jet peak, same as in pp
1668 //e_jet = fJetE1[0]+fJetE2[0]*e_gamma, simulation fit
1669 par[2] = fJetSigma1[0]; par[3] = fJetSigma2[0];
1670 //Sigma of the jet peak, same as in pp
1671 //sigma_jet = fJetSigma1[0]+fJetSigma2[0]*e_gamma, simulation fit
1672 par[4] = fBkgMean[1 + iTPC]; par[5] = fBkgRMS[1 + iTPC];
1673 //Mean value and RMS of HIJING Bkg
1674 xmax[0] = fJetXMax1[1 + iTPC]; xmax[1] = fJetXMax2[1 + iTPC];
1675 xmin[0] = fJetXMin1[1 + iTPC]; xmin[1] = fJetXMin2[1 + 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
1683 //Phythia + HIJING with pt_th > 0.5 GeV/c, r = 0.3
1684 par[0] = fJetE1[1]; par[1] = fJetE2[1];
1685 //Energy of the jet peak, pt_th > 2 GeV/c, r = 0.3
1686 //e_jet = fJetE1[0]+fJetE2[0]*e_gamma, simulation fit
1687 par[2] = fJetSigma1[1]; par[3] = fJetSigma2[1];
1688 //Sigma of the jet peak, pt_th > 2 GeV/c, r = 0.3
1689 //sigma_jet = fJetSigma1[0]+fJetSigma2[0]*e_gamma, simulation fit
1690 par[4] = fBkgMean[2 + iTPC]; par[5] = fBkgRMS[2 + iTPC];
1691 //Mean value and RMS of HIJING Bkg in a 0.3 cone, pt > 2 GeV.
1692 xmax[0] = fJetXMax1[2 + iTPC]; xmax[1] = fJetXMax2[2 + iTPC];
1693 xmin[0] = fJetXMin1[2 + iTPC]; xmin[1] = fJetXMin2[2 + iTPC];
1694 //Factor that multiplies sigma to obtain the best limits,
1695 //by observation, of mono jet ratios (ptjet/ptg) mixed with HIJING Bkg,
1696 //pt_th > 2 GeV, r = 0.3
1697 //X_jet = fJetX1[0]+fJetX2[0]*e_gamma
1699 }//If low pt jet in bkg
1702 //Calculate minimum and maximum limits of the jet ratio.
1703 Double_t min = CalculateJetRatioLimit(ptg, par, xmin);
1704 Double_t max = CalculateJetRatioLimit(ptg, par, xmax);
1706 //Info("IsJetSeleted","%s : Limits min %f, max %f, ptg / ptj %f",
1707 // type.Data(),min,max,ptj/ptg);
1708 if(( min < ptj/ptg ) && ( max > ptj/ptg))
1715 //____________________________________________________________________________
1716 void AliPHOSGammaJet::List() const
1720 Info("List", "%d histograms found", fListHistos->GetEntries() ) ;
1721 TIter next(fListHistos) ;
1723 while ( (h = dynamic_cast<TH2F*>(next())) )
1724 Info("List", "%s", h->GetName()) ;
1727 //____________________________________________________________________________
1728 Double_t AliPHOSGammaJet::MakeEnergy(const Double_t energy)
1730 // Smears the energy according to the energy dependent energy resolution.
1731 // A gaussian distribution is assumed
1733 Double_t sigma = SigmaE(energy) ;
1734 return fRan.Gaus(energy, sigma) ;
1738 //____________________________________________________________________________
1739 void AliPHOSGammaJet::MakeHistos()
1741 // Create histograms to be saved in output file and
1742 // stores them in a TObjectArray
1744 fOutputFile = new TFile(fOutputFileName, "recreate") ;
1746 fListHistos = new TObjArray(10000) ;
1748 // Histos gamma pt vs leading pt
1750 TH1F * hPtSpectra = new TH1F
1751 ("PtSpectra","p_{T i} vs p_{T #gamma}",200,0,200);
1752 hPtSpectra->SetXTitle("p_{T} (GeV/c)");
1753 fListHistos->Add(hPtSpectra) ;
1755 //Histos ratio charged leading pt / gamma pt vs pt
1757 TH2F * hChargeRatio = new TH2F
1758 ("ChargeRatio","p_{T leading charge} /p_{T #gamma} vs p_{T #gamma}",
1760 hChargeRatio->SetYTitle("p_{T lead charge} /p_{T #gamma}");
1761 hChargeRatio->SetXTitle("p_{T #gamma} (GeV/c)");
1762 fListHistos->Add(hChargeRatio) ;
1764 TH2F * hTPCRatio = new TH2F
1765 ("TPCRatio","p_{T leading charge} /p_{T #gamma} vs p_{T #gamma}",
1767 hTPCRatio->SetYTitle("p_{T lead charge} /p_{T #gamma}");
1768 hTPCRatio->SetXTitle("p_{T #gamma} (GeV/c)");
1769 fListHistos->Add(hTPCRatio) ;
1772 TH2F * hPi0Ratio = new TH2F
1773 ("Pi0Ratio","p_{T leading #pi^{0}} /p_{T #gamma} vs p_{T #gamma}",
1775 hPi0Ratio->SetYTitle("p_{T lead #pi^{0}} /p_{T #gamma}");
1776 hPi0Ratio->SetXTitle("p_{T #gamma} (GeV/c)");
1777 fListHistos->Add(hPi0Ratio) ;
1779 TH2F * hPhiGamma = new TH2F
1780 ("PhiGamma","#phi_{#gamma}",200,0,120,200,0,7);
1781 hPhiGamma->SetYTitle("#phi");
1782 hPhiGamma->SetXTitle("p_{T #gamma} (GeV/c)");
1783 fListHistos->Add(hPhiGamma) ;
1785 TH2F * hEtaGamma = new TH2F
1786 ("EtaGamma","#phi_{#gamma}",200,0,120,200,-0.8,0.8);
1787 hEtaGamma->SetYTitle("#eta");
1788 hEtaGamma->SetXTitle("p_{T #gamma} (GeV/c)");
1789 fListHistos->Add(hEtaGamma) ;
1791 // //Jet reconstruction check
1792 // TH2F * hDeltaPhiJet = new TH2F
1793 // ("DeltaPhiJet","#phi_{jet} - #phi_{pyth jet} vs p_{T #gamma}",
1794 // 200,0,120,200,-1.5,1.5);
1795 // hDeltaPhiJet->SetYTitle("#Delta #phi");
1796 // hDeltaPhiJet->SetXTitle("p_{T #gamma} (GeV/c)");
1797 // fListHistos->Add(hDeltaPhiJet) ;
1799 // TH2F * hDeltaPhiTPCJet = new TH2F
1800 // ("DeltaPhiTPCJet","#phi_{jet TPC} - #phi_{pyth jet} vs p_{T #gamma}",
1801 // 200,0,120,200,-1.5,1.5);
1802 // hDeltaPhiTPCJet->SetYTitle("#Delta #phi");
1803 // hDeltaPhiTPCJet->SetXTitle("p_{T #gamma} (GeV/c)");
1804 // fListHistos->Add(hDeltaPhiTPCJet) ;
1806 // TH2F * hDeltaEtaJet = new TH2F
1807 // ("DeltaEtaJet","#phi_{jet} - #phi_{pyth jet} vs p_{T #gamma}",
1808 // 200,0,120,200,-1.5,1.5);
1809 // hDeltaEtaJet->SetYTitle("#Delta #phi");
1810 // hDeltaEtaJet->SetXTitle("p_{T #gamma} (GeV/c)");
1811 // fListHistos->Add(hDeltaEtaJet) ;
1813 // TH2F * hDeltaEtaTPCJet = new TH2F
1814 // ("DeltaEtaTPCJet","#phi_{jet TPC} - #phi_{pyth jet} vs p_{T #gamma}",
1815 // 200,0,120,200,-1.5,1.5);
1816 // hDeltaEtaTPCJet->SetYTitle("#Delta #phi");
1817 // hDeltaEtaTPCJet->SetXTitle("p_{T #gamma} (GeV/c)");
1818 // fListHistos->Add(hDeltaEtaTPCJet) ;
1820 // TH2F * hDeltaPtJet = new TH2F
1821 // ("DeltaPtJet","#phi_{jet} - #phi_{pyth jet} vs p_{T #gamma}",
1822 // 200,0,120,200,0.,100.);
1823 // hDeltaPtJet->SetYTitle("#Delta #phi");
1824 // hDeltaPtJet->SetXTitle("p_{T #gamma} (GeV/c)");
1825 // fListHistos->Add(hDeltaPtJet) ;
1827 // TH2F * hDeltaPtTPCJet = new TH2F
1828 // ("DeltaPtTPCJet","#phi_{jet TPC} - #phi_{pyth jet} vs p_{T #gamma}",
1829 // 200,0,120,200,0.,100.);
1830 // hDeltaPtTPCJet->SetYTitle("#Delta #phi");
1831 // hDeltaPtTPCJet->SetXTitle("p_{T #gamma} (GeV/c)");
1832 // fListHistos->Add(hDeltaPtTPCJet) ;
1835 TH2F * hDeltaPhiCharge = new TH2F
1836 ("DeltaPhiCharge","#phi_{#gamma} - #phi_{charge} vs p_{T #gamma}",
1837 200,0,120,200,0,6.4);
1838 hDeltaPhiCharge->SetYTitle("#Delta #phi");
1839 hDeltaPhiCharge->SetXTitle("p_{T #gamma} (GeV/c)");
1840 fListHistos->Add(hDeltaPhiCharge) ;
1842 TH2F * hDeltaPhiTPC = new TH2F
1843 ("DeltaPhiTPC","#phi_{#gamma} - #phi_{charge} vs p_{T #gamma}",
1844 200,0,120,200,0,6.4);
1845 hDeltaPhiTPC->SetYTitle("#Delta #phi");
1846 hDeltaPhiTPC->SetXTitle("p_{T #gamma} (GeV/c)");
1847 fListHistos->Add(hDeltaPhiTPC) ;
1848 TH2F * hDeltaPhiPi0 = new TH2F
1849 ("DeltaPhiPi0","#phi_{#gamma} - #phi_{ #pi^{0}} vs p_{T #gamma}",
1850 200,0,120,200,0,6.4);
1851 hDeltaPhiPi0->SetYTitle("#Delta #phi");
1852 hDeltaPhiPi0->SetXTitle("p_{T #gamma} (GeV/c)");
1853 fListHistos->Add(hDeltaPhiPi0) ;
1855 TH2F * hDeltaEtaCharge = new TH2F
1856 ("DeltaEtaCharge","#eta_{#gamma} - #eta_{charge} vs p_{T #gamma}",
1857 200,0,120,200,-2,2);
1858 hDeltaEtaCharge->SetYTitle("#Delta #eta");
1859 hDeltaEtaCharge->SetXTitle("p_{T #gamma} (GeV/c)");
1860 fListHistos->Add(hDeltaEtaCharge) ;
1862 TH2F * hDeltaEtaTPC = new TH2F
1863 ("DeltaEtaTPC","#eta_{#gamma} - #eta_{charge} vs p_{T #gamma}",
1864 200,0,120,200,-2,2);
1865 hDeltaEtaTPC->SetYTitle("#Delta #eta");
1866 hDeltaEtaTPC->SetXTitle("p_{T #gamma} (GeV/c)");
1867 fListHistos->Add(hDeltaEtaTPC) ;
1869 TH2F * hDeltaEtaPi0 = new TH2F
1870 ("DeltaEtaPi0","#eta_{#gamma} - #eta_{ #pi^{0}} vs p_{T #gamma}",
1871 200,0,120,200,-2,2);
1872 hDeltaEtaPi0->SetYTitle("#Delta #eta");
1873 hDeltaEtaPi0->SetXTitle("p_{T #gamma} (GeV/c)");
1874 fListHistos->Add(hDeltaEtaPi0) ;
1878 TH2F * hAnglePair = new TH2F
1880 "Angle between #pi^{0} #gamma pair vs p_{T #pi^{0}}",
1881 200,0,50,200,0,0.2);
1882 hAnglePair->SetYTitle("Angle (rad)");
1883 hAnglePair->SetXTitle("E_{ #pi^{0}} (GeV/c)");
1884 fListHistos->Add(hAnglePair) ;
1886 TH2F * hAnglePairAccepted = new TH2F
1887 ("AnglePairAccepted",
1888 "Angle between #pi^{0} #gamma pair vs p_{T #pi^{0}}, both #gamma in eta<0.7, inside window",
1889 200,0,50,200,0,0.2);
1890 hAnglePairAccepted->SetYTitle("Angle (rad)");
1891 hAnglePairAccepted->SetXTitle("E_{ #pi^{0}} (GeV/c)");
1892 fListHistos->Add(hAnglePairAccepted) ;
1894 TH2F * hAnglePairNoCut = new TH2F
1896 "Angle between all #gamma pair vs p_{T #pi^{0}}",200,0,50,200,0,0.2);
1897 hAnglePairNoCut->SetYTitle("Angle (rad)");
1898 hAnglePairNoCut->SetXTitle("E_{ #pi^{0}} (GeV/c)");
1899 fListHistos->Add(hAnglePairNoCut) ;
1901 TH2F * hAnglePairLeadingCut = new TH2F
1902 ("AnglePairLeadingCut",
1903 "Angle between all #gamma pair that have a good phi and pt vs p_{T #pi^{0}}",
1904 200,0,50,200,0,0.2);
1905 hAnglePairLeadingCut->SetYTitle("Angle (rad)");
1906 hAnglePairLeadingCut->SetXTitle("E_{ #pi^{0}} (GeV/c)");
1907 fListHistos->Add(hAnglePairLeadingCut) ;
1909 TH2F * hAnglePairAngleCut = new TH2F
1910 ("AnglePairAngleCut",
1911 "Angle between all #gamma pair (angle + leading cut) vs p_{T #pi^{0}}"
1912 ,200,0,50,200,0,0.2);
1913 hAnglePairAngleCut->SetYTitle("Angle (rad)");
1914 hAnglePairAngleCut->SetXTitle("E_{ #pi^{0}} (GeV/c)");
1915 fListHistos->Add(hAnglePairAngleCut) ;
1917 TH2F * hAnglePairAllCut = new TH2F
1919 "Angle between all #gamma pair (angle + inv mass cut+leading) vs p_{T #pi^{0}}"
1920 ,200,0,50,200,0,0.2);
1921 hAnglePairAllCut->SetYTitle("Angle (rad)");
1922 hAnglePairAllCut->SetXTitle("E_{ #pi^{0}} (GeV/c)");
1923 fListHistos->Add(hAnglePairAllCut) ;
1925 TH2F * hAnglePairLeading = new TH2F
1926 ("AnglePairLeading",
1927 "Angle between all #gamma pair finally selected vs p_{T #pi^{0}}",
1928 200,0,50,200,0,0.2);
1929 hAnglePairLeading->SetYTitle("Angle (rad)");
1930 hAnglePairLeading->SetXTitle("E_{ #pi^{0}} (GeV/c)");
1931 fListHistos->Add(hAnglePairLeading) ;
1934 TH2F * hInvMassPairNoCut = new TH2F
1935 ("InvMassPairNoCut","Invariant Mass of all #gamma pair vs p_{T #gamma}",
1936 120,0,120,360,0,0.5);
1937 hInvMassPairNoCut->SetYTitle("Invariant Mass (GeV/c^{2})");
1938 hInvMassPairNoCut->SetXTitle("p_{T #gamma} (GeV/c)");
1939 fListHistos->Add(hInvMassPairNoCut) ;
1941 TH2F * hInvMassPairLeadingCut = new TH2F
1942 ("InvMassPairLeadingCut",
1943 "Invariant Mass of #gamma pair (leading cuts) vs p_{T #gamma}",
1944 120,0,120,360,0,0.5);
1945 hInvMassPairLeadingCut->SetYTitle("Invariant Mass (GeV/c^{2})");
1946 hInvMassPairLeadingCut->SetXTitle("p_{T #gamma} (GeV/c)");
1947 fListHistos->Add(hInvMassPairLeadingCut) ;
1949 TH2F * hInvMassPairAngleCut = new TH2F
1950 ("InvMassPairAngleCut",
1951 "Invariant Mass of #gamma pair (angle cut) vs p_{T #gamma}",
1952 120,0,120,360,0,0.5);
1953 hInvMassPairAngleCut->SetYTitle("Invariant Mass (GeV/c^{2})");
1954 hInvMassPairAngleCut->SetXTitle("p_{T #gamma} (GeV/c)");
1955 fListHistos->Add(hInvMassPairAngleCut) ;
1958 TH2F * hInvMassPairAllCut = new TH2F
1959 ("InvMassPairAllCut",
1960 "Invariant Mass of #gamma pair (angle+invmass cut+leading) vs p_{T #gamma}",
1961 120,0,120,360,0,0.5);
1962 hInvMassPairAllCut->SetYTitle("Invariant Mass (GeV/c^{2})");
1963 hInvMassPairAllCut->SetXTitle("p_{T #gamma} (GeV/c)");
1964 fListHistos->Add(hInvMassPairAllCut) ;
1966 TH2F * hInvMassPairLeading = new TH2F
1967 ("InvMassPairLeading",
1968 "Invariant Mass of #gamma pair selected vs p_{T #gamma}",
1969 120,0,120,360,0,0.5);
1970 hInvMassPairLeading->SetYTitle("Invariant Mass (GeV/c^{2})");
1971 hInvMassPairLeading->SetXTitle("p_{T #gamma} (GeV/c)");
1972 fListHistos->Add(hInvMassPairLeading) ;
1977 TH1F * hNGamma = new TH1F("NGamma","Number of #gamma over PHOS",240,0,120);
1978 hNGamma->SetYTitle("N");
1979 hNGamma->SetXTitle("p_{T #gamma}(GeV/c)");
1980 fListHistos->Add(hNGamma) ;
1982 TH1F * hNBkg = new TH1F("NBkg","bkg multiplicity",9000,0,9000);
1983 hNBkg->SetYTitle("counts");
1984 hNBkg->SetXTitle("N");
1985 fListHistos->Add(hNBkg) ;
1987 TH2F * hNLeading = new TH2F
1988 ("NLeading","Accepted Jet Leading", 240,0,120,240,0,120);
1989 hNLeading->SetYTitle("p_{T charge} (GeV/c)");
1990 hNLeading->SetXTitle("p_{T #gamma}(GeV/c)");
1991 fListHistos->Add(hNLeading) ;
1994 TH1F * hN = new TH1F("NJet","Accepted jets",240,0,120);
1996 hN->SetXTitle("p_{T #gamma}(GeV/c)");
1997 fListHistos->Add(hN) ;
2000 //Ratios and Pt dist of reconstructed (not selected) jets
2002 TH2F * hJetRatio = new TH2F
2003 ("JetRatio","p_{T jet lead}/p_{T #gamma} vs p_{T #gamma}",
2004 240,0,120,200,0,10);
2005 hJetRatio->SetYTitle("p_{T jet lead #pi^{0}}/p_{T #gamma}");
2006 hJetRatio->SetXTitle("p_{T #gamma} (GeV/c)");
2007 fListHistos->Add(hJetRatio) ;
2010 TH2F * hJetPt = new TH2F
2011 ("JetPt", "p_{T jet lead} vs p_{T #gamma}",240,0,120,400,0,200);
2012 hJetPt->SetYTitle("p_{T jet lead #pi^{0}}/p_{T #gamma}");
2013 hJetPt->SetXTitle("p_{T #gamma} (GeV/c)");
2014 fListHistos->Add(hJetPt) ;
2019 TH2F * hBkgRatio = new TH2F
2020 ("BkgRatio","p_{T bkg lead}/p_{T #gamma} vs p_{T #gamma}",
2021 240,0,120,200,0,10);
2022 hBkgRatio->SetYTitle("p_{T bkg lead charge}/p_{T #gamma}");
2023 hBkgRatio->SetXTitle("p_{T #gamma} (GeV/c)");
2024 fListHistos->Add(hBkgRatio) ;
2027 TH2F * hBkgPt = new TH2F
2028 ("BkgPt","p_{T jet lead} vs p_{T #gamma}",240,0,120,400,0,200);
2029 hBkgPt->SetYTitle("p_{T jet lead charge}/p_{T #gamma}");
2030 hBkgPt->SetXTitle("p_{T #gamma} (GeV/c)");
2031 fListHistos->Add(hBkgPt) ;
2036 TH2F * hJetFragment =
2037 new TH2F("JetFragment","x = p_{T i charged}/p_{T #gamma}",
2038 240,0.,120.,1000,0.,1.2);
2039 hJetFragment->SetYTitle("x_{T}");
2040 hJetFragment->SetXTitle("p_{T #gamma}");
2041 fListHistos->Add(hJetFragment) ;
2043 TH2F * hBkgFragment = new TH2F
2044 ("BkgFragment","x = p_{T i charged}/p_{T #gamma}",
2045 240,0.,120.,1000,0.,1.2);
2046 hBkgFragment->SetYTitle("x_{T}");
2047 hBkgFragment->SetXTitle("p_{T #gamma}");
2048 fListHistos->Add(hBkgFragment) ;
2051 new TH2F("JetPtDist","x = p_{T i charged}",240,0.,120.,400,0.,200.);
2052 hJetPtDist->SetXTitle("p_{T #gamma} (GeV/c)");
2053 fListHistos->Add(hJetPtDist) ;
2055 TH2F * hBkgPtDist = new TH2F
2056 ("BkgPtDist","x = p_{T i charged}",240,0.,120.,400,0.,200.);
2057 hBkgPtDist->SetXTitle("p_{T #gamma} (GeV/c)");
2058 fListHistos->Add(hBkgPtDist) ;
2063 TH1F * hNBkgTPC = new TH1F
2064 ("NBkgTPC","TPC bkg multiplicity ",9000,0,9000);
2065 hNBkgTPC->SetYTitle("counts");
2066 hNBkgTPC->SetXTitle("N");
2067 fListHistos->Add(hNBkgTPC) ;
2069 TH2F * hNTPCLeading = new TH2F
2070 ("NTPCLeading","Accepted TPC jet leading",240,0,120,240,0,120);
2071 hNTPCLeading->SetYTitle("p_{T charge} (GeV/c)");
2072 hNTPCLeading->SetXTitle("p_{T #gamma}(GeV/c)");
2073 fListHistos->Add(hNTPCLeading) ;
2075 TH1F * hNTPC = new TH1F("NTPCJet","Number of TPC jets",240,0,120);
2076 hNTPC->SetYTitle("N");
2077 hNTPC->SetXTitle("p_{T #gamma}(GeV/c)");
2078 fListHistos->Add(hNTPC) ;
2080 TH2F * hJetTPCRatio = new TH2F
2081 ("JetTPCRatio", "p_{T jet lead TPC}/p_{T #gamma} vs p_{T #gamma}",
2082 240,0,120,200,0,10);
2083 hJetTPCRatio->SetYTitle("p_{T jet lead TPC}/p_{T #gamma}");
2084 hJetTPCRatio->SetXTitle("p_{T #gamma} (GeV/c)");
2085 fListHistos->Add(hJetTPCRatio) ;
2087 TH2F * hBkgTPCRatio = new TH2F
2088 ("BkgTPCRatio","p_{T bkg lead TPC}/p_{T #gamma} vs p_{T #gamma}",
2089 240,0,120,200,0,10);
2090 hBkgTPCRatio->SetYTitle("p_{T bkg lead TPC}/p_{T #gamma}");
2091 hBkgTPCRatio->SetXTitle("p_{T #gamma} (GeV/c)");
2092 fListHistos->Add(hBkgTPCRatio) ;
2094 TH2F * hJetTPCPt = new TH2F
2095 ("JetTPCPt", "p_{T jet lead TPC} vs p_{T #gamma}",240,0,120,400,0,200);
2096 hJetTPCPt->SetYTitle("p_{T jet lead TPC}/p_{T #gamma}");
2097 hJetTPCPt->SetXTitle("p_{T #gamma} (GeV/c)");
2098 fListHistos->Add(hJetTPCPt) ;
2100 TH2F * hBkgTPCPt = new TH2F
2101 ("BkgTPCPt", "p_{T bkg lead TPC} vs p_{T #gamma}",240,0,120,400,0,200);
2102 hBkgTPCPt->SetYTitle("p_{T bkg lead TPC}/p_{T #gamma}");
2103 hBkgTPCPt->SetXTitle("p_{T #gamma} (GeV/c)");
2104 fListHistos->Add(hBkgTPCPt) ;
2108 TH2F * hJetTPCFragment =
2109 new TH2F("JetTPCFragment","x = p_{T i charged}/p_{T #gamma}",
2110 240,0.,120.,1000,0.,1.2);
2111 hJetTPCFragment->SetYTitle("x_{T}");
2112 hJetTPCFragment->SetXTitle("p_{T #gamma}");
2113 fListHistos->Add(hJetTPCFragment) ;
2115 TH2F * hBkgTPCFragment = new TH2F
2116 ("BkgTPCFragment","x = p_{T i charged}/p_{T #gamma}",
2117 240,0.,120.,1000,0.,1.2);
2118 hBkgTPCFragment->SetYTitle("x_{T}");
2119 hBkgTPCFragment->SetXTitle("p_{T #gamma}");
2120 fListHistos->Add(hBkgTPCFragment) ;
2123 TH2F * hJetTPCPtDist = new TH2F("JetTPCPtDist",
2124 "x = p_{T i charged}",240,0.,120.,400,0.,200.);
2125 hJetTPCPtDist->SetXTitle("p_{T #gamma} (GeV/c)");
2126 fListHistos->Add(hJetTPCPtDist) ;
2128 TH2F * hBkgTPCPtDist = new TH2F
2129 ("BkgTPCPtDist","x = p_{T i charged}",240,0.,120.,400,0.,200.);
2130 hBkgTPCPtDist->SetXTitle("p_{T #gamma} (GeV/c)");
2131 fListHistos->Add(hBkgTPCPtDist) ;
2137 //If we want to study the jet for different cones and pt. Old version
2139 TH2F * hJetRatios[5][5];
2140 TH2F * hJetTPCRatios[5][5];
2142 TH2F * hJetPts[5][5];
2143 TH2F * hJetTPCPts[5][5];
2145 TH2F * hBkgRatios[5][5];
2146 TH2F * hBkgTPCRatios[5][5];
2148 TH2F * hBkgPts[5][5];
2149 TH2F * hBkgTPCPts[5][5];
2151 TH2F * hNLeadings[5][5];
2152 TH2F * hNTPCLeadings[5][5];
2155 TH1F * hNTPCs[5][5];
2157 TH1F * hNBkgs[5][5];
2158 TH1F * hNBkgTPCs[5][5];
2160 TH2F * hJetFragments[5][5];
2161 TH2F * hBkgFragments[5][5];
2162 TH2F * hJetPtDists[5][5];
2163 TH2F * hBkgPtDists[5][5];
2165 TH2F * hJetTPCFragments[5][5];
2166 TH2F * hBkgTPCFragments[5][5];
2167 TH2F * hJetTPCPtDists[5][5];
2168 TH2F * hBkgTPCPtDists[5][5];
2171 for(Int_t icone = 0; icone<fNCone; icone++){
2172 for(Int_t ipt = 0; ipt<fNPt;ipt++){
2177 hJetRatios[icone][ipt] = new TH2F
2178 ("JetRatioCone"+fNameCones[icone]+"Pt"+fNamePtThres[ipt],
2179 "p_{T jet lead #pi^{0}}/p_{T #gamma} vs p_{T #gamma}, cone ="
2180 +fNameCones[icone]+", pt>" +fNamePtThres[ipt]+" GeV/c",
2181 240,0,120,200,0,10);
2182 hJetRatios[icone][ipt]->
2183 SetYTitle("p_{T jet lead #pi^{0}}/p_{T #gamma}");
2184 hJetRatios[icone][ipt]->SetXTitle("p_{T #gamma} (GeV/c)");
2185 fListHistos->Add(hJetRatios[icone][ipt]) ;
2188 hJetPts[icone][ipt] = new TH2F
2189 ("JetPtCone"+fNameCones[icone]+"Pt"+fNamePtThres[ipt],
2190 "p_{T jet lead #pi^{0}}/p_{T #gamma} vs p_{T #gamma}, cone ="
2191 +fNameCones[icone]+", pt>" +fNamePtThres[ipt]+" GeV/c",
2192 240,0,120,400,0,200);
2193 hJetPts[icone][ipt]->
2194 SetYTitle("p_{T jet lead #pi^{0}}/p_{T #gamma}");
2195 hJetPts[icone][ipt]->SetXTitle("p_{T #gamma} (GeV/c)");
2196 fListHistos->Add(hJetPts[icone][ipt]) ;
2200 hBkgRatios[icone][ipt] = new TH2F
2201 ("BkgRatioCone"+fNameCones[icone]+"Pt"+fNamePtThres[ipt],
2202 "p_{T bkg lead #pi^{0}}/p_{T #gamma} vs p_{T #gamma}, cone ="
2203 +fNameCones[icone]+", pt>" +fNamePtThres[ipt]+" GeV/c",
2204 240,0,120,200,0,10);
2205 hBkgRatios[icone][ipt]->
2206 SetYTitle("p_{T bkg lead #pi^{0}}/p_{T #gamma}");
2207 hBkgRatios[icone][ipt]->SetXTitle("p_{T #gamma} (GeV/c)");
2208 fListHistos->Add(hBkgRatios[icone][ipt]) ;
2212 hBkgPts[icone][ipt] = new TH2F
2213 ("BkgPtCone"+fNameCones[icone]+"Pt"+fNamePtThres[ipt],
2214 "p_{T jet lead #pi^{0}}/p_{T #gamma} vs p_{T #gamma}, cone ="
2215 +fNameCones[icone]+", pt>" +fNamePtThres[ipt]+" GeV/c",
2216 240,0,120,400,0,200);
2217 hBkgPts[icone][ipt]->
2218 SetYTitle("p_{T jet lead #pi^{0}}/p_{T #gamma}");
2219 hBkgPts[icone][ipt]->SetXTitle("p_{T #gamma} (GeV/c)");
2220 fListHistos->Add(hBkgPts[icone][ipt]) ;
2225 hJetTPCRatios[icone][ipt] = new TH2F
2226 ("JetTPCRatioCone"+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,200,0,10);
2230 hJetTPCRatios[icone][ipt]->SetYTitle("p_{T jet lead TPC}/p_{T #gamma}");
2231 hJetTPCRatios[icone][ipt]->SetXTitle("p_{T #gamma} (GeV/c)");
2232 fListHistos->Add(hJetTPCRatios[icone][ipt]) ;
2234 hBkgTPCRatios[icone][ipt] = new TH2F
2235 ("BkgTPCRatioCone"+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,200,0,10);
2239 hBkgTPCRatios[icone][ipt]->SetYTitle("p_{T bkg lead TPC}/p_{T #gamma}");
2240 hBkgTPCRatios[icone][ipt]->SetXTitle("p_{T #gamma} (GeV/c)");
2241 fListHistos->Add(hBkgTPCRatios[icone][ipt]) ;
2243 hJetTPCPts[icone][ipt] = new TH2F
2244 ("JetTPCPtCone"+fNameCones[icone]+"Pt"+fNamePtThres[ipt],
2245 "p_{T jet lead TPC}/p_{T #gamma} vs p_{T #gamma}, cone ="
2246 +fNameCones[icone]+", pt>" +fNamePtThres[ipt]+" GeV/c",
2247 240,0,120,400,0,200);
2248 hJetTPCPts[icone][ipt]->SetYTitle("p_{T jet lead TPC}/p_{T #gamma}");
2249 hJetTPCPts[icone][ipt]->SetXTitle("p_{T #gamma} (GeV/c)");
2250 fListHistos->Add(hJetTPCPts[icone][ipt]) ;
2252 hBkgTPCPts[icone][ipt] = new TH2F
2253 ("BkgTPCPtCone"+fNameCones[icone]+"Pt"+fNamePtThres[ipt],
2254 "p_{T bkg lead TPC}/p_{T #gamma} vs p_{T #gamma}, cone ="
2255 +fNameCones[icone]+", pt>" +fNamePtThres[ipt]+" GeV/c",
2256 240,0,120,400,0,200);
2257 hBkgTPCPts[icone][ipt]->SetYTitle("p_{T bkg lead TPC}/p_{T #gamma}");
2258 hBkgTPCPts[icone][ipt]->SetXTitle("p_{T #gamma} (GeV/c)");
2259 fListHistos->Add(hBkgTPCPts[icone][ipt]) ;
2264 hNBkgs[icone][ipt] = new TH1F
2265 ("NBkgCone"+fNameCones[icone]+"Pt"+fNamePtThres[ipt],
2266 "bkg multiplicity cone ="+fNameCones[icone]+", pt>"
2267 +fNamePtThres[ipt]+" GeV/c",9000,0,9000);
2268 hNBkgs[icone][ipt]->SetYTitle("counts");
2269 hNBkgs[icone][ipt]->SetXTitle("N");
2270 fListHistos->Add(hNBkgs[icone][ipt]) ;
2272 hNBkgTPCs[icone][ipt] = new TH1F
2273 ("NBkgTPCCone"+fNameCones[icone]+"Pt"+fNamePtThres[ipt],
2274 "bkg multiplicity cone ="+fNameCones[icone]+", pt>"
2275 +fNamePtThres[ipt]+" GeV/c",9000,0,9000);
2276 hNBkgTPCs[icone][ipt]->SetYTitle("counts");
2277 hNBkgTPCs[icone][ipt]->SetXTitle("N");
2278 fListHistos->Add(hNBkgTPCs[icone][ipt]) ;
2281 hNLeadings[icone][ipt] = new TH2F
2282 ("NLeadingCone"+fNameCones[icone]+"Pt"+fNamePtThres[ipt],
2283 "p_{T #gamma} vs p_{T #pi^{0}} cone ="+fNameCones[icone]+", pt>"
2284 +fNamePtThres[ipt]+" GeV/c",120,0,120,120,0,120);
2285 hNLeadings[icone][ipt]->SetYTitle("p_{T #pi^{0}}(GeV/c)");
2286 hNLeadings[icone][ipt]->SetXTitle("p_{T #gamma}(GeV/c)");
2287 fListHistos->Add(hNLeadings[icone][ipt]) ;
2289 hNs[icone][ipt] = new TH1F
2290 ("NJetCone"+fNameCones[icone]+"Pt"+fNamePtThres[ipt],
2291 "Number of neutral jets, cone ="+fNameCones[icone]+", pt>"
2292 +fNamePtThres[ipt]+" GeV/c",120,0,120);
2293 hNs[icone][ipt]->SetYTitle("N");
2294 hNs[icone][ipt]->SetXTitle("p_{T #gamma}(GeV/c)");
2295 fListHistos->Add(hNs[icone][ipt]) ;
2297 //Fragmentation Function
2298 hJetFragments[icone][ipt] = new TH2F
2299 ("JetFragmentCone"+fNameCones[icone]+"Pt"+fNamePtThres[ipt],
2300 "x_{T} = p_{T i}/p_{T #gamma}, cone ="+fNameCones[icone]+", pt>"
2301 +fNamePtThres[ipt]+" GeV/c",120,0.,120.,240,0.,1.2);
2302 hJetFragments[icone][ipt]->SetYTitle("x_{T}");
2303 hJetFragments[icone][ipt]->SetXTitle("p_{T #gamma}");
2304 fListHistos->Add(hJetFragments[icone][ipt]) ;
2306 hBkgFragments[icone][ipt] = new TH2F
2307 ("BkgFragmentCone"+fNameCones[icone]+"Pt"+fNamePtThres[ipt],
2308 "x_{T} = p_{T i}/p_{T #gamma}, cone ="+fNameCones[icone]+", pt>"
2309 +fNamePtThres[ipt]+" GeV/c",120,0.,120.,240,0.,1.2);
2310 hBkgFragments[icone][ipt]->SetYTitle("x_{T}");
2311 hBkgFragments[icone][ipt]->SetXTitle("p_{T #gamma}");
2312 fListHistos->Add(hBkgFragments[icone][ipt]) ;
2315 //Jet particle distribution
2317 hJetPtDists[icone][ipt] = new TH2F
2318 ("JetPtDistCone"+fNameCones[icone]+"Pt"+fNamePtThres[ipt],
2319 "p_{T i}, cone ="+fNameCones[icone]+", pt>" +fNamePtThres[ipt]+
2320 " GeV/c",120,0.,120.,120,0.,120.);
2321 hJetPtDists[icone][ipt]->SetXTitle("p_{T #gamma} (GeV/c)");
2322 fListHistos->Add(hJetPtDists[icone][ipt]) ;
2324 hBkgPtDists[icone][ipt] = new TH2F
2325 ("BkgPtDistCone"+fNameCones[icone]+"Pt"+fNamePtThres[ipt],
2326 "p_{T i}, cone ="+fNameCones[icone]+", pt>" +fNamePtThres[ipt]+
2327 " GeV/c",120,0.,120.,120,0.,120.);
2328 hBkgPtDists[icone][ipt]->SetXTitle("p_{T #gamma} (GeV/c)");
2329 fListHistos->Add(hBkgPtDists[icone][ipt]) ;
2333 hNTPCLeadings[icone][ipt] = new TH2F
2334 ("NTPCLeadingCone"+fNameCones[icone]+"Pt"+fNamePtThres[ipt],
2335 "p_{T #gamma} vs p_{T charge} cone ="+fNameCones[icone]+", pt>"
2336 +fNamePtThres[ipt]+" GeV/c",120,0,120,120,0,120);
2337 hNTPCLeadings[icone][ipt]->SetYTitle("p_{T charge} (GeV/c)");
2338 hNTPCLeadings[icone][ipt]->SetXTitle("p_{T #gamma}(GeV/c)");
2339 fListHistos->Add(hNTPCLeadings[icone][ipt]) ;
2341 hNTPCs[icone][ipt] = new TH1F
2342 ("NTPCJetCone"+fNameCones[icone]+"Pt"+fNamePtThres[ipt],
2343 "Number of charged jets, cone ="+fNameCones[icone]+", pt>"
2344 +fNamePtThres[ipt]+" GeV/c",120,0,120);
2345 hNTPCs[icone][ipt]->SetYTitle("N");
2346 hNTPCs[icone][ipt]->SetXTitle("p_{T #gamma}(GeV/c)");
2347 fListHistos->Add(hNTPCs[icone][ipt]) ;
2349 hJetTPCFragments[icone][ipt] = new TH2F
2350 ("JetTPCFragmentCone"+fNameCones[icone]+"Pt"+fNamePtThres[ipt],
2351 "x = p_{T i charged}/p_{T #gamma}, cone ="+fNameCones[icone]+", pt>"
2352 +fNamePtThres[ipt]+" GeV/c",120,0.,120.,240,0.,1.2);
2353 hJetTPCFragments[icone][ipt]->SetYTitle("x_{T}");
2354 hJetTPCFragments[icone][ipt]->SetXTitle("p_{T #gamma}");
2355 fListHistos->Add(hJetTPCFragments[icone][ipt]) ;
2357 hBkgTPCFragments[icone][ipt] = new TH2F
2358 ("BkgTPCFragmentCone"+fNameCones[icone]+"Pt"+fNamePtThres[ipt],
2359 "x = p_{T i charged}/p_{T #gamma}, cone ="+fNameCones[icone]+", pt>"
2360 +fNamePtThres[ipt]+" GeV/c",120,0.,120.,240,0.,1.2);
2361 hBkgTPCFragments[icone][ipt]->SetYTitle("x_{T}");
2362 hBkgTPCFragments[icone][ipt]->SetXTitle("p_{T #gamma}");
2363 fListHistos->Add(hBkgTPCFragments[icone][ipt]) ;
2365 hJetTPCPtDists[icone][ipt] = new TH2F
2366 ("JetTPCPtDistCone"+fNameCones[icone]+"Pt"+fNamePtThres[ipt],
2367 "x = p_{T i charged}, cone ="+fNameCones[icone]+", pt>"
2368 +fNamePtThres[ipt]+" GeV/c",120,0.,120.,120,0.,120.);
2369 hJetTPCPtDists[icone][ipt]->SetXTitle("p_{T #gamma} (GeV/c)");
2370 fListHistos->Add(hJetTPCPtDists[icone][ipt]) ;
2372 hBkgTPCPtDists[icone][ipt] = new TH2F
2373 ("BkgTPCPtDistCone"+fNameCones[icone]+"Pt"+fNamePtThres[ipt],
2374 "x = p_{T i charged}, cone ="+fNameCones[icone]+", pt>" +
2375 fNamePtThres[ipt]+" GeV/c",120,0.,120.,120,0.,120.);
2376 hBkgTPCPtDists[icone][ipt]->SetXTitle("p_{T #gamma} (GeV/c)");
2377 fListHistos->Add(hBkgTPCPtDists[icone][ipt]) ;
2381 }//If we want to study any cone or pt threshold
2385 //____________________________________________________________________________
2386 void AliPHOSGammaJet::MakeJet(TClonesArray * pl,
2387 Double_t ptg, Double_t phig,
2388 Double_t ptl, Double_t phil, Double_t etal,
2389 TString conf, TLorentzVector & jet)
2391 //Fill the jet with the particles around the leading particle with
2392 //R=fCone and pt_th = fPtThres. Calculate the energy of the jet and
2393 //check if we select it. Fill jet histograms
2394 Float_t ptcut = 0. ;
2396 if(ptg > fPtJetSelectionCut) ptcut = 2. ;
2400 TClonesArray * jetList = new TClonesArray("TParticle",1000);
2401 TClonesArray * bkgList = new TClonesArray("TParticle",1000);
2403 TLorentzVector bkg(0,0,0,0);
2404 TLorentzVector lv (0,0,0,0);
2406 Double_t ptjet = 0.0;
2407 Double_t ptbkg = 0.0;
2414 TParticle * particle = 0 ;
2416 while ( (particle = dynamic_cast<TParticle*>(next())) ) {
2422 SetJet(particle, b0, fCone, etal, phil) ;
2425 new((*jetList)[n0++]) TParticle(*particle) ;
2426 particle->Momentum(lv);
2427 if(particle->Pt() > ptcut ){
2429 ptjet+=particle->Pt();
2433 //Background around (phi_gamma-pi, eta_leading)
2434 SetJet(particle, b1, fCone,etal, phig) ;
2437 new((*bkgList)[n1++]) TParticle(*particle) ;
2438 if(particle->Pt() > ptcut ){
2440 ptbkg+=particle->Pt();
2448 if(strstr(fOptionGJ,"deb") || strstr(fOptionGJ,"deb all"))
2449 Info("MakeJet","Gamma pt %f, Jet pt %f, Bkg pt %f",ptg,ptjet,ptbkg);
2454 Double_t rat = ptjet/ptg ;
2455 Double_t ratbkg = ptbkg/ptg ;
2458 (fListHistos->FindObject("Jet"+conf+"Ratio"))->Fill(ptg,rat);
2460 (fListHistos->FindObject("Jet"+conf+"Pt")) ->Fill(ptg,ptjet);
2462 (fListHistos->FindObject("Bkg"+conf+"Ratio"))->Fill(ptg,ratbkg);
2464 (fListHistos->FindObject("Bkg"+conf+"Pt")) ->Fill(ptg,ptbkg);
2467 if(IsJetSelected(ptg,ptjet,conf) || fSelect){
2468 if(strstr(fOptionGJ,"deb") || strstr(fOptionGJ,"deb all"))
2469 Info("MakeJet","JetSelected");
2470 dynamic_cast<TH1F*>(fListHistos->FindObject("N"+conf+"Jet"))->
2472 dynamic_cast<TH2F*>(fListHistos->FindObject("N"+conf+"Leading"))
2474 FillJetHistos(jetList, ptg, conf, "Jet");
2475 FillJetHistos(bkgList, ptg, conf, "Bkg");
2482 //____________________________________________________________________________
2483 void AliPHOSGammaJet::MakeJetAnyConeOrPt(TClonesArray * pl, Double_t ptg,
2484 Double_t phig, Double_t ptl,
2485 Double_t phil, Double_t etal,
2488 //Fill the jet with the particles around the leading particle with
2489 //R=fCone(i) and pt_th = fPtThres(i). Calculate the energy of the jet and
2490 //check if we select it. Fill jet i histograms
2492 TClonesArray * jetList = new TClonesArray("TParticle",1000);
2493 TClonesArray * bkgList = new TClonesArray("TParticle",1000);
2495 Double_t ptjet = 0.0;
2496 Double_t ptbkg = 0.0;
2503 //Create as many jets as cones and pt thresholds are defined
2504 Double_t maxcut = fJetRatioMaxCut;
2505 Double_t mincut = fJetRatioMinCut;
2508 maxcut = fJetTPCRatioMaxCut;
2509 mincut = fJetTPCRatioMinCut;
2512 Double_t ratjet = 0;
2513 Double_t ratbkg = 0;
2515 for(Int_t icone = 0; icone<fNCone; icone++) {
2516 for(Int_t ipt = 0; ipt<fNPt;ipt++) {
2518 TString cone = fNameCones[icone] ;
2519 TString ptcut = fNamePtThres[ipt] ;
2522 TParticle * particle = 0 ;
2527 while ( (particle = dynamic_cast<TParticle*>(next())) ) {
2531 SetJet(particle, b0, fCones[icone],etal, phil) ;
2532 SetJet(particle, b1, fCones[icone],etal, phig) ;
2535 new((*jetList)[n0++]) TParticle(*particle) ;
2536 if(particle->Pt() > fPtThres[ipt] )
2537 ptjet+=particle->Pt();
2540 new((*bkgList)[n1++]) TParticle(*particle) ;
2541 if(particle->Pt() > fPtThres[ipt] )
2542 ptbkg+=particle->Pt();
2550 if(strstr(fOptionGJ,"deb")){
2551 Info("MakeJetAnyPt","cone %f, ptcut %f",fCones[icone],fPtThres[ipt]);
2552 Info("MakeJetAnyPt","pT: Gamma %f, Jet %f, Bkg %f",ptg,ptjet,ptbkg);
2555 ratjet = ptjet /ptg;
2559 (fListHistos->FindObject("Jet"+conf+"RatioCone"+cone+"Pt"+ptcut))
2562 (fListHistos->FindObject("Jet"+conf+"PtCone"+cone+"Pt"+ptcut))
2566 (fListHistos->FindObject("Bkg"+conf+"RatioCone"+cone+"Pt"+ptcut))
2570 (fListHistos->FindObject("Bkg"+conf+"PtCone"+cone+"Pt"+ptcut))
2575 if((ratjet < maxcut) && (ratjet > mincut)){
2578 (fListHistos->FindObject("N"+conf+"JetCone"+cone+"Pt"+ptcut))->
2581 (fListHistos->FindObject("N"+conf+"LeadingCone"+cone+"Pt"+ptcut))
2584 FillJetHistosAnyConeOrPt
2585 (jetList,ptg,conf,"Jet",fNameCones[icone],fNamePtThres[ipt]);
2586 FillJetHistosAnyConeOrPt
2587 (bkgList,ptg,conf,"Bkg",fNameCones[icone],fNamePtThres[ipt]);
2597 //____________________________________________________________________________
2598 void AliPHOSGammaJet::MakePhoton(TLorentzVector & particle)
2600 //Fast reconstruction for photons
2601 Double_t energy = particle.E() ;
2602 Double_t modenergy = MakeEnergy(energy) ;
2603 //Info("MakePhoton","Energy %f, Modif %f",energy,modenergy);
2605 // get the detected direction
2606 TVector3 pos = particle.Vect();
2608 TVector3 modpos = MakePosition(energy, pos) ;
2609 modpos *= modenergy / 460.;
2611 Float_t modtheta = modpos.Theta();
2612 Float_t modphi = modpos.Phi();
2614 // Set the modified 4-momentum of the reconstructed particle
2615 Float_t py = modenergy*TMath::Sin(modphi)*TMath::Sin(modtheta);
2616 Float_t px = modenergy*TMath::Cos(modphi)*TMath::Sin(modtheta);
2617 Float_t pz = modenergy*TMath::Cos(modtheta);
2619 particle.SetPxPyPzE(px,py,pz,modenergy);
2623 //____________________________________________________________________________
2624 TVector3 AliPHOSGammaJet::MakePosition(const Double_t energy, const TVector3 pos)
2626 // Smears the impact position according to the energy dependent position resolution
2627 // A gaussian position distribution is assumed
2631 Double_t sigma = SigmaP(energy) ;
2632 Double_t x = fRan.Gaus( pos.X(), sigma ) ;
2633 Double_t z = fRan.Gaus( pos.Z(), sigma ) ;
2634 Double_t y = pos.Y() ;
2640 // Info("MakePosition","Theta dif %f",pos.Theta()-newpos.Theta());
2641 // Info("MakePosition","Phi dif %f",pos.Phi()-newpos.Phi());
2645 //____________________________________________________________________________
2646 void AliPHOSGammaJet::Pi0Decay(Double_t mPi0, TLorentzVector &p0,
2647 TLorentzVector &p1, TLorentzVector &p2, Double_t &angle) {
2648 // Perform isotropic decay pi0 -> 2 photons
2649 // p0 is pi0 4-momentum (inut)
2650 // p1 and p2 are photon 4-momenta (output)
2651 // cout<<"Boost vector"<<endl;
2652 TVector3 b = p0.BoostVector();
2653 //cout<<"Parameters"<<endl;
2654 //Double_t mPi0 = p0.M();
2655 Double_t phi = TMath::TwoPi() * gRandom->Rndm();
2656 Double_t cosThe = 2 * gRandom->Rndm() - 1;
2657 Double_t cosPhi = TMath::Cos(phi);
2658 Double_t sinPhi = TMath::Sin(phi);
2659 Double_t sinThe = TMath::Sqrt(1-cosThe*cosThe);
2660 Double_t ePi0 = mPi0/2.;
2661 //cout<<"ePi0 "<<ePi0<<endl;
2662 //cout<<"Components"<<endl;
2663 p1.SetPx(+ePi0*cosPhi*sinThe);
2664 p1.SetPy(+ePi0*sinPhi*sinThe);
2665 p1.SetPz(+ePi0*cosThe);
2667 //cout<<"p1: "<<p1.Px()<<" "<<p1.Py()<<" "<<p1.Pz()<<" "<<p1.E()<<endl;
2668 //cout<<"p1 Mass: "<<p1.Px()*p1.Px()+p1.Py()*p1.Py()+p1.Pz()*p1.Pz()-p1.E()*p1.E()<<endl;
2669 p2.SetPx(-ePi0*cosPhi*sinThe);
2670 p2.SetPy(-ePi0*sinPhi*sinThe);
2671 p2.SetPz(-ePi0*cosThe);
2673 //cout<<"p2: "<<p2.Px()<<" "<<p2.Py()<<" "<<p2.Pz()<<" "<<p2.E()<<endl;
2674 //cout<<"p2 Mass: "<<p2.Px()*p2.Px()+p2.Py()*p2.Py()+p2.Pz()*p2.Pz()-p2.E()*p2.E()<<endl;
2675 //cout<<"Boost "<<b.X()<<" "<<b.Y()<<" "<<b.Z()<<endl;
2677 //cout<<"p1: "<<p1.Px()<<" "<<p1.Py()<<" "<<p1.Pz()<<" "<<p1.E()<<endl;
2679 //cout<<"p2: "<<p2.Px()<<" "<<p2.Py()<<" "<<p2.Pz()<<" "<<p2.E()<<endl;
2680 //cout<<"angle"<<endl;
2681 angle = p1.Angle(p2.Vect());
2682 //cout<<angle<<endl;
2685 //____________________________________________________________________________
2686 void AliPHOSGammaJet::Plot(TString what, Option_t * option) const
2688 //Plot some relevant histograms of the analysis
2689 TH2F * h = dynamic_cast<TH2F*>(fOutputFile->Get(what));
2693 else if (what == "all") {
2694 TCanvas * can = new TCanvas("GammaJet", "Gamma-Jet Study",10,40,1600,1200);
2696 can->Range(0,0,22,20);
2697 TPaveLabel *pl1 = new TPaveLabel(1,18,20,19.5,"Titre","br");
2698 pl1->SetFillColor(18);
2699 pl1->SetTextFont(32);
2700 pl1->SetTextColor(49);
2704 Float_t begx = -0.29, begy = 0., endx = 0., endy = 0.30 ;
2705 for (index = 0 ; index < fListHistos->GetEntries() ; index++) {
2706 TString name("pad") ;
2710 if (begx >= 1.0 || endx >= 1.0) {
2716 printf("%f %f %f %f \n", begx, begy, endx, endy) ;
2717 pad = new TPad(name,"This is a pad",begx,begy,endx,endy,33);
2720 fListHistos->At(index)->Draw(option) ;
2728 Info("Draw", "Histogram %s does not exist or unknown option", what.Data()) ;
2733 //____________________________________________________________________________
2734 void AliPHOSGammaJet::Print(const Option_t * opt) const
2737 //Print some relevant parameters set for the analysis
2741 Info("Print", "%s %s", GetName(), GetTitle() ) ;
2743 printf("Eta cut : %f\n", fEtaCut) ;
2744 printf("D phi max cut : %f\n", fPhiMaxCut) ;
2745 printf("D phi min cut : %f\n", fPhiMinCut) ;
2746 printf("Leading Ratio max cut : %f\n", fRatioMaxCut) ;
2747 printf("Leading Ratio min cut : %f\n", fRatioMinCut) ;
2748 printf("Jet Ratio max cut : %f\n", fJetRatioMaxCut) ;
2749 printf("Jet Ratio min cut : %f\n", fJetRatioMinCut) ;
2750 printf("Jet TPC Ratio max cut : %f\n", fJetTPCRatioMaxCut) ;
2751 printf("Jet TPC Ratio min cut : %f\n", fJetTPCRatioMinCut) ;
2752 printf("Fast recons : %d\n", fOptFast);
2753 printf("Inv Mass max cut : %f\n", fInvMassMaxCut) ;
2754 printf("Inv Mass min cut : %f\n", fInvMassMinCut) ;
2758 //__________________________________________________________________________
2759 TChain * AliPHOSGammaJet::ReadESDfromdisk(const UInt_t eventsToRead,
2760 const TString dirName,
2761 const TString esdTreeName,
2762 const char * pattern)
2764 // Reads ESDs from Disk
2767 AliInfo( Form("\nReading files in %s \nESD tree name is %s \nReading %d events",
2768 dirName.Data(), esdTreeName.Data(), eventsToRead) ) ;
2770 // create a TChain of all the files
2771 TChain * cESDTree = new TChain(esdTreeName) ;
2773 // read from the directory file until the require number of events are collected
2774 void * from = gSystem->OpenDirectory(dirName) ;
2776 AliError( Form("Directory %s does not exist") ) ;
2779 else{ // reading file names from directory
2780 const char * subdir ;
2781 // search all subdirectories witch matching pattern
2782 while( (subdir = gSystem->GetDirEntry(from)) &&
2783 (cESDTree->GetEntries() < eventsToRead)) {
2784 if ( strstr(subdir, pattern) != 0 ) {
2786 sprintf(file, "%s%s/AliESDs.root", dirName.Data(), subdir);
2787 AliInfo( Form("Adding %s\n", file) );
2788 cESDTree->Add(file) ;
2792 AliInfo( Form(" %d events read", cESDTree->GetEntriesFast()) ) ;
2795 } // reading file names from directory
2799 //__________________________________________________________________________
2800 TChain * AliPHOSGammaJet::ReadESD(const UInt_t eventsToRead,
2801 const TString dirName,
2802 const TString esdTreeName,
2803 const char * pattern)
2805 // Read AliESDs files and return a Chain of events
2807 if ( dirName == "" ) {
2808 AliError("Give the name of the DIR where to find files") ;
2811 if ( esdTreeName == "" )
2812 return ReadESDfromdisk(eventsToRead, dirName) ;
2813 else if ( strcmp(pattern, "") == 0 )
2814 return ReadESDfromdisk(eventsToRead, dirName, esdTreeName) ;
2816 return ReadESDfromdisk(eventsToRead, dirName, esdTreeName, pattern) ;
2819 //___________________________________________________________________
2820 void AliPHOSGammaJet::SetJet(TParticle * part, Bool_t & b, Float_t cone,
2821 Double_t eta, Double_t phi)
2824 //Check if the particle is inside the cone defined by the leading particle
2827 if(phi > TMath::TwoPi())
2828 phi-=TMath::TwoPi();
2830 phi+=TMath::TwoPi();
2832 Double_t rad = 10000 + cone;
2834 if(TMath::Abs(part->Phi()-phi) <= (TMath::TwoPi() - cone))
2835 rad = TMath::Sqrt(TMath::Power(part->Eta()-eta,2)+
2836 TMath::Power(part->Phi()-phi,2));
2838 if(part->Phi()-phi > TMath::TwoPi() - cone)
2839 rad = TMath::Sqrt(TMath::Power(part->Eta()-eta,2)+
2840 TMath::Power((part->Phi()-TMath::TwoPi())-phi,2));
2841 if(part->Phi()-phi < -(TMath::TwoPi() - cone))
2842 rad = TMath::Sqrt(TMath::Power(part->Eta()-eta,2)+
2843 TMath::Power((part->Phi()+TMath::TwoPi())-phi,2));
2851 //____________________________________________________________________________
2852 Double_t AliPHOSGammaJet::SigmaE(Double_t energy)
2854 // Calculates the energy dependent energy resolution
2858 rv = TMath::Sqrt( TMath::Power(fResPara1/energy, 2)
2859 + TMath::Power(fResPara2/TMath::Sqrt(energy), 2)
2860 + TMath::Power(fResPara3, 2) ) ;
2862 return rv * energy ;
2865 //____________________________________________________________________________
2866 Double_t AliPHOSGammaJet::SigmaP(Double_t energy)
2868 // Calculates the energy dependent position resolution
2870 Double_t sigma = TMath::Sqrt(TMath::Power(fPosParaA,2) +
2871 TMath::Power(fPosParaB,2) / energy) ;
2874 return sigma ; // in cm