]> git.uio.no Git - u/mrichter/AliRoot.git/blame - PHOS/AliPHOSGammaJet.cxx
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
AliTrackPoint and AliTrackPointArray moved from STEER to ESD library
[u/mrichter/AliRoot.git] / PHOS / AliPHOSGammaJet.cxx
CommitLineData
77414c90 1/**************************************************************************
2 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
3 * *
4 * Author: The ALICE Off-line Project. *
5 * Contributors are mentioned in the code where appropriate. *
6 * *
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 **************************************************************************/
15/* $Id$ */
16
702ab87e 17/* History of cvs commits:
18 *
19 * $Log$
eb0b1051 20 * Revision 1.6 2005/05/28 14:19:04 schutz
21 * Compilation warnings fixed by T.P.
22 *
702ab87e 23 */
24
77414c90 25//_________________________________________________________________________
26// Class for the analysis of gamma-jet correlations
e8daeb92 27// Basically it seaches for a prompt photon in the PHOS acceptance,
28// if so we construct a jet around the highest pt particle in the opposite
29// side in azimuth, inside the TPC and EMCAL acceptances. First the leading
30// particle and then the jet have to fullfill several conditions
31// (energy, direction ..) to be accepted. Then the fragmentation function
32// of this jet is constructed
77414c90 33//
34//*-- Author: Gustavo Conesa & Yves Schutz (IFIC, CERN)
35//////////////////////////////////////////////////////////////////////////////
36
37
38// --- ROOT system ---
39
77414c90 40#include "TParticle.h"
1305bc01 41#include "TCanvas.h"
42#include "TPaveLabel.h"
43#include "TPad.h"
77414c90 44#include "AliPHOSGammaJet.h"
45#include "AliPHOSGetter.h"
e8daeb92 46#include "TH2.h"
1305bc01 47#include "AliPHOSGeometry.h"
e8daeb92 48#include "AliPHOSFastGlobalReconstruction.h"
dffafdc2 49#include "AliESD.h"
50#include "AliESDtrack.h"
51#include "Riostream.h"
52#include "macros/AliReadESD.C"
e8daeb92 53//#include "../PYTHIA6/AliGenPythia.h"
77414c90 54
55ClassImp(AliPHOSGammaJet)
56
57//____________________________________________________________________________
1305bc01 58AliPHOSGammaJet::AliPHOSGammaJet() : TTask()
59{
77414c90 60 // ctor
1305bc01 61 fAngleMaxParam.Set(4) ;
62 fAngleMaxParam.Reset(0.);
e8daeb92 63 fAnyConeOrPt = 0 ;
1305bc01 64 fEtaCut = 0. ;
dffafdc2 65 fESDdata = 0 ;
e8daeb92 66 fFastRec = 0 ;
1305bc01 67 fInvMassMaxCut = 0. ;
68 fInvMassMinCut = 0. ;
69 fJetRatioMaxCut = 0. ;
70 fJetRatioMinCut = 0. ;
71 fJetTPCRatioMaxCut = 0. ;
72 fJetTPCRatioMinCut = 0. ;
73 fMinDistance = 0. ;
e8daeb92 74 fNEvent = 0 ;
75 fNCone = 0 ;
76 fNPt = 0 ;
77 fOnlyCharged = 0 ;
78 fOptFast = 0 ;
1305bc01 79 fPhiMaxCut = 0. ;
80 fPhiMinCut = 0. ;
81 fPtCut = 0. ;
82 fNeutralPtCut = 0. ;
83 fChargedPtCut = 0. ;
84 fRatioMaxCut = 0. ;
85 fRatioMinCut = 0. ;
86 fCone = 0 ;
87 fPtThreshold = 0 ;
e8daeb92 88 fTPCCutsLikeEMCAL = 0 ;
89 fSelect = 0 ;
dffafdc2 90
91 fDirName = "" ;
92 fESDTree = "" ;
93 fPattern = "" ;
94
1305bc01 95 for(Int_t i = 0; i<10; i++){
96 fCones[i] = 0.0 ;
97 fNameCones[i] = "" ;
98 fPtThres[i] = 0.0 ;
99 fPtJetSelectionCut = 0.0 ;
100 fNamePtThres[i] = "" ;
101 if( i < 6 ){
102 fJetXMin1[i] = 0.0 ;
103 fJetXMin2[i] = 0.0 ;
104 fJetXMax1[i] = 0.0 ;
105 fJetXMax2[i] = 0.0 ;
106 fBkgMean[i] = 0.0 ;
107 fBkgRMS[i] = 0.0 ;
108 if( i < 2 ){
109 fJetE1[i] = 0.0 ;
110 fJetE2[i] = 0.0 ;
111 fJetSigma1[i] = 0.0 ;
112 fJetSigma2[i] = 0.0 ;
113 fPhiEMCALCut[i] = 0.0 ;
114 }
115 }
116 }
117
e8daeb92 118 fOptionGJ = "" ;
1305bc01 119 fOutputFile = new TFile(gDirectory->GetName()) ;
120 fInputFileName = gDirectory->GetName() ;
121 fOutputFileName = gDirectory->GetName() ;
122 fHIJINGFileName = gDirectory->GetName() ;
123 fHIJING = 0 ;
124 fPosParaA = 0. ;
125 fPosParaB = 0. ;
126 fRan = 0 ;
127 fResPara1 = 0. ;
128 fResPara2 = 0. ;
129 fResPara3 = 0. ;
130
131 fListHistos = new TObjArray(100) ;
132 TList * list = gDirectory->GetListOfKeys() ;
133 TIter next(list) ;
134 TH2F * h = 0 ;
135 Int_t index ;
136 for (index = 0 ; index < list->GetSize()-1 ; index++) {
137 //-1 to avoid GammaJet Task
138 h = dynamic_cast<TH2F*>(gDirectory->Get(list->At(index)->GetName())) ;
139 fListHistos->Add(h) ;
140 }
141 List() ;
77414c90 142}
143
144//____________________________________________________________________________
1305bc01 145AliPHOSGammaJet::AliPHOSGammaJet(const TString inputfilename) :
146 TTask("GammaJet","Analysis of gamma-jet correlations")
147{
dffafdc2 148 // ctor
149// gSystem->Load("$ALICE_ROOT/PHOS/macros/AliReadESD_C.so");
1305bc01 150 fInputFileName = inputfilename;
151 fFastRec = new AliPHOSFastGlobalReconstruction(fInputFileName);
152 AliPHOSGetter * gime = AliPHOSGetter::Instance(fInputFileName) ;
153 fNEvent = gime->MaxEvent();
154 InitParameters();
155 fListHistos = 0 ;
77414c90 156}
157
158//____________________________________________________________________________
1305bc01 159AliPHOSGammaJet::AliPHOSGammaJet(const AliPHOSGammaJet & gj) : TTask(gj)
160{
161 // cpy ctor
162 fAngleMaxParam = gj.fAngleMaxParam;
163 fAnyConeOrPt = gj.fAnyConeOrPt;
dffafdc2 164 fESDdata = gj.fESDdata;
1305bc01 165 fEtaCut = gj.fEtaCut ;
166 fInvMassMaxCut = gj.fInvMassMaxCut ;
167 fInvMassMinCut = gj.fInvMassMinCut ;
168 fFastRec = gj.fFastRec ;
e8daeb92 169 fOptionGJ = gj.fOptionGJ ;
1305bc01 170 fMinDistance = gj.fMinDistance ;
171 fOptFast = gj.fOptFast ;
172 fOnlyCharged = gj.fOnlyCharged ;
173 fOutputFile = gj.fOutputFile ;
174 fInputFileName = gj.fInputFileName ;
175 fOutputFileName = gj.fOutputFileName ;
176 fHIJINGFileName = gj.fHIJINGFileName ;
177 fHIJING = gj.fHIJING ;
178 fRatioMaxCut = gj.fRatioMaxCut ;
179 fRatioMinCut = gj.fRatioMinCut ;
180 fJetRatioMaxCut = gj.fJetRatioMaxCut ;
181 fJetRatioMinCut = gj.fJetRatioMinCut ;
182 fJetTPCRatioMaxCut = gj.fJetRatioMaxCut ;
183 fJetTPCRatioMinCut = gj.fJetRatioMinCut ;
184 fNEvent = gj.fNEvent ;
185 fNCone = gj.fNCone ;
186 fNPt = gj.fNPt ;
187 fResPara1 = gj.fResPara1 ;
188 fResPara2 = gj.fResPara2 ;
189 fResPara3 = gj.fResPara3 ;
190 fPtCut = gj.fPtCut ;
191 fNeutralPtCut = gj.fNeutralPtCut ;
192 fChargedPtCut = gj.fChargedPtCut ;
193 fPtJetSelectionCut = gj.fPtJetSelectionCut ;
194 fPhiMaxCut = gj.fPhiMaxCut ;
195 fPhiMinCut = gj.fPhiMinCut ;
196 fPosParaA = gj.fPosParaA ;
197 fPosParaB = gj.fPosParaB ;
198 fSelect = gj.fSelect ;
199 fTPCCutsLikeEMCAL = gj.fTPCCutsLikeEMCAL ;
200 fCone = gj.fCone ;
201 fPtThreshold = gj.fPtThreshold ;
202
dffafdc2 203 fDirName = gj.fDirName ;
204 fESDTree = gj.fESDTree ;
205 fPattern = gj.fPattern ;
206
1305bc01 207 SetName (gj.GetName()) ;
208 SetTitle(gj.GetTitle()) ;
209
210 for(Int_t i = 0; i<10; i++){
211 fCones[i] = gj.fCones[i] ;
212 fNameCones[i] = gj.fNameCones[i] ;
213 fPtThres[i] = gj.fPtThres[i] ;
214 fNamePtThres[i] = gj.fNamePtThres[i] ;
215 if( i < 6 ){
216 fJetXMin1[i] = gj.fJetXMin1[i] ;
217 fJetXMin2[i] = gj.fJetXMin2[i] ;
218 fJetXMax1[i] = gj.fJetXMax1[i] ;
219 fJetXMax2[i] = gj.fJetXMax2[i] ;
220 fBkgMean[i] = gj.fBkgMean[i] ;
221 fBkgRMS[i] = gj.fBkgRMS[i] ;
222 if( i < 2 ){
223 fJetE1[i] = gj.fJetE1[i] ;
224 fJetE2[i] = gj.fJetE2[i] ;
225 fJetSigma1[i] = gj.fJetSigma1[i] ;
226 fJetSigma2[i] = gj.fJetSigma2[i] ;
227 fPhiEMCALCut[i] = gj.fPhiEMCALCut[i] ;
228 }
229 }
230 }
77414c90 231}
232
233//____________________________________________________________________________
1305bc01 234AliPHOSGammaJet::~AliPHOSGammaJet()
235{
236 fOutputFile->Close() ;
237
77414c90 238}
239
240//____________________________________________________________________________
1305bc01 241void AliPHOSGammaJet::AddHIJINGToList(Int_t iEvent, TClonesArray * particleList,
242 TClonesArray * plCh,
243 TClonesArray * plNe,
dffafdc2 244 TClonesArray * plNePHOS)
77414c90 245{
e957fea8 246
1305bc01 247 // List of particles copied to a root file.
248
249// Char_t sb[] = "bgrd/";
250// // cout<<sb<<endl;
251// Char_t si[10];
252// Char_t sf[] = "/list.root";
253// // cout<<sf<<endl;
254// sprintf(si,"%d",iEvent);
255// strcat(sb,si) ;
256// // cout<<si<<endl;
257// strcat(sb,sf) ;
258// // cout<<si<<endl;
259// TFile * f = TFile::Open(sb) ;
260// //cout<<f->GetName()<<endl;
261
262 Char_t fi[100];
263 sprintf(fi,"bgrd/%d/list.root",iEvent);
264 TFile * f = TFile::Open(fi) ;
e8daeb92 265 //cout<<f->GetName()<<endl;
1305bc01 266
267 TParticle *particle = new TParticle();
e8daeb92 268 TTree *t = (TTree*) f->Get("T");
269 TBranch *branch = t->GetBranch("primaries");
1305bc01 270 branch->SetAddress(&particle);
271
272 Int_t index = particleList->GetEntries() ;
273 Int_t indexNe = plNe->GetEntries() ;
274 Int_t indexCh = plCh->GetEntries() ;
275 Int_t indexNePHOS = plNePHOS->GetEntries() ;
276 Double_t charge = 0.;
277 Int_t iParticle = 0 ;
278 Int_t m = 0;
279 Double_t x = 0., z = 0.;
e8daeb92 280 // cout<<"bkg entries "<<t->GetEntries()<<endl;
dffafdc2 281
282 AliPHOSGetter * gime = AliPHOSGetter::Instance() ;
283 const AliPHOSGeometry * geom = gime->PHOSGeometry() ;
284
1305bc01 285 if(!fOptFast){
e8daeb92 286 for (iParticle=0 ; iParticle < t->GetEntries() ; iParticle++) {
287 t->GetEvent(iParticle) ;
1305bc01 288
289 m = 0 ;
290 x = 0. ;
291 z = 0. ;
292
293 charge = TDatabasePDG::Instance()
294 ->GetParticle(particle->GetPdgCode())->Charge();
295
296 if((charge != 0) && (particle->Pt() > fChargedPtCut)){
297 if(TMath::Abs(particle->Eta())<fEtaCut){
298 new((*particleList)[index]) TParticle(*particle) ;
299 (dynamic_cast<TParticle*>(particleList->At(index)))
300 ->SetStatusCode(0) ;
301 index++ ;
302
303 new((*plCh)[indexCh]) TParticle(*particle) ;
304 (dynamic_cast<TParticle*>(plCh->At(indexCh)))->SetStatusCode(0) ;
305 indexCh++ ;
306
e8daeb92 307 if(strstr(fOptionGJ,"deb all")||strstr(fOptionGJ,"deb")){
1305bc01 308 dynamic_cast<TH1F*>(fListHistos->FindObject("PtSpectra"))->
309 Fill(particle->Pt());
310 }
311 }
312 else if((charge == 0) && (particle->Pt() > fNeutralPtCut) ){
313 geom->ImpactOnEmc(particle->Theta(),particle->Phi(), m,z,x);
314
315 if(m != 0)
316 {//Is in PHOS
e8daeb92 317 if(strstr(fOptionGJ,"deb all")|| strstr(fOptionGJ,"deb"))
1305bc01 318 dynamic_cast<TH1F*>(fListHistos->FindObject("PtSpectra"))->
319 Fill(particle->Pt());
320
321 new((*plNePHOS)[indexNePHOS]) TParticle(*particle) ;
322 (dynamic_cast<TParticle*>(plNePHOS->At(indexNePHOS)))->SetStatusCode(0) ;
323 indexNePHOS++ ;
324 }
325
326 if((particle->Phi()>fPhiEMCALCut[0]) &&
327 (particle->Phi()<fPhiEMCALCut[1]) && m == 0)
328 {//Is in EMCAL
e8daeb92 329 if(strstr(fOptionGJ,"deb all")|| strstr(fOptionGJ,"deb"))
1305bc01 330 dynamic_cast<TH1F*>(fListHistos->FindObject("PtSpectra"))->
331 Fill(particle->Pt());
332
333 new((*particleList)[index]) TParticle(*particle) ;
334 (dynamic_cast<TParticle*>(particleList->At(index)))
335 ->SetStatusCode(0) ;
336 index++ ;
337 new((*plNe)[indexNe]) TParticle(*particle) ;
338 (dynamic_cast<TParticle*>(plNe->At(indexNe)))->SetStatusCode(0) ;
339 indexNe++ ;
340 }
341 }
342 }
343 }
344 } //No OptFast
345 else{
346
347 Double_t mass = TDatabasePDG::Instance()->GetParticle(111)->Mass();
77414c90 348 TLorentzVector pPi0, pGamma1, pGamma2 ;
1305bc01 349 Double_t angle = 0, cellDistance = 0.;
350 Bool_t p1 = kFALSE;
351
352 // fFastRec = new AliPHOSFastGlobalReconstruction(fHIJINGFileName);
353 // fFastRec->FastReconstruction(iiEvent);
354
e8daeb92 355 for (iParticle=0 ; iParticle < t->GetEntries() ; iParticle++) {
356 t->GetEvent(iParticle) ;
1305bc01 357 m = 0 ;
358 x = 0. ;
359 z = 0. ;
360 charge = TDatabasePDG::Instance()
361 ->GetParticle(particle->GetPdgCode())->Charge();
362
363 if((charge != 0) && (particle->Pt() > fChargedPtCut)
364 && (TMath::Abs(particle->Eta())<fEtaCut)){
365
366 new((*particleList)[index]) TParticle(*particle) ;
367 (dynamic_cast<TParticle*>(particleList->At(index)))
368 ->SetStatusCode(0) ;
369 index++ ;
370
371 new((*plCh)[indexCh]) TParticle(*particle) ;
372 (dynamic_cast<TParticle*>(plCh->At(indexCh)))->SetStatusCode(0) ;
373 indexCh++ ;
374 }
375 else if(charge == 0){
376 geom->ImpactOnEmc(particle->Theta(),particle->Phi(), m,z,x);
377 if((particle->GetPdgCode() != 111) && (particle->Pt() > fNeutralPtCut) &&
378 (TMath::Abs(particle->Eta())<fEtaCut) ){
379 TLorentzVector part(particle->Px(),particle->Py(),
380 particle->Pz(),particle->Energy());
381 MakePhoton(part);
382 if(part.Pt() > fNeutralPtCut){
383
384 if(particle->Phi()>fPhiEMCALCut[0] &&
385 particle->Phi()<fPhiEMCALCut[1] && m == 0)
386 {
387 new((*particleList)[index]) TParticle(*particle) ;
388 (dynamic_cast<TParticle*>(particleList->At(index)))
389 ->SetStatusCode(0) ;
390 (dynamic_cast<TParticle*>(particleList->At(index)))
391 ->SetMomentum(part.Px(),part.Py(),part.Pz(),part.E());
392 index++ ;
393
394
395 new((*plNe)[indexNe]) TParticle(*particle) ;
396 (dynamic_cast<TParticle*>(plNe->At(indexNe)))
397 ->SetMomentum(part.Px(),part.Py(),part.Pz(),part.E());
398 (dynamic_cast<TParticle*>(plNe->At(indexNe)))->SetStatusCode(0) ;
399 indexNe++ ;
400 }
401 if(m != 0)
402 {
403 new((*plNePHOS)[indexNePHOS]) TParticle(*particle) ;
404 (dynamic_cast<TParticle*>(plNePHOS->At(indexNePHOS)))
405 ->SetMomentum(part.Px(),part.Py(),part.Pz(),part.E());
406 (dynamic_cast<TParticle*>(plNePHOS->At(indexNePHOS)))->SetStatusCode(0) ;
407 indexNePHOS++ ;
408 }
409 }
410 }
411 if((particle->GetPdgCode() == 111) && (particle->Pt() > fNeutralPtCut) &&
412 (TMath::Abs(particle->Eta())<fEtaCut+1))
413 {
414
415 pPi0.SetPxPyPzE(particle->Px(),particle->Py(),particle->Pz(),
416 particle->Energy());
417
418 //Decay
419
420 Pi0Decay(mass,pPi0,pGamma1,pGamma2,angle);
421 //Check if decay photons are too close for PHOS
422 cellDistance = angle*460; //cm
423 if (cellDistance < fMinDistance) {
e8daeb92 424 if(strstr(fOptionGJ,"deb all")|| strstr(fOptionGJ,"deb"))
1305bc01 425 dynamic_cast<TH1F*>(fListHistos->FindObject("PtSpectra"))->
426 Fill(particle->Pt());
427
428 //Pi0 inside phi EMCAL acceptance
429
430 TLorentzVector part(particle->Px(),particle->Py(),
431 particle->Pz(),particle->Energy());
432 MakePhoton(part);
433 if(part.Pt() > fNeutralPtCut){
434 if(particle->Phi()>fPhiEMCALCut[0] &&
435 particle->Phi()<fPhiEMCALCut[1] && m == 0){
436
437 new((*particleList)[index]) TParticle(*particle) ;
438 (dynamic_cast<TParticle*>(particleList->At(index)))->SetStatusCode(0) ;
439 (dynamic_cast<TParticle*>(particleList->At(index)))
440 ->SetMomentum(part.Px(),part.Py(),part.Pz(),part.E());
441 index++ ;
442
443 new((*plNe)[indexNe]) TParticle(*particle) ;
444 (dynamic_cast<TParticle*>(plNe->At(indexNe))) ->SetStatusCode(0) ;
445 (dynamic_cast<TParticle*>(plNe->At(indexNe)))
446 ->SetMomentum(part.Px(),part.Py(),part.Pz(),part.E());
447 indexNe++ ;
448 }
449 if(m != 0){
e8daeb92 450 if(strstr(fOptionGJ,"deb all")|| strstr(fOptionGJ,"deb"))
1305bc01 451 dynamic_cast<TH1F*>(fListHistos->FindObject("PtSpectra"))->
452 Fill(particle->Pt());
453 new((*plNePHOS)[indexNePHOS]) TParticle(*particle) ;
454 (dynamic_cast<TParticle*>(plNePHOS->At(indexNePHOS))) ->SetStatusCode(0) ;
455 (dynamic_cast<TParticle*>(plNePHOS->At(indexNePHOS)))
456 ->SetMomentum(part.Px(),part.Py(),part.Pz(),part.E());
457 indexNePHOS++;
458 }//In PHOS
459 }
460 }// if cell<distance
461 else {
462
463 dynamic_cast<TH2F*>(fListHistos->FindObject("AnglePair"))
464 ->Fill(pPi0.E(),angle);
465
466 p1 = kFALSE;
467 if(pGamma1.Pt() > 0. && TMath::Abs(pGamma1.Eta())<fEtaCut){
468
469 MakePhoton(pGamma1);
470
471 if(pGamma1.Pt() > fNeutralPtCut ){
472
473 TParticle * photon1 =
474 new TParticle(22,1,0,0,0,0,pGamma1.Px(),pGamma1.Py(),
475 pGamma1.Pz(),pGamma1.E(),0,0,0,0);
476 geom->ImpactOnEmc(photon1->Theta(),photon1->Phi(), m,z,x);
477 if( photon1->Phi()>fPhiEMCALCut[0] && photon1->Phi()<fPhiEMCALCut[1]
478 && m == 0){
e8daeb92 479 if(strstr(fOptionGJ,"deb all") || strstr(fOptionGJ,"deb"))
1305bc01 480 dynamic_cast<TH1F*>(fListHistos->FindObject("PtSpectra"))->
481 Fill(photon1->Pt());
482 new((*particleList)[index]) TParticle(*photon1) ;
483 (dynamic_cast<TParticle*>(particleList->At(index)))->SetStatusCode(0) ;
484 index++ ;
485
486 new((*plNe)[indexNe]) TParticle(*photon1) ;
487 (dynamic_cast<TParticle*>(plNe->At(indexNe))) ->SetStatusCode(0) ;
488 indexNe++ ;
489 p1 = kTRUE;
490 }
491 if(m != 0){
e8daeb92 492 if(strstr(fOptionGJ,"deb all") || strstr(fOptionGJ,"deb"))
1305bc01 493 dynamic_cast<TH1F*>(fListHistos->FindObject("PtSpectra"))->
494 Fill(photon1->Pt());
495 new((*plNePHOS)[indexNePHOS]) TParticle(*photon1) ;
496 (dynamic_cast<TParticle*>(plNePHOS->At(indexNePHOS)))->SetStatusCode(0) ;
497 indexNePHOS++;
498 p1 = kTRUE;
499 }//Is in PHOS
500 }
501 }
502 if(pGamma2.Pt() > 0. && TMath::Abs(pGamma2.Eta())<fEtaCut){
503
504 MakePhoton(pGamma2);
505 if(pGamma2.Pt() > fNeutralPtCut){
506
507 TParticle * photon2 =
508 new TParticle(22,1,0,0,0,0,pGamma2.Px(), pGamma2.Py(),
509 pGamma2.Pz(),pGamma2.E(),0,0,0,0);
510 geom->ImpactOnEmc(photon2->Theta(),photon2->Phi(), m,z,x);
511 if(photon2->Phi()>fPhiEMCALCut[0] &&
512 photon2->Phi()<fPhiEMCALCut[1] && m == 0){
e8daeb92 513 if(strstr(fOptionGJ,"deb all") || strstr(fOptionGJ,"deb"))
1305bc01 514 dynamic_cast<TH1F*>(fListHistos->FindObject("PtSpectra"))->
515 Fill(photon2->Pt());
516 new((*particleList)[index]) TParticle(*photon2) ;
517 (dynamic_cast<TParticle*>(particleList->At(index)))->SetStatusCode(0) ;
518 index++ ;
519
520 new((*plNe)[indexNe]) TParticle(*photon2) ;
521 (dynamic_cast<TParticle*>(plNe->At(indexNe))) ->SetStatusCode(0) ;
522 indexNe++ ;
523 }
524 if(m != 0){
e8daeb92 525 if(strstr(fOptionGJ,"deb all") || strstr(fOptionGJ,"deb"))
1305bc01 526 dynamic_cast<TH1F*>(fListHistos->FindObject("PtSpectra"))->
527 Fill(photon2->Pt());
528 new((*plNePHOS)[indexNePHOS]) TParticle(*photon2) ;
529 (dynamic_cast<TParticle*>(plNePHOS->At(indexNePHOS))) ->SetStatusCode(0) ;
530 indexNePHOS++;
531 }
532 if(p1){
533 // e = (pGamma1+pGamma2).E();
534 // if(IsAngleInWindow(angle,e))
535 dynamic_cast<TH2F*>
536 (fListHistos->FindObject("AnglePairAccepted"))->
537 Fill(pPi0.E(),angle);
538 }
539 }
540 }//photon2 in acceptance
541 }//if angle > mindist
542 }//if pi0
77414c90 543 }
77414c90 544 }//for (iParticle<nParticle)
1305bc01 545 }
546
547 //Info("AddHIJINGToList","End HIJING");
548}
549
550//____________________________________________________________________________
dffafdc2 551Double_t AliPHOSGammaJet::CalculateJetRatioLimit(const Double_t ptg,
1305bc01 552 const Double_t *par,
553 const Double_t *x) {
554
555 //Info("CalculateLimit","x1 %f, x2%f",x[0],x[1]);
e8daeb92 556 Double_t epp = par[0] + par[1] * ptg ;
557 Double_t spp = par[2] + par[3] * ptg ;
1305bc01 558 Double_t f = x[0] + x[1] * ptg ;
e8daeb92 559 Double_t epb = epp + par[4] ;
560 Double_t spb = TMath::Sqrt(spp*spp+ par[5]*par[5]) ;
561 Double_t rat = (epb - spb * f) / ptg ;
562 //Info("CalculateLimit","epp %f, spp %f, f %f", epp, spp, f);
563 //Info("CalculateLimit","epb %f, spb %f, rat %f", epb, spb, rat);
1305bc01 564 return rat ;
565}
566
567//____________________________________________________________________________
568void AliPHOSGammaJet::CreateParticleList(Int_t iEvent,
569 TClonesArray * particleList,
570 TClonesArray * plCh,
571 TClonesArray * plNe,
dffafdc2 572 TClonesArray * plNePHOS)
1305bc01 573{
574 //Info("CreateParticleList","Inside");
dffafdc2 575 AliPHOSGetter * gime = AliPHOSGetter::Instance(fInputFileName) ;
576 const AliPHOSGeometry * geom = gime->PHOSGeometry() ;
1305bc01 577 gime->Event(iEvent, "X") ;
578
dffafdc2 579
1305bc01 580 Int_t index = particleList->GetEntries() ;
581 Int_t indexCh = plCh->GetEntries() ;
582 Int_t indexNe = plNe->GetEntries() ;
583 Int_t indexNePHOS = plNePHOS->GetEntries() ;
584 Int_t iParticle = 0 ;
585 Double_t charge = 0.;
586 Int_t m = 0;
587 Double_t x = 0., z = 0.;
588 if(!fOptFast){
589
590 for (iParticle=0 ; iParticle < gime->NPrimaries() ; iParticle++) {
591 const TParticle * particle = gime->Primary(iParticle) ;
592
593 m = 0 ;
594 x = 0. ;
595 z = 0. ;
596
597 //Keep Stable particles within eta range
598 if((particle->GetStatusCode() == 1) &&
599 (particle->Pt() > 0)){
600 if(TMath::Abs(particle->Eta())<fEtaCut){
601 //Fill lists
602
603 charge = TDatabasePDG::Instance()
604 ->GetParticle(particle->GetPdgCode())->Charge();
605 if((charge != 0) && (particle->Pt() > fChargedPtCut)){
606
e8daeb92 607 if(strstr(fOptionGJ,"deb all")|| strstr(fOptionGJ,"deb"))
1305bc01 608 dynamic_cast<TH1F*>(fListHistos->FindObject("PtSpectra"))->
609 Fill(particle->Pt());
610 new((*plCh)[indexCh++]) TParticle(*particle) ;
611 new((*particleList)[index++]) TParticle(*particle) ;
612 }
613 else if((charge == 0) && (particle->Pt() > fNeutralPtCut)){
614 geom->ImpactOnEmc(particle->Theta(),particle->Phi(), m,z,x);
615 if(m != 0)
616 {//Is in PHOS
e8daeb92 617 if(strstr(fOptionGJ,"deb all")|| strstr(fOptionGJ,"deb"))
1305bc01 618 dynamic_cast<TH1F*>(fListHistos->FindObject("PtSpectra"))->
619 Fill(particle->Pt());
620
621 new((*plNePHOS)[indexNePHOS++]) TParticle(*particle) ;
622 }
623 if((particle->Phi()>fPhiEMCALCut[0]) &&
624 (particle->Phi()<fPhiEMCALCut[1]) && m == 0)
625 {//Is in EMCAL
e8daeb92 626 if(strstr(fOptionGJ,"deb all")|| strstr(fOptionGJ,"deb"))
1305bc01 627 dynamic_cast<TH1F*>(fListHistos->FindObject("PtSpectra"))->
628 Fill(particle->Pt());
629 new((*plNe)[indexNe++]) TParticle(*particle) ;
630 new((*particleList)[index++]) TParticle(*particle) ;
631 }
632 }
633 }
634 }//final particle, etacut
635 }//for (iParticle<nParticle)
636 }// No Op
637 else
638 {
639 Double_t mass = TDatabasePDG::Instance()->GetParticle(111)->Mass();
640 TLorentzVector pPi0, pGamma1, pGamma2 ;
641 Double_t angle = 0, cellDistance = 0.;
642
643 fFastRec = new AliPHOSFastGlobalReconstruction(fInputFileName);
644 fFastRec->FastReconstruction(iEvent);
645
646 for (iParticle=0 ; iParticle < gime->NPrimaries() ; iParticle++) {
647 const TParticle * particle = gime->Primary(iParticle) ;
648 m = 0 ;
649 x = 0. ;
650 z = 0. ;
651 //Keep Stable particles within eta range
652 if((particle->GetStatusCode() == 1) && (particle->Pt() > 0)){
653
654 //Fill lists
655
656 charge = TDatabasePDG::Instance()
657 ->GetParticle(particle->GetPdgCode())->Charge();
658 if((charge != 0) && (particle->Pt() > fChargedPtCut) && (TMath::Abs(particle->Eta())<fEtaCut)){
659 new((*plCh)[indexCh++]) TParticle(*particle) ;
660 new((*particleList)[index++]) TParticle(*particle) ;
661 }
662 else if(charge == 0) {
663 geom->ImpactOnEmc(particle->Theta(),particle->Phi(), m,z,x);
664 if((particle->GetPdgCode() != 111) && particle->Pt() > 0 &&
665 (TMath::Abs(particle->Eta())<fEtaCut))
666 {
667
668 TLorentzVector part(particle->Px(),particle->Py(),
669 particle->Pz(),particle->Energy());
670
671 MakePhoton(part);
672
673 if(part.Pt() > fNeutralPtCut){
674 if(particle->Phi()>fPhiEMCALCut[0] &&
675 particle->Phi()<fPhiEMCALCut[1] && m == 0)
676 {
677 new((*particleList)[index]) TParticle(*particle) ;
678 (dynamic_cast<TParticle*>(particleList->At(index)))
679 ->SetMomentum(part.Px(),part.Py(),part.Pz(),part.E());
680 index++ ;
681
682 new((*plNe)[indexNe]) TParticle(*particle) ;
683 (dynamic_cast<TParticle*>(plNe->At(indexNe)))
684 ->SetMomentum(part.Px(),part.Py(),part.Pz(),part.E());
685 indexNe++ ;
686 }
687 if(m != 0)
688 {
689 new((*plNePHOS)[indexNePHOS]) TParticle(*particle) ;
690 (dynamic_cast<TParticle*>(plNePHOS->At(indexNePHOS)))
691 ->SetMomentum(part.Px(),part.Py(),part.Pz(),part.E());
692 indexNePHOS++ ;
693 }
694 }// Small pt
695 } //No Pi0
696 if((particle->GetPdgCode() == 111) && (particle->Pt() > fNeutralPtCut) &&
697 (TMath::Abs(particle->Eta())<fEtaCut+1))
698 {
699
700 pPi0.SetPxPyPzE(particle->Px(),particle->Py(),particle->Pz(),
701 particle->Energy());
702
703 //Decay
704
705 Pi0Decay(mass,pPi0,pGamma1,pGamma2,angle);
706 //Check if decay photons are too close for PHOS
707 cellDistance = angle*460; //cm
708
709 if (cellDistance < fMinDistance) {
710
711 //Pi0 inside phi EMCAL acceptance
712
713
714 TLorentzVector part(particle->Px(),particle->Py(),
715 particle->Pz(),particle->Energy());
716 MakePhoton(part);
717
718 if(part.Pt() > fNeutralPtCut){
719 if(particle->Phi()>fPhiEMCALCut[0] &&
720 particle->Phi()<fPhiEMCALCut[1] && m == 0){
e8daeb92 721 if(strstr(fOptionGJ,"deb all")|| strstr(fOptionGJ,"deb"))
1305bc01 722 dynamic_cast<TH1F*>(fListHistos->FindObject("PtSpectra"))->
723 Fill(particle->Pt());
724
725 new((*plNe)[indexNe]) TParticle(*particle) ;
726 (dynamic_cast<TParticle*>(plNe->At(indexNe)))
727 ->SetMomentum(part.Px(),part.Py(),part.Pz(),part.E());
728 new((*particleList)[index]) TParticle(*particle) ;
729 (dynamic_cast<TParticle*>(particleList->At(index)))
730 ->SetMomentum(part.Px(),part.Py(),part.Pz(),part.E());
731 index++;
732 indexNe++;
733 }//InEMCAL
734 if(m != 0){
e8daeb92 735 if(strstr(fOptionGJ,"deb all")|| strstr(fOptionGJ,"deb"))
1305bc01 736 dynamic_cast<TH1F*>(fListHistos->FindObject("PtSpectra"))->
737 Fill(particle->Pt());
738 new((*plNePHOS)[indexNePHOS]) TParticle(*particle) ;
739 (dynamic_cast<TParticle*>(plNePHOS->At(indexNePHOS)))
740 ->SetMomentum(part.Px(),part.Py(),part.Pz(),part.E());
741 indexNePHOS++;
742 }//In PHOS
743 }//Small Pt
744 }// if cell<distance
745 else {
746
747 dynamic_cast<TH2F*>(fListHistos->FindObject("AnglePair"))
748 ->Fill(pPi0.E(),angle);
749
750 Bool_t p1 = kFALSE;
751
752 if(pGamma1.Pt() > 0 && TMath::Abs(pGamma1.Eta())<fEtaCut){
753 MakePhoton(pGamma1);
754
755 if(pGamma1.Pt() > fNeutralPtCut){
756 TParticle * photon1 =
757 new TParticle(22,1,0,0,0,0,pGamma1.Px(),pGamma1.Py(),
758 pGamma1.Pz(),pGamma1.E(),0,0,0,0);
759 geom->ImpactOnEmc(photon1->Theta(),photon1->Phi(), m,z,x);
760 if( photon1->Phi()>fPhiEMCALCut[0] && photon1->Phi()<fPhiEMCALCut[1]
761 && m == 0){
e8daeb92 762 if(strstr(fOptionGJ,"deb all") || strstr(fOptionGJ,"deb"))
1305bc01 763 dynamic_cast<TH1F*>(fListHistos->FindObject("PtSpectra"))->
764 Fill(photon1->Pt());
765 new((*plNe)[indexNe++]) TParticle(*photon1) ;
766 new((*particleList)[index++]) TParticle(*photon1) ;
767 //photon1->Print();
768 p1 = kTRUE;
769 }
770 if(m != 0){
e8daeb92 771 if(strstr(fOptionGJ,"deb all") || strstr(fOptionGJ,"deb"))
1305bc01 772 dynamic_cast<TH1F*>(fListHistos->FindObject("PtSpectra"))->
773 Fill(photon1->Pt());
774 new((*plNePHOS)[indexNePHOS++]) TParticle(*photon1) ;
775 p1 = kTRUE;
776 }
777 }
778 }
779 if(pGamma2.Pt() > 0 && TMath::Abs(pGamma2.Eta())<fEtaCut){
780 MakePhoton(pGamma2);
781
782 if(pGamma2.Pt() > fNeutralPtCut ){
783 TParticle * photon2 =
784 new TParticle(22,1,0,0,0,0,pGamma2.Px(), pGamma2.Py(),
785 pGamma2.Pz(),pGamma2.E(),0,0,0,0);
786 geom->ImpactOnEmc(photon2->Theta(),photon2->Phi(), m,z,x);
787 if(photon2->Phi()>fPhiEMCALCut[0] &&
788 photon2->Phi()<fPhiEMCALCut[1] && m == 0){
e8daeb92 789 if(strstr(fOptionGJ,"deb all") || strstr(fOptionGJ,"deb"))
1305bc01 790 dynamic_cast<TH1F*>(fListHistos->FindObject("PtSpectra"))->
791 Fill(photon2->Pt());
792 new((*plNe)[indexNe++]) TParticle(*photon2) ;
793 new((*particleList)[index++]) TParticle(*photon2) ;
794 }
795 if(m != 0){
e8daeb92 796 if(strstr(fOptionGJ,"deb all") || strstr(fOptionGJ,"deb"))
1305bc01 797 dynamic_cast<TH1F*>(fListHistos->FindObject("PtSpectra"))->
798 Fill(photon2->Pt());
799 new((*plNePHOS)[indexNePHOS++]) TParticle(*photon2) ;
800 }
801
802 if(p1){
803 // Float_t e = (pGamma1+pGamma2).E();
804 // if(IsAngleInWindow(angle,e))
805 dynamic_cast<TH2F*>
806 (fListHistos->FindObject("AnglePairAccepted"))->
807 Fill(pPi0.E(),angle);
808 }
809 }
810 }//photon2 in acceptance
811 }//if angle > mindist
812 }//if pi0
813 }//If neutral
814 }//final particle, etacut
815 }//for (iParticle<nParticle)
816 }//OptFast
817 //gime->Delete() ;
818}
dffafdc2 819//____________________________________________________________________________
820void AliPHOSGammaJet::CreateParticleListFromESD(TClonesArray * pl,
821 TClonesArray * plCh,
822 TClonesArray * plNe,
823 TClonesArray * plNePHOS,
824 const AliESD * esd){
825 //, Int_t iEvent){
826 //Info("CreateParticleListFromESD","Inside");
827 //AliPHOSGetter * gime = AliPHOSGetter::Instance(fInputFileName) ;
828 //gime->Event(iEvent, "P") ;
829
830 Double_t v[3] ; //vertex ;
831 esd->GetVertex()->GetXYZ(v) ;
832 Int_t index = pl->GetEntries() ;
833
834 //PHOS
835 Int_t begphos = esd->GetFirstPHOSParticle();
836 Int_t endphos = esd->GetFirstPHOSParticle() + esd->GetNumberOfPHOSParticles() ;
837
838 Int_t indexNePHOS = plNePHOS->GetEntries() ;
839 Int_t nphP ;
840 //cout<<"Fill PHOS"<<endl;
841 //cout<<"begin "<<begphos<<" end "<<endphos<<endl;
842 for (nphP = begphos; nphP < endphos; nphP++) {//////////////PHOS track loop
843 AliESDtrack * track = esd->GetTrack(nphP) ; // retrieve track from esd
844
845 // PID of this track
846 Double_t pid[AliPID::kSPECIESN];
847 track->GetPHOSpid(pid);
848 TParticle * particle = new TParticle() ;
849 Double_t en = track->GetPHOSsignal() ;
850 Double_t * p = new Double_t();
851 track->GetPHOSposition(p) ;
852 TVector3 pos(p[0],p[1],p[2]) ;
853 Double_t phi = pos.Phi();
854 Double_t theta= pos.Theta();
855 Double_t px = en*TMath::Cos(phi)*TMath::Sin(theta);;
856 Double_t py = en*TMath::Sin(phi)*TMath::Sin(theta);
857 Double_t pz = en*TMath::Cos(theta);
858 particle->SetMomentum(px,py,pz,en) ;
859 //cout<<"pid "<<pid[AliPID::kPhoton]<<endl ;
860 if( pid[AliPID::kPhoton] > 0.75)
861 new((*plNePHOS)[indexNePHOS++]) TParticle(*particle) ;
862 }
863 //TPC
864 Int_t begtpc = 0 ;
865 Int_t endtpc = esd->GetNumberOfTracks() ;
866 Int_t indexCh = plCh->GetEntries() ;
867 //cout<<"Fill TPC"<<endl;
868 //cout<<"begin "<<begtpc<<" end "<<endtpc<<endl;
869 for (nphP = begtpc; nphP < endtpc; nphP++) {////////////// track loop
870 AliESDtrack * track = esd->GetTrack(nphP) ; // retrieve track from esd
871
872 // PID of this track
873 //Double_t pid[AliPID::kSPECIESN];
874 //track->GetEMCALpid(pid);
875 TParticle * particle = new TParticle() ;
876 Double_t en = track->GetTPCsignal() ;
877 //Double_t *p = 0;
878 //track->GetEMCALposition(p) ;
879 //TVector3 pos(p[0],p[1],p[2]) ;
880 //Double_t phi = pos.Phi();
881 //Double_t theta= pos.Theta();
882
883 TVector3 mom = track->P3() ;
884 Double_t px = mom.Px();
885 Double_t py = mom.Py();
886 Double_t pz = mom.Pz();
887
888 particle->SetMomentum(px,py,pz,en) ;
889 //if( pid[AliPID::kPhoton] > 0.75)
890 new((*plCh)[indexCh++]) TParticle(*particle) ;
891 new((*pl)[index++]) TParticle(*particle) ;
892
893 }
894
895 //EMCAL
896 Int_t begem = esd->GetFirstEMCALParticle();
897 Int_t endem = esd->GetFirstEMCALParticle() + esd->GetNumberOfEMCALParticles() ;
898 Int_t indexNe = plNe->GetEntries() ;
899 //cout<<"Fill EMCAL"<<endl;
900 //cout<<"begin "<<begem<<" end "<<endem<<endl;
901 for (nphP = begem; nphP < endem; nphP++) {//////////////PHOS track loop
902 AliESDtrack * track = esd->GetTrack(nphP) ; // retrieve track from esd
903
904 // PID of this track
905 Double_t pid[AliPID::kSPECIESN];
906 track->GetEMCALpid(pid);
907 Double_t en = track->GetEMCALsignal() ;
908 //cout<<"EMCAL signal "<<en<<endl;
909 Double_t *p = new Double_t();
910 track->GetEMCALposition(p) ;
911 TVector3 pos(p[0],p[1],p[2]) ;
912 Double_t phi = pos.Phi();
913 Double_t theta= pos.Theta();
914
915 Double_t px = en*TMath::Cos(phi)*TMath::Sin(theta);;
916 Double_t py = en*TMath::Sin(phi)*TMath::Sin(theta);
917 Double_t pz = en*TMath::Cos(theta);
918 Int_t pdg = 0;
919 //particle->SetMomentum(px,py,pz,en) ;
920 if( pid[AliPID::kPhoton] > 0.75) //This has to be fixen.
921 pdg = 22;
922 //particle->SetPdgCode(22);
923 else if( pid[AliPID::kPi0] > 0.75)
924 pdg = 111;
925 // particle->SetPdgCode(111);
926 pdg = 22 ; //Assume all photons, no PID for EMCAL available.
927 TParticle * particle = new TParticle(pdg, 1, -1, -1, -1, -1,
928 px, py, pz, en, v[0], v[1], v[2], 0);
929 //cout<<"pid; ph "<<pid[AliPID::kPhoton]<<" pi0 "<<pid[AliPID::kPi0]<<endl;
930 new((*plNe)[indexNe++]) TParticle(*particle) ;
931 new((*pl)[index++]) TParticle(*particle) ;
932 }
933 //Info("CreateParticleListFromESD","End Inside");
934
935}
1305bc01 936
937
938
939//____________________________________________________________________________
940void AliPHOSGammaJet::Exec(Option_t *option)
941{
942 // does the job
e8daeb92 943 fOptionGJ = option;
1305bc01 944 MakeHistos() ;
dffafdc2 945
946 AliESD * esd = 0;
947 TChain * t = 0 ;
948
949 if(fESDdata){
950 // Create chain of esd trees
951 const UInt_t kNevent = static_cast<const UInt_t>(GetNEvent()) ;
952 t = AliReadESD(kNevent, fDirName, fESDTree, fPattern) ;
953 if(!t) {
954 AliError("Could not create the TChain") ;
955 //break ;
956 }
957
958 // ESD
959 t->SetBranchAddress("ESD",&esd); // point to the container esd where to put the event from the esdTree
1305bc01 960
dffafdc2 961 }
962
1305bc01 963
964// AliGenPythia* pyth = (AliGenPythia*) gAlice->Generator();
965// pyth->Init();
966
967 TClonesArray * particleList = new TClonesArray("TParticle",1000);
968 TClonesArray * plCh = new TClonesArray("TParticle",1000);
969 TClonesArray * plNe = new TClonesArray("TParticle",1000);
970 TClonesArray * plNePHOS = new TClonesArray("TParticle",1000);
971
972 for (Int_t iEvent = 0 ; iEvent < fNEvent ; iEvent++) {
e8daeb92 973 if(strstr(fOptionGJ,"deb")||strstr(fOptionGJ,"deb all"))
1305bc01 974 Info("Exec", "Event %d", iEvent) ;
975
976 fRan.SetSeed(0);
977
978 Double_t phig = 0., phil = 0., phich = 0 , phipi = 0;
979 Double_t etag = 0., etal = 0., etach = 0., etapi = 0. ;
980 Double_t ptg = 0., ptl = 0., ptch = 0., ptpi = 0.;
981
982 TLorentzVector jet (0,0,0,0);
983 TLorentzVector jettpc(0,0,0,0);
1305bc01 984
dffafdc2 985 if(fESDdata){
1305bc01 986
dffafdc2 987 Int_t iNbytes = t->GetEntry(iEvent); // store event in esd
988 //cout<<"nbytes "<<iNbytes<<endl;
989 if ( iNbytes == 0 ) {
990 AliError("Empty TChain") ;
991 break ;
992 }
993 CreateParticleListFromESD(particleList, plCh,plNe,plNePHOS, esd); //,iEvent);
994 }
995 else{
996 CreateParticleList(iEvent, particleList, plCh,plNe,plNePHOS);
997
998 // TLorentzVector pyjet(0,0,0,0);
999
1000 // Int_t nJ, nJT;
1001 // Float_t jets[4][10];
1002 // pyth->SetJetReconstructionMode(1);
1003 // pyth->LoadEvent();
1004 // pyth->GetJets(nJ, nJT, jets);
1005
1006 // Float_t pxJ = jets[0][0];
1007 // Float_t pyJ = jets[1][0];
1008 // Float_t pzJ = jets[2][0];
1009 // Float_t eJ = jets[3][0];
1010 // pyjet.SetPxPyPzE(pxJ,pyJ,pzJ,eJ ) ;
1305bc01 1011
dffafdc2 1012 // if(nJT > 1){
1013 // //Info("Exec",">>>>>>>>>>Number of jets !!!! %d",nJT);
1014 // for (Int_t iJ = 1; iJ < nJT; iJ++) {
1015 // Float_t pxJ = jets[0][iJ];
1016 // Float_t pyJ = jets[1][iJ];
1017 // Float_t pzJ = jets[2][iJ];
1018 // Float_t eJ = jets[3][iJ];
1019 // pyjet.SetPxPyPzE(pxJ,pyJ,pzJ,eJ ) ;
1020 // //Info("Exec",">>>>>Pythia Jet: %d, Phi %f, Eta %f, Pt %f",
1021 // // iJ,pyjet.Phi(),pyjet.Eta(),pyjet.Pt());
1022 // }
1023
1024 // }
1025
1026 if(fHIJING)
1027 AddHIJINGToList(iEvent, particleList, plCh,plNe, plNePHOS);
1028 }
1305bc01 1029
e8daeb92 1030 Bool_t iIsInPHOS = kFALSE ;
1031 GetGammaJet(plNePHOS, ptg, phig, etag, iIsInPHOS) ;
1305bc01 1032
e8daeb92 1033 if(iIsInPHOS){
1305bc01 1034
1035 //Info("Exec"," In PHOS") ;
1036 dynamic_cast<TH1F*>(fListHistos->FindObject("NGamma"))->Fill(ptg);
1037 dynamic_cast<TH2F*>(fListHistos->FindObject("PhiGamma"))
1038 ->Fill(ptg,phig);
1039 dynamic_cast<TH2F*>(fListHistos->FindObject("EtaGamma"))
1040 ->Fill(ptg,etag);
e8daeb92 1041 if(strstr(fOptionGJ,"deb")||strstr(fOptionGJ,"deb all"))
1305bc01 1042 Info("Exec", "Gamma: pt %f, phi %f, eta %f", ptg,
1043 phig,etag) ;
1044
1045// cout<<"n charged "<<plCh->GetEntries()<<endl;
1046// cout<<"n neutral "<<plNe->GetEntries()<<endl;
1047// cout<<"n All "<<particleList->GetEntries()<<endl;
1048
1049 GetLeadingCharge(plCh, ptg, phig, ptch, etach, phich) ;
1050 GetLeadingPi0 (plNe, ptg, phig, ptpi, etapi, phipi) ;
1051
1052// cout<<"n2 charged "<<plCh->GetEntries()<<endl;
1053// cout<<"n2 neutral "<<plNe->GetEntries()<<endl;
1054// cout<<"n2 All "<<particleList->GetEntries()<<endl;
1055
1056
1057 //TPC+EMCAL
1058
1059 //Is the leading cone inside EMCAL?
1060 Bool_t insidech = kFALSE ;
1061 if((phich - fCone) > fPhiEMCALCut[0] &&
1062 (phich + fCone) < fPhiEMCALCut[1]){
1063 insidech = kTRUE ;
1064 }
1065 Bool_t insidepi = kFALSE ;
1066 if((phipi - fCone) > fPhiEMCALCut[0] &&
1067 (phipi + fCone) < fPhiEMCALCut[1]){
1068 insidepi = kTRUE ;
1069 }
1070
1071 if ((ptch > 0 || ptpi > 0)){
1072 if((ptch > ptpi) && insidech){
1073 phil = phich ;
1074 etal = etach ;
1075 ptl = ptch ;
1076 dynamic_cast<TH2F*>(fListHistos->FindObject("ChargeRatio"))
1077 ->Fill(ptg,ptch/ptg);
1078 dynamic_cast<TH2F*>(fListHistos->FindObject("DeltaPhiCharge"))
1079 ->Fill(ptg,phig-phich);
1080 dynamic_cast<TH2F*>(fListHistos->FindObject("DeltaEtaCharge"))
1081 ->Fill(ptg,etag-etach);
e8daeb92 1082 if(strstr(fOptionGJ,"deb"))
1305bc01 1083 Info("Exec"," Charged Leading") ;
1084 }
1085 if((ptpi > ptch) && insidepi){
1086 phil = phipi ;
1087 etal = etapi ;
1088 ptl = ptpi ;
1089
1090 dynamic_cast<TH2F*>(fListHistos->FindObject("Pi0Ratio"))
1091 ->Fill(ptg,ptpi/ptg);
1092 dynamic_cast<TH2F*>(fListHistos->FindObject("DeltaPhiPi0"))
1093 ->Fill(ptg,phig-phipi);
1094 dynamic_cast<TH2F*>(fListHistos->FindObject("DeltaEtaPi0"))
1095 ->Fill(ptg,etag-etapi);
1096
e8daeb92 1097 if(ptpi > 0. && strstr(fOptionGJ,"deb"))
1305bc01 1098 Info("Exec"," Pi0 Leading") ;
1099 }
1100
e8daeb92 1101 if(strstr(fOptionGJ,"deb"))
1305bc01 1102 Info("Exec","Leading pt %f, phi %f",ptl,phil);
1103 if(insidech || insidepi){
1104 if(!fAnyConeOrPt){
1105
1106 MakeJet(particleList, ptg, phig, ptl, phil, etal, "", jet);
1107
e8daeb92 1108 if(strstr(fOptionGJ,"deb")){
1305bc01 1109// Info("Exec","Pythia Jet: Phi %f, Eta %f, Pt %f",
1110// pyjet.Phi(),pyjet.Eta(),pyjet.Pt());
1111 Info("Exec","TPC+EMCAL Jet: Phi %f, Eta %f, Pt %f",
1112 jet.Phi(),jet.Eta(),jet.Pt());
1113 }
1114// dynamic_cast<TH2F*>(fListHistos->FindObject("DeltaPhiJet"))
1115// ->Fill(ptg,pyjet.Phi()-jet.Phi());
1116// dynamic_cast<TH2F*>(fListHistos->FindObject("DeltaEtaJet"))
1117// ->Fill(ptg,pyjet.Eta()-jet.Eta());
1118// dynamic_cast<TH2F*>(fListHistos->FindObject("DeltaPtJet"))
1119// ->Fill(ptg,pyjet.Pt()-jet.Pt());
1120 }
1121 else
1122 MakeJetAnyConeOrPt(particleList, ptg, phig, ptl, phil, etal, "");
1123 }
1124
1125 //TPC
1126 if(fOnlyCharged && ptch > 0.)
1127 {
e8daeb92 1128 if(strstr(fOptionGJ,"deb"))
1305bc01 1129 Info("Exec","Leading TPC pt %f, phi %f",ptch,phich);
1130
1131 dynamic_cast<TH2F*>(fListHistos->FindObject("TPCRatio"))
1132 ->Fill(ptg,ptch/ptg);
1133 dynamic_cast<TH2F*>(fListHistos->FindObject("DeltaPhiTPC"))
1134 ->Fill(ptg,phig-phich);
1135 dynamic_cast<TH2F*>(fListHistos->FindObject("DeltaEtaTPC"))
1136 ->Fill(ptg,etag-etach);
1137
1138 if(!fAnyConeOrPt){
1139
1140 MakeJet(plCh, ptg, phig, ptch, phich, etach, "TPC",jettpc);
1141
e8daeb92 1142 if(strstr(fOptionGJ,"deb")){
1305bc01 1143// Info("Exec","Pythia Jet: Phi %f, Eta %f, Pt %f",
1144// pyjet.Phi(),pyjet.Eta(),pyjet.Pt());
1145 Info("Exec","TPC Jet: Phi %f, Eta %f, Pt %f",
1146 jettpc.Phi(),jettpc.Eta(),jettpc.Pt());
1147 }
1148// dynamic_cast<TH2F*>(fListHistos->FindObject("DeltaPhiTPCJet"))
1149// ->Fill(ptg,pyjet.Phi()-jettpc.Phi());
1150// dynamic_cast<TH2F*>(fListHistos->FindObject("DeltaEtaTPCJet"))
1151// ->Fill(ptg,pyjet.Eta()-jettpc.Eta());
1152// dynamic_cast<TH2F*>(fListHistos->FindObject("DeltaPtTPCJet"))
1153// ->Fill(ptg,pyjet.Pt()-jettpc.Pt());
1154 }
1155 else
1156 MakeJetAnyConeOrPt(plCh, ptg, phig, ptch, phich, etach, "TPC");
1157
1158 }
1159 }
1160 }
1161
1162 particleList->Delete() ;
1163 plCh->Delete() ;
1164 plNe->Delete() ;
1165 plNePHOS->Delete() ;
77414c90 1166 }//loop: events
1305bc01 1167
1168 delete plNe ;
1169 delete plCh ;
1170 delete particleList ;
1171
1172 fOutputFile->Write() ;
1173 fOutputFile->cd();
1174 this->Write();
77414c90 1175}
1176
1305bc01 1177//____________________________________________________________________________
1178void AliPHOSGammaJet::FillJetHistos(TClonesArray * pl, Double_t ptg,
1179 TString conf, TString type)
1180{
e8daeb92 1181 //Fill jet fragmentation histograms if !fAnyCone,
1182 //only for fCone and fPtThres
1305bc01 1183 TParticle * particle = 0 ;
1184 Int_t ipr = -1 ;
1185 Float_t charge = 0;
1186
1187 TIter next(pl) ;
1188 while ( (particle = dynamic_cast<TParticle*>(next())) ) {
1189 ipr++ ;
1190 Double_t pt = particle->Pt();
1191
1192 charge = TDatabasePDG::Instance()
1193 ->GetParticle(particle->GetPdgCode())->Charge();
1194 if(charge != 0){//Only jet Charged particles
1195 dynamic_cast<TH2F*>
1196 (fListHistos->FindObject(type+conf+"Fragment"))
1197 ->Fill(ptg,pt/ptg);
1198 dynamic_cast<TH2F*>
1199 (fListHistos->FindObject(type+conf+"PtDist"))
1200 ->Fill(ptg,pt);
1201 }
1202 }
1203 if(type == "Bkg"){
1204 dynamic_cast<TH1F*>
1205 (fListHistos->FindObject("NBkg"+conf))->Fill(ipr);
1206 }
1207}
1208//____________________________________________________________________________
1209void AliPHOSGammaJet::FillJetHistosAnyConeOrPt(TClonesArray * pl, Double_t ptg,
1210 TString conf, TString type,
1211 TString cone, TString ptcut)
1212{
e8daeb92 1213 //Fill jet fragmentation histograms if fAnyCone,
1214 //for several cones and pt thresholds
1305bc01 1215 TParticle *particle = 0;
1216 Int_t ipr=-1;
1217 Float_t charge = 0;
1218
1219 TIter next(pl) ;
1220 while ( (particle = dynamic_cast<TParticle*>(next())) ) {
1221 ipr++;
1222 Double_t pt = particle->Pt();
1223 charge = TDatabasePDG::Instance()
1224 ->GetParticle(particle->GetPdgCode())->Charge();
1225 if(charge != 0){//Only jet Charged particles
1226 dynamic_cast<TH2F*>
1227 (fListHistos->FindObject(type+conf+"FragmentCone"+cone+"Pt"+ptcut))
1228 ->Fill(ptg,pt/ptg);
1229 dynamic_cast<TH2F*>
1230 (fListHistos->FindObject(type+conf+"PtDistCone"+cone+"Pt"+ptcut))
1231 ->Fill(ptg,pt);
1232 }
1233 }//while
1234
1235 if(type == "Bkg"){
1236 dynamic_cast<TH1F*>
1237 (fListHistos->FindObject("NBkg"+conf+"Cone"+cone+"Pt"+ptcut))
1238 ->Fill(ipr);
1239 }
1240}
1241
1242//____________________________________________________________________________
1243void AliPHOSGammaJet::GetGammaJet(TClonesArray * pl, Double_t &pt,
e8daeb92 1244 Double_t &phi, Double_t &eta, Bool_t &Is) const
1305bc01 1245{
e8daeb92 1246 //Search for the prompt photon in PHOS with pt > fPtCut
1305bc01 1247 pt = -10.;
1248 eta = -10.;
1249 phi = -10.;
1250
1251 for(Int_t ipr = 0;ipr < pl->GetEntries() ; ipr ++ ){
1252 TParticle * particle = dynamic_cast<TParticle *>(pl->At(ipr)) ;
1253
dffafdc2 1254 if((particle->Pt() > fPtCut) && (particle->Pt() > pt)){
1255
1256 pt = particle->Pt();
1257 phi = particle->Phi() ;
1258 eta = particle->Eta() ;
1259 Is = kTRUE;
1305bc01 1260 }
1261 }
1262}
1263
1264//____________________________________________________________________________
1265void AliPHOSGammaJet::GetLeadingCharge(TClonesArray * pl,
1266 Double_t ptg, Double_t phig,
e8daeb92 1267 Double_t &pt, Double_t &eta, Double_t &phi) const
1305bc01 1268{
e8daeb92 1269 //Search for the charged particle with highest with
1270 //Phi=Phi_gamma-Pi and pT=0.1E_gamma
1305bc01 1271 pt = -100.;
1272 eta = -100;
1273 phi = -100;
1274
1275 for(Int_t ipr = 0;ipr < pl->GetEntries() ; ipr ++ ){
1276
1277 TParticle * particle = dynamic_cast<TParticle *>(pl->At(ipr)) ;
1278
1279 Double_t ptl = particle->Pt();
1280 Double_t rat = ptl/ptg ;
1281 Double_t phil = particle->Phi() ;
1282
1283 if(((phig-phil)> fPhiMinCut) && ((phig-phil)<fPhiMaxCut) &&
1284 (rat > fRatioMinCut) && (rat < fRatioMaxCut) && (ptl > pt)) {
1285 eta = particle->Eta() ;
1286 phi = phil ;
1287 pt = ptl ;
1288 //printf("GetLeadingCharge: %f %f %f %f \n", pt, eta, phi,rat) ;
1289 }
1290 }
1291 //printf("GetLeadingCharge: %f %f %f \n", pt, eta, phi) ;
1292
1293}
1294
1295
1296//____________________________________________________________________________
1297void AliPHOSGammaJet::GetLeadingPi0(TClonesArray * pl,
1298 Double_t ptg, Double_t phig,
e8daeb92 1299 Double_t &pt, Double_t &eta, Double_t &phi)
1305bc01 1300{
e8daeb92 1301
1302 //Search for the neutral pion with highest with
1303 //Phi=Phi_gamma-Pi and pT=0.1E_gamma
1305bc01 1304 pt = -100.;
1305 eta = -100.;
1306 phi = -100.;
1307 Double_t ptl = -100.;
1308 Double_t rat = -100.;
1309 Double_t phil = -100. ;
1310
1311 TIter next(pl);
1312 TParticle * particle = 0;
1313 Float_t ef = 0;
1314 if(!fOptFast){
1315 Float_t e = 0;
1316 while ( (particle = (TParticle*)next()) ) {
1317 if( particle->GetPdgCode() == 111){
1318 ptl = particle->Pt();
1319 rat = ptl/ptg ;
1320 phil = particle->Phi() ;
1321 e = particle->Energy();
1322 dynamic_cast<TH2F*>
1323 (fListHistos->FindObject("AnglePairNoCut"))->
1324 Fill(e,0.1);
1325 dynamic_cast<TH2F*>
1326 (fListHistos->FindObject("InvMassPairNoCut"))->
1327 Fill(ptg,1.35);
1328
1329 if(((phig-phil)> fPhiMinCut) && ((phig-phil)<fPhiMaxCut) &&
1330 (rat > fRatioMinCut) && (rat < fRatioMaxCut)) {
1331
1332 dynamic_cast<TH2F*>
1333 (fListHistos->FindObject("AnglePairLeadingCut"))->
1334 Fill(e,0.1);
1335 dynamic_cast<TH2F*>
1336 (fListHistos->FindObject("InvMassPairLeadingCut"))->
1337 Fill(ptg,1.35);
1338
1339 dynamic_cast<TH2F*>
1340 (fListHistos->FindObject("AnglePairAngleCut"))->
1341 Fill(e,0.15);
1342 dynamic_cast<TH2F*>
1343 (fListHistos->FindObject("InvMassPairAngleCut"))->
1344 Fill(ptg,1.36);
1345
1346 dynamic_cast<TH2F*>
1347 (fListHistos->FindObject("InvMassPairAllCut"))->
1348 Fill(ptg,0.27);
1349 dynamic_cast<TH2F*>
1350 (fListHistos->FindObject("AnglePairAllCut"))->
1351 Fill(e,1.34);
1352
1353
1354 if(ptl > pt){
1355 eta = particle->Eta() ;
1356 phi = phil ;
1357 pt = ptl ;
1358 ef = e;
1359 //printf("GetLeadingPi0: %f %f %f %f %f \n", pt, eta, phi, rat, ptg) ;
1360 }
1361 }
1362 }
1363
1364 dynamic_cast<TH2F*>
1365 (fListHistos->FindObject("InvMassPairLeading"))->
1366 Fill(ptg,1.35);
1367 dynamic_cast<TH2F*>
1368 (fListHistos->FindObject("AnglePairLeading"))->
1369 Fill(ef,0.1);
1370 }
1371 }//No fOptfast
1372 else{
1373 Int_t iPrimary = -1;
1374 TLorentzVector gammai,gammaj;
1375 Double_t angle = 0., e = 0., invmass = 0.;
1376 Double_t anglef = 0., ef = 0., invmassf = 0.;
1377 Int_t ksPdg = 0;
1378 Int_t jPrimary=-1;
1379
1380 while ( (particle = (TParticle*)next()) ) {
1381 iPrimary++;
dffafdc2 1382
1305bc01 1383 ksPdg = particle->GetPdgCode();
1384 ptl = particle->Pt();
1385 if(ksPdg == 111){ //2 gamma
1386 rat = ptl/ptg ;
1387 phil = particle->Phi() ;
1388 if((ptl> pt)&& (rat > fRatioMinCut) && (rat < fRatioMaxCut) &&
1389 ((phig-phil)>fPhiMinCut)&&((phig-phil)<fPhiMaxCut)){
1390 eta = particle->Eta() ;
1391 phi = phil ;
1392 pt = ptl ;
1393 }// cuts
1394 }// pdg = 111
1395 if(ksPdg == 22){//1 gamma
1396 particle->Momentum(gammai);
1397 jPrimary=-1;
1398 TIter next2(pl);
1399 while ( (particle = (TParticle*)next2()) ) {
1400 jPrimary++;
1401 if(jPrimary>iPrimary){
1402 ksPdg = particle->GetPdgCode();
1403 if(ksPdg == 22){
1404 particle->Momentum(gammaj);
1405 //Info("GetLeadingPi0","Egammai %f, Egammaj %f",
1406 //gammai.Pt(),gammaj.Pt());
1407
1408 ptl = (gammai+gammaj).Pt();
1409 phil = (gammai+gammaj).Phi();
1410 if(phil < 0)
1411 phil+=TMath::TwoPi();
1412 rat = ptl/ptg ;
1413 invmass = (gammai+gammaj).M();
1414 angle = gammaj.Angle(gammai.Vect());
1415 //Info("GetLeadingPi0","Angle %f", angle);
1416 e = (gammai+gammaj).E();
1417
1418 dynamic_cast<TH2F*>
1419 (fListHistos->FindObject("AnglePairNoCut"))->
1420 Fill(e,angle);
1421 dynamic_cast<TH2F*>
1422 (fListHistos->FindObject("InvMassPairNoCut"))->
1423 Fill(ptg,invmass);
1424
1425 if((rat > fRatioMinCut) && (rat < fRatioMaxCut) &&
1426 ((phig-phil)>fPhiMinCut)&&((phig-phil)<fPhiMaxCut)){
1427
1428 dynamic_cast<TH2F*>
1429 (fListHistos->FindObject("AnglePairLeadingCut"))->
1430 Fill(e,angle);
1431 dynamic_cast<TH2F*>
1432 (fListHistos->FindObject("InvMassPairLeadingCut"))->
1433 Fill(ptg,invmass);
1434
1435 if(IsAngleInWindow(angle,e)){
1436 dynamic_cast<TH2F*>
1437 (fListHistos->FindObject("AnglePairAngleCut"))->
1438 Fill(e,angle);
1439 dynamic_cast<TH2F*>
1440 (fListHistos->FindObject("InvMassPairAngleCut"))->
1441 Fill(ptg,invmass);
1442
1443 //Info("GetLeadingPi0","InvMass %f", invmass);
1444 if((invmass>fInvMassMinCut) && (invmass<fInvMassMaxCut)){
1445 dynamic_cast<TH2F*>
1446 (fListHistos->FindObject("InvMassPairAllCut"))->
1447 Fill(ptg,invmass);
1448 dynamic_cast<TH2F*>
1449 (fListHistos->FindObject("AnglePairAllCut"))->
1450 Fill(e,angle);
1451 if(ptl > pt ){
1452 pt = ptl;
1453 eta = particle->Eta() ;
1454 phi = phil ;
1455 ef = e ;
1456 anglef = angle ;
1457 invmassf = invmass ;
1458
1459 }
1460 }//cuts
1461 }//(invmass>0.125) && (invmass<0.145)
1462 }//gammaj.Angle(gammai.Vect())<0.04
1463 }//if pdg = 22
1464 }//iprimary<jprimary
1465 }//while
1466 }// if pdg = 22
1467 // }
1468 }//while
1469
1470 if(ef > 0.){
1471 dynamic_cast<TH2F*>
1472 (fListHistos->FindObject("InvMassPairLeading"))->
1473 Fill(ptg,invmassf);
1474 dynamic_cast<TH2F*>
1475 (fListHistos->FindObject("AnglePairLeading"))->
1476 Fill(ef,anglef);
1477 }
1478 }//fOptFast
1479 // printf("GetLeadingPi0: %f %f %f \n", pt, eta, phi) ;
1480}
1481
1482
1483//____________________________________________________________________________
1484void AliPHOSGammaJet::InitParameters()
1485{
e8daeb92 1486
1487 //Initialize the parameters of the analysis.
1488
1305bc01 1489 fAngleMaxParam.Set(4) ;
1490 fAngleMaxParam.AddAt(0.4,0);//={0.4,-0.25,0.025,-2e-4};
1491 fAngleMaxParam.AddAt(-0.25,1) ;
1492 fAngleMaxParam.AddAt(0.025,2) ;
1493 fAngleMaxParam.AddAt(-2e-4,3) ;
1494 fAnyConeOrPt = kFALSE ;
1495 fOutputFileName = "GammaJet.root" ;
e8daeb92 1496 fOptionGJ = "";
1305bc01 1497 fHIJINGFileName = "galice.root" ;
1498 fHIJING = kFALSE ;
1499 fMinDistance = 3.6 ;
1500 fEtaCut = 0.7 ;
1501 fInvMassMaxCut = 0.15 ;
1502 fInvMassMinCut = 0.12 ;
1503 fOnlyCharged = kFALSE ;
1504 fOptFast = kFALSE ;
dffafdc2 1505 fESDdata = kTRUE ;
1305bc01 1506 fPhiEMCALCut[0] = 60. *TMath::Pi()/180.;
1507 fPhiEMCALCut[1] = 180.*TMath::Pi()/180.;
1508 fPhiMaxCut = 3.4 ;
1509 fPhiMinCut = 2.9 ;
1510 fPtCut = 10. ;
1511 fNeutralPtCut = 0.5 ;
1512 fChargedPtCut = 0.5 ;
1513 fTPCCutsLikeEMCAL = kFALSE ;
1514 //Jet selection parameters
1515 //Fixed cut (old)
1516 fRatioMaxCut = 1.0 ;
1517 fRatioMinCut = 0.1 ;
1518 fJetRatioMaxCut = 1.2 ;
1519 fJetRatioMinCut = 0.8 ;
1520 fJetTPCRatioMaxCut = 1.2 ;
1521 fJetTPCRatioMinCut = 0.3 ;
1522 fSelect = kFALSE ;
1523
dffafdc2 1524 fDirName = "./" ;
1525 fESDTree = "esdTree" ;
1526 fPattern = "." ;
1527
1305bc01 1528 //Cut depending on gamma energy
1529
1530 fPtJetSelectionCut = 20.; //For Low pt jets+BKG, another limits applyed
1531 //Reconstructed jet energy dependence parameters
1532 //e_jet = a1+e_gamma b2.
1533 //Index 0-> Pt>2 GeV r = 0.3; Index 1-> Pt>0.5 GeV r = 0.3
1534 fJetE1[0] = -5.75; fJetE1[1] = -4.1;
1535 fJetE2[0] = 1.005; fJetE2[1] = 1.05;
1536
1537 //Reconstructed sigma of jet energy dependence parameters
1538 //s_jet = a1+e_gamma b2.
1539 //Index 0-> Pt>2 GeV r = 0.3; Index 1-> Pt>0.5 GeV r = 0.3
1540 fJetSigma1[0] = 2.65; fJetSigma1[1] = 2.75;
1541 fJetSigma2[0] = 0.0018; fJetSigma2[1] = 0.033;
1542
1543 //Background mean energy and RMS
1544 //Index 0-> No BKG; Index 1-> BKG > 2 GeV;
1545 //Index 2-> (low pt jets)BKG > 0.5 GeV;
1546 //Index > 2, same for TPC conf
1547 fBkgMean[0] = 0.; fBkgMean[1] = 8.8 ; fBkgMean[2] = 69.5;
1548 fBkgMean[3] = 0.; fBkgMean[4] = 6.4; fBkgMean[5] = 48.6;
1549 fBkgRMS[0] = 0.; fBkgRMS[1] = 7.5; fBkgRMS[2] = 22.0;
1550 fBkgRMS[3] = 0.; fBkgRMS[4] = 5.4; fBkgRMS[5] = 13.2;
1551
1552 //Factor x of min/max = E -+ x * sigma. Obtained after selecting the
1553 //limits for monoenergetic jets.
1554 //Index 0-> No BKG; Index 1-> BKG > 2 GeV;
1555 //Index 2-> (low pt jets) BKG > 0.5 GeV;
1556 //Index > 2, same for TPC conf
1557
1558 fJetXMin1[0] =-0.69 ; fJetXMin1[1] = 0.39 ; fJetXMin1[2] =-0.88 ;
1559 fJetXMin1[3] =-2.0 ; fJetXMin1[4] =-0.442 ; fJetXMin1[5] =-1.1 ;
1560 fJetXMin2[0] = 0.066; fJetXMin2[1] = 0.038; fJetXMin2[2] = 0.034;
1561 fJetXMin2[3] = 0.25 ; fJetXMin2[4] = 0.113; fJetXMin2[5] = 0.077 ;
1562 fJetXMax1[0] =-3.8 ; fJetXMax1[1] =-0.76 ; fJetXMax1[2] =-3.6 ;
1563 fJetXMax1[3] =-2.7 ; fJetXMax1[4] =-1.21 ; fJetXMax1[5] =-3.7 ;
1564 fJetXMax2[0] =-0.012; fJetXMax2[1] =-0.022; fJetXMax2[2] = 0.016;
1565 fJetXMax2[3] =-0.024; fJetXMax2[4] =-0.008; fJetXMax2[5] = 0.027;
1566
1567
1568 //Photon fast reconstruction
1569 fResPara1 = 0.0255 ; // GeV
1570 fResPara2 = 0.0272 ;
1571 fResPara3 = 0.0129 ;
1572
1573 fPosParaA = 0.096 ; // cm
1574 fPosParaB = 0.229 ;
1575
1576 //Different cones and pt thresholds to construct the jet
1577
1578 fCone = 0.3 ;
1579 fPtThreshold = 0. ;
1580 fNCone = 1 ;
1581 fNPt = 1 ;
1582 fCones[0] = 0.3 ; fNameCones[0] = "03" ;
1583 fPtThres[0] = 0.5 ; fNamePtThres[0] = "05" ;
1584
1585}
1586
1587//__________________________________________________________________________-
dffafdc2 1588Bool_t AliPHOSGammaJet::IsAngleInWindow(const Float_t angle,const Float_t e) {
e8daeb92 1589 //Check if the opening angle of the candidate pairs is inside
1590 //our selection windowd
1305bc01 1591 Bool_t result = kFALSE;
1592 Double_t mpi0 = 0.1349766;
1593 Double_t max = fAngleMaxParam.At(0)*TMath::Exp(fAngleMaxParam.At(1)*e)
1594 +fAngleMaxParam.At(2)+fAngleMaxParam.At(3)*e;
1595 Double_t arg = (e*e-2*mpi0*mpi0)/(e*e);
1596 Double_t min = 100. ;
1597 if(arg>0.)
1598 min = TMath::ACos(arg);
1599
1600 if((angle<max)&&(angle>=min))
1601 result = kTRUE;
1602
1603 return result;
1604}
1605
1606//__________________________________________________________________________-
dffafdc2 1607Bool_t AliPHOSGammaJet::IsJetSelected(const Double_t ptg, const Double_t ptj,
1305bc01 1608 const TString type ){
e8daeb92 1609 //Check if the energy of the reconstructed jet is within an energy window
1305bc01 1610
1611 Double_t par[6];
1612 Double_t xmax[2];
1613 Double_t xmin[2];
1614
1615 Int_t iTPC = 0;
1616
1617 if(type == "TPC" && !fTPCCutsLikeEMCAL){
1618 iTPC = 3 ;//If(fTPCCutsLikeEMCAL) take jet energy cuts like EMCAL
1619 }
1620
1621
1622 if(!fHIJING){
1623 //Phythia alone, jets with pt_th > 0.2, r = 0.3
1624 par[0] = fJetE1[0]; par[1] = fJetE2[0];
1625 //Energy of the jet peak
1626 //e_jet = fJetE1[0]+fJetE2[0]*e_gamma, simulation fit
1627 par[2] = fJetSigma1[0]; par[3] = fJetSigma2[0];
1628 //Sigma of the jet peak
1629 //sigma_jet = fJetSigma1[0]+fJetSigma2[0]*e_gamma, simulation fit
1630 par[4] = fBkgMean[0 + iTPC]; par[5] = fBkgRMS[0 + iTPC];
1631 //Parameters reserved for HIJING bkg.
1632 xmax[0] = fJetXMax1[0 + iTPC]; xmax[1] = fJetXMax2[0 + iTPC];
1633 xmin[0] = fJetXMin1[0 + iTPC]; xmin[1] = fJetXMin2[0 + iTPC];
1634 //Factor that multiplies sigma to obtain the best limits,
1635 //by observation, of mono jet ratios (ptjet/ptg)
1636 //X_jet = fJetX1[0]+fJetX2[0]*e_gamma
1637
1638 }
1639 else{
1640 if(ptg > fPtJetSelectionCut){
1641 //Phythia +HIJING with pt_th > 2 GeV/c, r = 0.3
1642 par[0] = fJetE1[0]; par[1] = fJetE2[0];
1643 //Energy of the jet peak, same as in pp
1644 //e_jet = fJetE1[0]+fJetE2[0]*e_gamma, simulation fit
1645 par[2] = fJetSigma1[0]; par[3] = fJetSigma2[0];
1646 //Sigma of the jet peak, same as in pp
1647 //sigma_jet = fJetSigma1[0]+fJetSigma2[0]*e_gamma, simulation fit
1648 par[4] = fBkgMean[1 + iTPC]; par[5] = fBkgRMS[1 + iTPC];
1649 //Mean value and RMS of HIJING Bkg
1650 xmax[0] = fJetXMax1[1 + iTPC]; xmax[1] = fJetXMax2[1 + iTPC];
1651 xmin[0] = fJetXMin1[1 + iTPC]; xmin[1] = fJetXMin2[1 + iTPC];
1652 //Factor that multiplies sigma to obtain the best limits,
1653 //by observation, of mono jet ratios (ptjet/ptg) mixed with HIJING Bkg,
1654 //pt_th > 2 GeV, r = 0.3
1655 //X_jet = fJetX1[0]+fJetX2[0]*e_gamma
1656
1657 }
1658 else{
1659 //Phythia + HIJING with pt_th > 0.5 GeV/c, r = 0.3
1660 par[0] = fJetE1[1]; par[1] = fJetE2[1];
1661 //Energy of the jet peak, pt_th > 2 GeV/c, r = 0.3
1662 //e_jet = fJetE1[0]+fJetE2[0]*e_gamma, simulation fit
1663 par[2] = fJetSigma1[1]; par[3] = fJetSigma2[1];
1664 //Sigma of the jet peak, pt_th > 2 GeV/c, r = 0.3
1665 //sigma_jet = fJetSigma1[0]+fJetSigma2[0]*e_gamma, simulation fit
1666 par[4] = fBkgMean[2 + iTPC]; par[5] = fBkgRMS[2 + iTPC];
1667 //Mean value and RMS of HIJING Bkg in a 0.3 cone, pt > 2 GeV.
1668 xmax[0] = fJetXMax1[2 + iTPC]; xmax[1] = fJetXMax2[2 + iTPC];
1669 xmin[0] = fJetXMin1[2 + iTPC]; xmin[1] = fJetXMin2[2 + iTPC];
1670 //Factor that multiplies sigma to obtain the best limits,
1671 //by observation, of mono jet ratios (ptjet/ptg) mixed with HIJING Bkg,
1672 //pt_th > 2 GeV, r = 0.3
1673 //X_jet = fJetX1[0]+fJetX2[0]*e_gamma
1674
1675 }//If low pt jet in bkg
1676 }//if Bkg
1677
1678 //Calculate minimum and maximum limits of the jet ratio.
1679 Double_t min = CalculateJetRatioLimit(ptg, par, xmin);
1680 Double_t max = CalculateJetRatioLimit(ptg, par, xmax);
1681
1682 //Info("IsJetSeleted","%s : Limits min %f, max %f, ptg / ptj %f",
1683 // type.Data(),min,max,ptj/ptg);
1684 if(( min < ptj/ptg ) && ( max > ptj/ptg))
1685 return kTRUE;
1686 else
1687 return kFALSE;
1688
1689}
1690
1691//____________________________________________________________________________
1692void AliPHOSGammaJet::List() const
1693{
1694 // List the histos
1695
1696 Info("List", "%d histograms found", fListHistos->GetEntries() ) ;
1697 TIter next(fListHistos) ;
1698 TH2F * h = 0 ;
1699 while ( (h = dynamic_cast<TH2F*>(next())) )
1700 Info("List", "%s", h->GetName()) ;
1701}
1702
1703//____________________________________________________________________________
dffafdc2 1704Double_t AliPHOSGammaJet::MakeEnergy(const Double_t energy)
1305bc01 1705{
1706 // Smears the energy according to the energy dependent energy resolution.
1707 // A gaussian distribution is assumed
1708
1709 Double_t sigma = SigmaE(energy) ;
1710 return fRan.Gaus(energy, sigma) ;
1711
1712
1713}
1714//____________________________________________________________________________
1715void AliPHOSGammaJet::MakeHistos()
1716{
1717 // Create histograms to be saved in output file and
1718 // stores them in a TObjectArray
1719
1720 fOutputFile = new TFile(fOutputFileName, "recreate") ;
1721
1722 fListHistos = new TObjArray(10000) ;
1723
1724 // Histos gamma pt vs leading pt
1725
1726 TH1F * hPtSpectra = new TH1F
1727 ("PtSpectra","p_{T i} vs p_{T #gamma}",200,0,200);
1728 hPtSpectra->SetXTitle("p_{T} (GeV/c)");
1729 fListHistos->Add(hPtSpectra) ;
1730
1731 //Histos ratio charged leading pt / gamma pt vs pt
1732
1733 TH2F * hChargeRatio = new TH2F
1734 ("ChargeRatio","p_{T leading charge} /p_{T #gamma} vs p_{T #gamma}",
1735 120,0,120,120,0,1);
1736 hChargeRatio->SetYTitle("p_{T lead charge} /p_{T #gamma}");
1737 hChargeRatio->SetXTitle("p_{T #gamma} (GeV/c)");
1738 fListHistos->Add(hChargeRatio) ;
1739
1740 TH2F * hTPCRatio = new TH2F
1741 ("TPCRatio","p_{T leading charge} /p_{T #gamma} vs p_{T #gamma}",
1742 120,0,120,120,0,1);
1743 hTPCRatio->SetYTitle("p_{T lead charge} /p_{T #gamma}");
1744 hTPCRatio->SetXTitle("p_{T #gamma} (GeV/c)");
1745 fListHistos->Add(hTPCRatio) ;
1746
1747
1748 TH2F * hPi0Ratio = new TH2F
1749 ("Pi0Ratio","p_{T leading #pi^{0}} /p_{T #gamma} vs p_{T #gamma}",
1750 120,0,120,120,0,1);
1751 hPi0Ratio->SetYTitle("p_{T lead #pi^{0}} /p_{T #gamma}");
1752 hPi0Ratio->SetXTitle("p_{T #gamma} (GeV/c)");
1753 fListHistos->Add(hPi0Ratio) ;
1754
1755 TH2F * hPhiGamma = new TH2F
1756 ("PhiGamma","#phi_{#gamma}",200,0,120,200,0,7);
1757 hPhiGamma->SetYTitle("#phi");
1758 hPhiGamma->SetXTitle("p_{T #gamma} (GeV/c)");
1759 fListHistos->Add(hPhiGamma) ;
1760
1761 TH2F * hEtaGamma = new TH2F
1762 ("EtaGamma","#phi_{#gamma}",200,0,120,200,-0.8,0.8);
1763 hEtaGamma->SetYTitle("#eta");
1764 hEtaGamma->SetXTitle("p_{T #gamma} (GeV/c)");
1765 fListHistos->Add(hEtaGamma) ;
1766
1767// //Jet reconstruction check
1768// TH2F * hDeltaPhiJet = new TH2F
1769// ("DeltaPhiJet","#phi_{jet} - #phi_{pyth jet} vs p_{T #gamma}",
1770// 200,0,120,200,-1.5,1.5);
1771// hDeltaPhiJet->SetYTitle("#Delta #phi");
1772// hDeltaPhiJet->SetXTitle("p_{T #gamma} (GeV/c)");
1773// fListHistos->Add(hDeltaPhiJet) ;
1774
1775// TH2F * hDeltaPhiTPCJet = new TH2F
1776// ("DeltaPhiTPCJet","#phi_{jet TPC} - #phi_{pyth jet} vs p_{T #gamma}",
1777// 200,0,120,200,-1.5,1.5);
1778// hDeltaPhiTPCJet->SetYTitle("#Delta #phi");
1779// hDeltaPhiTPCJet->SetXTitle("p_{T #gamma} (GeV/c)");
1780// fListHistos->Add(hDeltaPhiTPCJet) ;
1781
1782// TH2F * hDeltaEtaJet = new TH2F
1783// ("DeltaEtaJet","#phi_{jet} - #phi_{pyth jet} vs p_{T #gamma}",
1784// 200,0,120,200,-1.5,1.5);
1785// hDeltaEtaJet->SetYTitle("#Delta #phi");
1786// hDeltaEtaJet->SetXTitle("p_{T #gamma} (GeV/c)");
1787// fListHistos->Add(hDeltaEtaJet) ;
1788
1789// TH2F * hDeltaEtaTPCJet = new TH2F
1790// ("DeltaEtaTPCJet","#phi_{jet TPC} - #phi_{pyth jet} vs p_{T #gamma}",
1791// 200,0,120,200,-1.5,1.5);
1792// hDeltaEtaTPCJet->SetYTitle("#Delta #phi");
1793// hDeltaEtaTPCJet->SetXTitle("p_{T #gamma} (GeV/c)");
1794// fListHistos->Add(hDeltaEtaTPCJet) ;
1795
1796// TH2F * hDeltaPtJet = new TH2F
1797// ("DeltaPtJet","#phi_{jet} - #phi_{pyth jet} vs p_{T #gamma}",
1798// 200,0,120,200,0.,100.);
1799// hDeltaPtJet->SetYTitle("#Delta #phi");
1800// hDeltaPtJet->SetXTitle("p_{T #gamma} (GeV/c)");
1801// fListHistos->Add(hDeltaPtJet) ;
1802
1803// TH2F * hDeltaPtTPCJet = new TH2F
1804// ("DeltaPtTPCJet","#phi_{jet TPC} - #phi_{pyth jet} vs p_{T #gamma}",
1805// 200,0,120,200,0.,100.);
1806// hDeltaPtTPCJet->SetYTitle("#Delta #phi");
1807// hDeltaPtTPCJet->SetXTitle("p_{T #gamma} (GeV/c)");
1808// fListHistos->Add(hDeltaPtTPCJet) ;
1809
1810 //
1811 TH2F * hDeltaPhiCharge = new TH2F
1812 ("DeltaPhiCharge","#phi_{#gamma} - #phi_{charge} vs p_{T #gamma}",
1813 200,0,120,200,0,6.4);
1814 hDeltaPhiCharge->SetYTitle("#Delta #phi");
1815 hDeltaPhiCharge->SetXTitle("p_{T #gamma} (GeV/c)");
1816 fListHistos->Add(hDeltaPhiCharge) ;
1817
1818 TH2F * hDeltaPhiTPC = new TH2F
1819 ("DeltaPhiTPC","#phi_{#gamma} - #phi_{charge} vs p_{T #gamma}",
1820 200,0,120,200,0,6.4);
1821 hDeltaPhiTPC->SetYTitle("#Delta #phi");
1822 hDeltaPhiTPC->SetXTitle("p_{T #gamma} (GeV/c)");
1823 fListHistos->Add(hDeltaPhiTPC) ;
1824 TH2F * hDeltaPhiPi0 = new TH2F
1825 ("DeltaPhiPi0","#phi_{#gamma} - #phi_{ #pi^{0}} vs p_{T #gamma}",
1826 200,0,120,200,0,6.4);
1827 hDeltaPhiPi0->SetYTitle("#Delta #phi");
1828 hDeltaPhiPi0->SetXTitle("p_{T #gamma} (GeV/c)");
1829 fListHistos->Add(hDeltaPhiPi0) ;
1830
1831 TH2F * hDeltaEtaCharge = new TH2F
1832 ("DeltaEtaCharge","#eta_{#gamma} - #eta_{charge} vs p_{T #gamma}",
1833 200,0,120,200,-2,2);
1834 hDeltaEtaCharge->SetYTitle("#Delta #eta");
1835 hDeltaEtaCharge->SetXTitle("p_{T #gamma} (GeV/c)");
1836 fListHistos->Add(hDeltaEtaCharge) ;
1837
1838 TH2F * hDeltaEtaTPC = new TH2F
1839 ("DeltaEtaTPC","#eta_{#gamma} - #eta_{charge} vs p_{T #gamma}",
1840 200,0,120,200,-2,2);
1841 hDeltaEtaTPC->SetYTitle("#Delta #eta");
1842 hDeltaEtaTPC->SetXTitle("p_{T #gamma} (GeV/c)");
1843 fListHistos->Add(hDeltaEtaTPC) ;
1844
1845 TH2F * hDeltaEtaPi0 = new TH2F
1846 ("DeltaEtaPi0","#eta_{#gamma} - #eta_{ #pi^{0}} vs p_{T #gamma}",
1847 200,0,120,200,-2,2);
1848 hDeltaEtaPi0->SetYTitle("#Delta #eta");
1849 hDeltaEtaPi0->SetXTitle("p_{T #gamma} (GeV/c)");
1850 fListHistos->Add(hDeltaEtaPi0) ;
1851
1852 if(fOptFast){
1853
1854 TH2F * hAnglePair = new TH2F
1855 ("AnglePair",
1856 "Angle between #pi^{0} #gamma pair vs p_{T #pi^{0}}",
1857 200,0,50,200,0,0.2);
1858 hAnglePair->SetYTitle("Angle (rad)");
1859 hAnglePair->SetXTitle("E_{ #pi^{0}} (GeV/c)");
1860 fListHistos->Add(hAnglePair) ;
1861
1862 TH2F * hAnglePairAccepted = new TH2F
1863 ("AnglePairAccepted",
1864 "Angle between #pi^{0} #gamma pair vs p_{T #pi^{0}}, both #gamma in eta<0.7, inside window",
1865 200,0,50,200,0,0.2);
1866 hAnglePairAccepted->SetYTitle("Angle (rad)");
1867 hAnglePairAccepted->SetXTitle("E_{ #pi^{0}} (GeV/c)");
1868 fListHistos->Add(hAnglePairAccepted) ;
1869
1870 TH2F * hAnglePairNoCut = new TH2F
1871 ("AnglePairNoCut",
1872 "Angle between all #gamma pair vs p_{T #pi^{0}}",200,0,50,200,0,0.2);
1873 hAnglePairNoCut->SetYTitle("Angle (rad)");
1874 hAnglePairNoCut->SetXTitle("E_{ #pi^{0}} (GeV/c)");
1875 fListHistos->Add(hAnglePairNoCut) ;
1876
1877 TH2F * hAnglePairLeadingCut = new TH2F
1878 ("AnglePairLeadingCut",
1879 "Angle between all #gamma pair that have a good phi and pt vs p_{T #pi^{0}}",
1880 200,0,50,200,0,0.2);
1881 hAnglePairLeadingCut->SetYTitle("Angle (rad)");
1882 hAnglePairLeadingCut->SetXTitle("E_{ #pi^{0}} (GeV/c)");
1883 fListHistos->Add(hAnglePairLeadingCut) ;
1884
1885 TH2F * hAnglePairAngleCut = new TH2F
1886 ("AnglePairAngleCut",
1887 "Angle between all #gamma pair (angle + leading cut) vs p_{T #pi^{0}}"
1888 ,200,0,50,200,0,0.2);
1889 hAnglePairAngleCut->SetYTitle("Angle (rad)");
1890 hAnglePairAngleCut->SetXTitle("E_{ #pi^{0}} (GeV/c)");
1891 fListHistos->Add(hAnglePairAngleCut) ;
1892
1893 TH2F * hAnglePairAllCut = new TH2F
1894 ("AnglePairAllCut",
1895 "Angle between all #gamma pair (angle + inv mass cut+leading) vs p_{T #pi^{0}}"
1896 ,200,0,50,200,0,0.2);
1897 hAnglePairAllCut->SetYTitle("Angle (rad)");
1898 hAnglePairAllCut->SetXTitle("E_{ #pi^{0}} (GeV/c)");
1899 fListHistos->Add(hAnglePairAllCut) ;
1900
1901 TH2F * hAnglePairLeading = new TH2F
1902 ("AnglePairLeading",
1903 "Angle between all #gamma pair finally selected vs p_{T #pi^{0}}",
1904 200,0,50,200,0,0.2);
1905 hAnglePairLeading->SetYTitle("Angle (rad)");
1906 hAnglePairLeading->SetXTitle("E_{ #pi^{0}} (GeV/c)");
1907 fListHistos->Add(hAnglePairLeading) ;
1908
1909
1910 TH2F * hInvMassPairNoCut = new TH2F
1911 ("InvMassPairNoCut","Invariant Mass of all #gamma pair vs p_{T #gamma}",
1912 120,0,120,360,0,0.5);
1913 hInvMassPairNoCut->SetYTitle("Invariant Mass (GeV/c^{2})");
1914 hInvMassPairNoCut->SetXTitle("p_{T #gamma} (GeV/c)");
1915 fListHistos->Add(hInvMassPairNoCut) ;
1916
1917 TH2F * hInvMassPairLeadingCut = new TH2F
1918 ("InvMassPairLeadingCut",
1919 "Invariant Mass of #gamma pair (leading cuts) vs p_{T #gamma}",
1920 120,0,120,360,0,0.5);
1921 hInvMassPairLeadingCut->SetYTitle("Invariant Mass (GeV/c^{2})");
1922 hInvMassPairLeadingCut->SetXTitle("p_{T #gamma} (GeV/c)");
1923 fListHistos->Add(hInvMassPairLeadingCut) ;
1924
1925 TH2F * hInvMassPairAngleCut = new TH2F
1926 ("InvMassPairAngleCut",
1927 "Invariant Mass of #gamma pair (angle cut) vs p_{T #gamma}",
1928 120,0,120,360,0,0.5);
1929 hInvMassPairAngleCut->SetYTitle("Invariant Mass (GeV/c^{2})");
1930 hInvMassPairAngleCut->SetXTitle("p_{T #gamma} (GeV/c)");
1931 fListHistos->Add(hInvMassPairAngleCut) ;
1932
1933
1934 TH2F * hInvMassPairAllCut = new TH2F
1935 ("InvMassPairAllCut",
1936 "Invariant Mass of #gamma pair (angle+invmass cut+leading) vs p_{T #gamma}",
1937 120,0,120,360,0,0.5);
1938 hInvMassPairAllCut->SetYTitle("Invariant Mass (GeV/c^{2})");
1939 hInvMassPairAllCut->SetXTitle("p_{T #gamma} (GeV/c)");
1940 fListHistos->Add(hInvMassPairAllCut) ;
1941
1942 TH2F * hInvMassPairLeading = new TH2F
1943 ("InvMassPairLeading",
1944 "Invariant Mass of #gamma pair selected vs p_{T #gamma}",
1945 120,0,120,360,0,0.5);
1946 hInvMassPairLeading->SetYTitle("Invariant Mass (GeV/c^{2})");
1947 hInvMassPairLeading->SetXTitle("p_{T #gamma} (GeV/c)");
1948 fListHistos->Add(hInvMassPairLeading) ;
1949 }
1950
1951 //Count
1952
1953 TH1F * hNGamma = new TH1F("NGamma","Number of #gamma over PHOS",240,0,120);
1954 hNGamma->SetYTitle("N");
1955 hNGamma->SetXTitle("p_{T #gamma}(GeV/c)");
1956 fListHistos->Add(hNGamma) ;
1957
1958 TH1F * hNBkg = new TH1F("NBkg","bkg multiplicity",9000,0,9000);
1959 hNBkg->SetYTitle("counts");
1960 hNBkg->SetXTitle("N");
1961 fListHistos->Add(hNBkg) ;
1962
1963 TH2F * hNLeading = new TH2F
1964 ("NLeading","Accepted Jet Leading", 240,0,120,240,0,120);
1965 hNLeading->SetYTitle("p_{T charge} (GeV/c)");
1966 hNLeading->SetXTitle("p_{T #gamma}(GeV/c)");
1967 fListHistos->Add(hNLeading) ;
1968
1969
1970 TH1F * hN = new TH1F("NJet","Accepted jets",240,0,120);
1971 hN->SetYTitle("N");
1972 hN->SetXTitle("p_{T #gamma}(GeV/c)");
1973 fListHistos->Add(hN) ;
1974
1975
1976 //Ratios and Pt dist of reconstructed (not selected) jets
1977 //Jet
1978 TH2F * hJetRatio = new TH2F
1979 ("JetRatio","p_{T jet lead}/p_{T #gamma} vs p_{T #gamma}",
1980 240,0,120,200,0,10);
1981 hJetRatio->SetYTitle("p_{T jet lead #pi^{0}}/p_{T #gamma}");
1982 hJetRatio->SetXTitle("p_{T #gamma} (GeV/c)");
1983 fListHistos->Add(hJetRatio) ;
1984
1985
1986 TH2F * hJetPt = new TH2F
1987 ("JetPt", "p_{T jet lead} vs p_{T #gamma}",240,0,120,400,0,200);
1988 hJetPt->SetYTitle("p_{T jet lead #pi^{0}}/p_{T #gamma}");
1989 hJetPt->SetXTitle("p_{T #gamma} (GeV/c)");
1990 fListHistos->Add(hJetPt) ;
1991
1992
1993 //Bkg
1994
1995 TH2F * hBkgRatio = new TH2F
1996 ("BkgRatio","p_{T bkg lead}/p_{T #gamma} vs p_{T #gamma}",
1997 240,0,120,200,0,10);
1998 hBkgRatio->SetYTitle("p_{T bkg lead charge}/p_{T #gamma}");
1999 hBkgRatio->SetXTitle("p_{T #gamma} (GeV/c)");
2000 fListHistos->Add(hBkgRatio) ;
2001
2002
2003 TH2F * hBkgPt = new TH2F
2004 ("BkgPt","p_{T jet lead} vs p_{T #gamma}",240,0,120,400,0,200);
2005 hBkgPt->SetYTitle("p_{T jet lead charge}/p_{T #gamma}");
2006 hBkgPt->SetXTitle("p_{T #gamma} (GeV/c)");
2007 fListHistos->Add(hBkgPt) ;
2008
2009
2010 //Jet Distributions
2011
2012 TH2F * hJetFragment =
2013 new TH2F("JetFragment","x = p_{T i charged}/p_{T #gamma}",
2014 240,0.,120.,1000,0.,1.2);
2015 hJetFragment->SetYTitle("x_{T}");
2016 hJetFragment->SetXTitle("p_{T #gamma}");
2017 fListHistos->Add(hJetFragment) ;
2018
2019 TH2F * hBkgFragment = new TH2F
2020 ("BkgFragment","x = p_{T i charged}/p_{T #gamma}",
2021 240,0.,120.,1000,0.,1.2);
2022 hBkgFragment->SetYTitle("x_{T}");
2023 hBkgFragment->SetXTitle("p_{T #gamma}");
2024 fListHistos->Add(hBkgFragment) ;
2025
2026 TH2F * hJetPtDist =
2027 new TH2F("JetPtDist","x = p_{T i charged}",240,0.,120.,400,0.,200.);
2028 hJetPtDist->SetXTitle("p_{T #gamma} (GeV/c)");
2029 fListHistos->Add(hJetPtDist) ;
2030
2031 TH2F * hBkgPtDist = new TH2F
2032 ("BkgPtDist","x = p_{T i charged}",240,0.,120.,400,0.,200.);
2033 hBkgPtDist->SetXTitle("p_{T #gamma} (GeV/c)");
2034 fListHistos->Add(hBkgPtDist) ;
2035
2036
2037 if(fOnlyCharged){
2038 //Counts
2039 TH1F * hNBkgTPC = new TH1F
2040 ("NBkgTPC","TPC bkg multiplicity ",9000,0,9000);
2041 hNBkgTPC->SetYTitle("counts");
2042 hNBkgTPC->SetXTitle("N");
2043 fListHistos->Add(hNBkgTPC) ;
2044
2045 TH2F * hNTPCLeading = new TH2F
2046 ("NTPCLeading","Accepted TPC jet leading",240,0,120,240,0,120);
2047 hNTPCLeading->SetYTitle("p_{T charge} (GeV/c)");
2048 hNTPCLeading->SetXTitle("p_{T #gamma}(GeV/c)");
2049 fListHistos->Add(hNTPCLeading) ;
2050
2051 TH1F * hNTPC = new TH1F("NTPCJet","Number of TPC jets",240,0,120);
2052 hNTPC->SetYTitle("N");
2053 hNTPC->SetXTitle("p_{T #gamma}(GeV/c)");
2054 fListHistos->Add(hNTPC) ;
2055
2056 TH2F * hJetTPCRatio = new TH2F
2057 ("JetTPCRatio", "p_{T jet lead TPC}/p_{T #gamma} vs p_{T #gamma}",
2058 240,0,120,200,0,10);
2059 hJetTPCRatio->SetYTitle("p_{T jet lead TPC}/p_{T #gamma}");
2060 hJetTPCRatio->SetXTitle("p_{T #gamma} (GeV/c)");
2061 fListHistos->Add(hJetTPCRatio) ;
2062
2063 TH2F * hBkgTPCRatio = new TH2F
2064 ("BkgTPCRatio","p_{T bkg lead TPC}/p_{T #gamma} vs p_{T #gamma}",
2065 240,0,120,200,0,10);
2066 hBkgTPCRatio->SetYTitle("p_{T bkg lead TPC}/p_{T #gamma}");
2067 hBkgTPCRatio->SetXTitle("p_{T #gamma} (GeV/c)");
2068 fListHistos->Add(hBkgTPCRatio) ;
2069
2070 TH2F * hJetTPCPt = new TH2F
2071 ("JetTPCPt", "p_{T jet lead TPC} vs p_{T #gamma}",240,0,120,400,0,200);
2072 hJetTPCPt->SetYTitle("p_{T jet lead TPC}/p_{T #gamma}");
2073 hJetTPCPt->SetXTitle("p_{T #gamma} (GeV/c)");
2074 fListHistos->Add(hJetTPCPt) ;
2075
2076 TH2F * hBkgTPCPt = new TH2F
2077 ("BkgTPCPt", "p_{T bkg lead TPC} vs p_{T #gamma}",240,0,120,400,0,200);
2078 hBkgTPCPt->SetYTitle("p_{T bkg lead TPC}/p_{T #gamma}");
2079 hBkgTPCPt->SetXTitle("p_{T #gamma} (GeV/c)");
2080 fListHistos->Add(hBkgTPCPt) ;
2081
2082 //JetDistributions
2083
2084 TH2F * hJetTPCFragment =
2085 new TH2F("JetTPCFragment","x = p_{T i charged}/p_{T #gamma}",
2086 240,0.,120.,1000,0.,1.2);
2087 hJetTPCFragment->SetYTitle("x_{T}");
2088 hJetTPCFragment->SetXTitle("p_{T #gamma}");
2089 fListHistos->Add(hJetTPCFragment) ;
2090
2091 TH2F * hBkgTPCFragment = new TH2F
2092 ("BkgTPCFragment","x = p_{T i charged}/p_{T #gamma}",
2093 240,0.,120.,1000,0.,1.2);
2094 hBkgTPCFragment->SetYTitle("x_{T}");
2095 hBkgTPCFragment->SetXTitle("p_{T #gamma}");
2096 fListHistos->Add(hBkgTPCFragment) ;
2097
2098
2099 TH2F * hJetTPCPtDist = new TH2F("JetTPCPtDist",
2100 "x = p_{T i charged}",240,0.,120.,400,0.,200.);
2101 hJetTPCPtDist->SetXTitle("p_{T #gamma} (GeV/c)");
2102 fListHistos->Add(hJetTPCPtDist) ;
2103
2104 TH2F * hBkgTPCPtDist = new TH2F
2105 ("BkgTPCPtDist","x = p_{T i charged}",240,0.,120.,400,0.,200.);
2106 hBkgTPCPtDist->SetXTitle("p_{T #gamma} (GeV/c)");
2107 fListHistos->Add(hBkgTPCPtDist) ;
2108
2109 }
2110
2111
2112 if(fAnyConeOrPt){
2113 //If we want to study the jet for different cones and pt. Old version
2114
2115 TH2F * hJetRatios[5][5];
2116 TH2F * hJetTPCRatios[5][5];
2117
2118 TH2F * hJetPts[5][5];
2119 TH2F * hJetTPCPts[5][5];
2120
2121 TH2F * hBkgRatios[5][5];
2122 TH2F * hBkgTPCRatios[5][5];
2123
2124 TH2F * hBkgPts[5][5];
2125 TH2F * hBkgTPCPts[5][5];
2126
2127 TH2F * hNLeadings[5][5];
2128 TH2F * hNTPCLeadings[5][5];
2129
2130 TH1F * hNs[5][5];
2131 TH1F * hNTPCs[5][5];
2132
2133 TH1F * hNBkgs[5][5];
2134 TH1F * hNBkgTPCs[5][5];
2135
2136 TH2F * hJetFragments[5][5];
2137 TH2F * hBkgFragments[5][5];
2138 TH2F * hJetPtDists[5][5];
2139 TH2F * hBkgPtDists[5][5];
2140
2141 TH2F * hJetTPCFragments[5][5];
2142 TH2F * hBkgTPCFragments[5][5];
2143 TH2F * hJetTPCPtDists[5][5];
2144 TH2F * hBkgTPCPtDists[5][5];
2145
2146
2147 for(Int_t icone = 0; icone<fNCone; icone++){
2148 for(Int_t ipt = 0; ipt<fNPt;ipt++){
2149 //Jet Pt / Gamma Pt
2150
2151 //Jet
2152
2153 hJetRatios[icone][ipt] = new TH2F
2154 ("JetRatioCone"+fNameCones[icone]+"Pt"+fNamePtThres[ipt],
2155 "p_{T jet lead #pi^{0}}/p_{T #gamma} vs p_{T #gamma}, cone ="
2156 +fNameCones[icone]+", pt>" +fNamePtThres[ipt]+" GeV/c",
2157 240,0,120,200,0,10);
2158 hJetRatios[icone][ipt]->
2159 SetYTitle("p_{T jet lead #pi^{0}}/p_{T #gamma}");
2160 hJetRatios[icone][ipt]->SetXTitle("p_{T #gamma} (GeV/c)");
2161 fListHistos->Add(hJetRatios[icone][ipt]) ;
2162
2163
2164 hJetPts[icone][ipt] = new TH2F
2165 ("JetPtCone"+fNameCones[icone]+"Pt"+fNamePtThres[ipt],
2166 "p_{T jet lead #pi^{0}}/p_{T #gamma} vs p_{T #gamma}, cone ="
2167 +fNameCones[icone]+", pt>" +fNamePtThres[ipt]+" GeV/c",
2168 240,0,120,400,0,200);
2169 hJetPts[icone][ipt]->
2170 SetYTitle("p_{T jet lead #pi^{0}}/p_{T #gamma}");
2171 hJetPts[icone][ipt]->SetXTitle("p_{T #gamma} (GeV/c)");
2172 fListHistos->Add(hJetPts[icone][ipt]) ;
2173
2174
2175 //Bkg
2176 hBkgRatios[icone][ipt] = new TH2F
2177 ("BkgRatioCone"+fNameCones[icone]+"Pt"+fNamePtThres[ipt],
2178 "p_{T bkg lead #pi^{0}}/p_{T #gamma} vs p_{T #gamma}, cone ="
2179 +fNameCones[icone]+", pt>" +fNamePtThres[ipt]+" GeV/c",
2180 240,0,120,200,0,10);
2181 hBkgRatios[icone][ipt]->
2182 SetYTitle("p_{T bkg lead #pi^{0}}/p_{T #gamma}");
2183 hBkgRatios[icone][ipt]->SetXTitle("p_{T #gamma} (GeV/c)");
2184 fListHistos->Add(hBkgRatios[icone][ipt]) ;
2185
2186
2187
2188 hBkgPts[icone][ipt] = new TH2F
2189 ("BkgPtCone"+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 hBkgPts[icone][ipt]->
2194 SetYTitle("p_{T jet lead #pi^{0}}/p_{T #gamma}");
2195 hBkgPts[icone][ipt]->SetXTitle("p_{T #gamma} (GeV/c)");
2196 fListHistos->Add(hBkgPts[icone][ipt]) ;
2197
2198
2199 if(fOnlyCharged){
2200
2201 hJetTPCRatios[icone][ipt] = new TH2F
2202 ("JetTPCRatioCone"+fNameCones[icone]+"Pt"+fNamePtThres[ipt],
2203 "p_{T jet lead TPC}/p_{T #gamma} vs p_{T #gamma}, cone ="
2204 +fNameCones[icone]+", pt>" +fNamePtThres[ipt]+" GeV/c",
2205 240,0,120,200,0,10);
2206 hJetTPCRatios[icone][ipt]->SetYTitle("p_{T jet lead TPC}/p_{T #gamma}");
2207 hJetTPCRatios[icone][ipt]->SetXTitle("p_{T #gamma} (GeV/c)");
2208 fListHistos->Add(hJetTPCRatios[icone][ipt]) ;
2209
2210 hBkgTPCRatios[icone][ipt] = new TH2F
2211 ("BkgTPCRatioCone"+fNameCones[icone]+"Pt"+fNamePtThres[ipt],
2212 "p_{T bkg lead TPC}/p_{T #gamma} vs p_{T #gamma}, cone ="
2213 +fNameCones[icone]+", pt>" +fNamePtThres[ipt]+" GeV/c",
2214 240,0,120,200,0,10);
2215 hBkgTPCRatios[icone][ipt]->SetYTitle("p_{T bkg lead TPC}/p_{T #gamma}");
2216 hBkgTPCRatios[icone][ipt]->SetXTitle("p_{T #gamma} (GeV/c)");
2217 fListHistos->Add(hBkgTPCRatios[icone][ipt]) ;
2218
2219 hJetTPCPts[icone][ipt] = new TH2F
2220 ("JetTPCPtCone"+fNameCones[icone]+"Pt"+fNamePtThres[ipt],
2221 "p_{T jet lead TPC}/p_{T #gamma} vs p_{T #gamma}, cone ="
2222 +fNameCones[icone]+", pt>" +fNamePtThres[ipt]+" GeV/c",
2223 240,0,120,400,0,200);
2224 hJetTPCPts[icone][ipt]->SetYTitle("p_{T jet lead TPC}/p_{T #gamma}");
2225 hJetTPCPts[icone][ipt]->SetXTitle("p_{T #gamma} (GeV/c)");
2226 fListHistos->Add(hJetTPCPts[icone][ipt]) ;
2227
2228 hBkgTPCPts[icone][ipt] = new TH2F
2229 ("BkgTPCPtCone"+fNameCones[icone]+"Pt"+fNamePtThres[ipt],
2230 "p_{T bkg lead TPC}/p_{T #gamma} vs p_{T #gamma}, cone ="
2231 +fNameCones[icone]+", pt>" +fNamePtThres[ipt]+" GeV/c",
2232 240,0,120,400,0,200);
2233 hBkgTPCPts[icone][ipt]->SetYTitle("p_{T bkg lead TPC}/p_{T #gamma}");
2234 hBkgTPCPts[icone][ipt]->SetXTitle("p_{T #gamma} (GeV/c)");
2235 fListHistos->Add(hBkgTPCPts[icone][ipt]) ;
2236
2237 }
2238
2239 //Counts
2240 hNBkgs[icone][ipt] = new TH1F
2241 ("NBkgCone"+fNameCones[icone]+"Pt"+fNamePtThres[ipt],
2242 "bkg multiplicity cone ="+fNameCones[icone]+", pt>"
2243 +fNamePtThres[ipt]+" GeV/c",9000,0,9000);
2244 hNBkgs[icone][ipt]->SetYTitle("counts");
2245 hNBkgs[icone][ipt]->SetXTitle("N");
2246 fListHistos->Add(hNBkgs[icone][ipt]) ;
2247
2248 hNBkgTPCs[icone][ipt] = new TH1F
2249 ("NBkgTPCCone"+fNameCones[icone]+"Pt"+fNamePtThres[ipt],
2250 "bkg multiplicity cone ="+fNameCones[icone]+", pt>"
2251 +fNamePtThres[ipt]+" GeV/c",9000,0,9000);
2252 hNBkgTPCs[icone][ipt]->SetYTitle("counts");
2253 hNBkgTPCs[icone][ipt]->SetXTitle("N");
2254 fListHistos->Add(hNBkgTPCs[icone][ipt]) ;
2255
2256
2257 hNLeadings[icone][ipt] = new TH2F
2258 ("NLeadingCone"+fNameCones[icone]+"Pt"+fNamePtThres[ipt],
2259 "p_{T #gamma} vs p_{T #pi^{0}} cone ="+fNameCones[icone]+", pt>"
2260 +fNamePtThres[ipt]+" GeV/c",120,0,120,120,0,120);
2261 hNLeadings[icone][ipt]->SetYTitle("p_{T #pi^{0}}(GeV/c)");
2262 hNLeadings[icone][ipt]->SetXTitle("p_{T #gamma}(GeV/c)");
2263 fListHistos->Add(hNLeadings[icone][ipt]) ;
2264
2265 hNs[icone][ipt] = new TH1F
2266 ("NJetCone"+fNameCones[icone]+"Pt"+fNamePtThres[ipt],
2267 "Number of neutral jets, cone ="+fNameCones[icone]+", pt>"
2268 +fNamePtThres[ipt]+" GeV/c",120,0,120);
2269 hNs[icone][ipt]->SetYTitle("N");
2270 hNs[icone][ipt]->SetXTitle("p_{T #gamma}(GeV/c)");
2271 fListHistos->Add(hNs[icone][ipt]) ;
2272
2273 //Fragmentation Function
2274 hJetFragments[icone][ipt] = new TH2F
2275 ("JetFragmentCone"+fNameCones[icone]+"Pt"+fNamePtThres[ipt],
2276 "x_{T} = p_{T i}/p_{T #gamma}, cone ="+fNameCones[icone]+", pt>"
2277 +fNamePtThres[ipt]+" GeV/c",120,0.,120.,240,0.,1.2);
2278 hJetFragments[icone][ipt]->SetYTitle("x_{T}");
2279 hJetFragments[icone][ipt]->SetXTitle("p_{T #gamma}");
2280 fListHistos->Add(hJetFragments[icone][ipt]) ;
2281
2282 hBkgFragments[icone][ipt] = new TH2F
2283 ("BkgFragmentCone"+fNameCones[icone]+"Pt"+fNamePtThres[ipt],
2284 "x_{T} = p_{T i}/p_{T #gamma}, cone ="+fNameCones[icone]+", pt>"
2285 +fNamePtThres[ipt]+" GeV/c",120,0.,120.,240,0.,1.2);
2286 hBkgFragments[icone][ipt]->SetYTitle("x_{T}");
2287 hBkgFragments[icone][ipt]->SetXTitle("p_{T #gamma}");
2288 fListHistos->Add(hBkgFragments[icone][ipt]) ;
2289
2290
2291 //Jet particle distribution
2292
2293 hJetPtDists[icone][ipt] = new TH2F
2294 ("JetPtDistCone"+fNameCones[icone]+"Pt"+fNamePtThres[ipt],
2295 "p_{T i}, cone ="+fNameCones[icone]+", pt>" +fNamePtThres[ipt]+
2296 " GeV/c",120,0.,120.,120,0.,120.);
2297 hJetPtDists[icone][ipt]->SetXTitle("p_{T #gamma} (GeV/c)");
2298 fListHistos->Add(hJetPtDists[icone][ipt]) ;
2299
2300 hBkgPtDists[icone][ipt] = new TH2F
2301 ("BkgPtDistCone"+fNameCones[icone]+"Pt"+fNamePtThres[ipt],
2302 "p_{T i}, cone ="+fNameCones[icone]+", pt>" +fNamePtThres[ipt]+
2303 " GeV/c",120,0.,120.,120,0.,120.);
2304 hBkgPtDists[icone][ipt]->SetXTitle("p_{T #gamma} (GeV/c)");
2305 fListHistos->Add(hBkgPtDists[icone][ipt]) ;
2306
2307
2308 if(fOnlyCharged){
2309 hNTPCLeadings[icone][ipt] = new TH2F
2310 ("NTPCLeadingCone"+fNameCones[icone]+"Pt"+fNamePtThres[ipt],
2311 "p_{T #gamma} vs p_{T charge} cone ="+fNameCones[icone]+", pt>"
2312 +fNamePtThres[ipt]+" GeV/c",120,0,120,120,0,120);
2313 hNTPCLeadings[icone][ipt]->SetYTitle("p_{T charge} (GeV/c)");
2314 hNTPCLeadings[icone][ipt]->SetXTitle("p_{T #gamma}(GeV/c)");
2315 fListHistos->Add(hNTPCLeadings[icone][ipt]) ;
2316
2317 hNTPCs[icone][ipt] = new TH1F
2318 ("NTPCJetCone"+fNameCones[icone]+"Pt"+fNamePtThres[ipt],
2319 "Number of charged jets, cone ="+fNameCones[icone]+", pt>"
2320 +fNamePtThres[ipt]+" GeV/c",120,0,120);
2321 hNTPCs[icone][ipt]->SetYTitle("N");
2322 hNTPCs[icone][ipt]->SetXTitle("p_{T #gamma}(GeV/c)");
2323 fListHistos->Add(hNTPCs[icone][ipt]) ;
2324
2325 hJetTPCFragments[icone][ipt] = new TH2F
2326 ("JetTPCFragmentCone"+fNameCones[icone]+"Pt"+fNamePtThres[ipt],
2327 "x = p_{T i charged}/p_{T #gamma}, cone ="+fNameCones[icone]+", pt>"
2328 +fNamePtThres[ipt]+" GeV/c",120,0.,120.,240,0.,1.2);
2329 hJetTPCFragments[icone][ipt]->SetYTitle("x_{T}");
2330 hJetTPCFragments[icone][ipt]->SetXTitle("p_{T #gamma}");
2331 fListHistos->Add(hJetTPCFragments[icone][ipt]) ;
2332
2333 hBkgTPCFragments[icone][ipt] = new TH2F
2334 ("BkgTPCFragmentCone"+fNameCones[icone]+"Pt"+fNamePtThres[ipt],
2335 "x = p_{T i charged}/p_{T #gamma}, cone ="+fNameCones[icone]+", pt>"
2336 +fNamePtThres[ipt]+" GeV/c",120,0.,120.,240,0.,1.2);
2337 hBkgTPCFragments[icone][ipt]->SetYTitle("x_{T}");
2338 hBkgTPCFragments[icone][ipt]->SetXTitle("p_{T #gamma}");
2339 fListHistos->Add(hBkgTPCFragments[icone][ipt]) ;
2340
2341 hJetTPCPtDists[icone][ipt] = new TH2F
2342 ("JetTPCPtDistCone"+fNameCones[icone]+"Pt"+fNamePtThres[ipt],
2343 "x = p_{T i charged}, cone ="+fNameCones[icone]+", pt>"
2344 +fNamePtThres[ipt]+" GeV/c",120,0.,120.,120,0.,120.);
2345 hJetTPCPtDists[icone][ipt]->SetXTitle("p_{T #gamma} (GeV/c)");
2346 fListHistos->Add(hJetTPCPtDists[icone][ipt]) ;
2347
2348 hBkgTPCPtDists[icone][ipt] = new TH2F
2349 ("BkgTPCPtDistCone"+fNameCones[icone]+"Pt"+fNamePtThres[ipt],
2350 "x = p_{T i charged}, cone ="+fNameCones[icone]+", pt>" +
2351 fNamePtThres[ipt]+" GeV/c",120,0.,120.,120,0.,120.);
2352 hBkgTPCPtDists[icone][ipt]->SetXTitle("p_{T #gamma} (GeV/c)");
2353 fListHistos->Add(hBkgTPCPtDists[icone][ipt]) ;
2354 }
2355 }//ipt
2356 } //icone
2357 }//If we want to study any cone or pt threshold
2358}
2359
2360
2361//____________________________________________________________________________
2362void AliPHOSGammaJet::MakeJet(TClonesArray * pl,
2363 Double_t ptg, Double_t phig,
2364 Double_t ptl, Double_t phil, Double_t etal,
2365 TString conf, TLorentzVector & jet)
2366{
e8daeb92 2367 //Fill the jet with the particles around the leading particle with
2368 //R=fCone and pt_th = fPtThres. Calculate the energy of the jet and
2369 //check if we select it. Fill jet histograms
1305bc01 2370 Float_t ptcut = 0. ;
2371 if(fHIJING){
2372 if(ptg > fPtJetSelectionCut) ptcut = 2. ;
2373 else ptcut = 0.5;
2374 }
2375
2376 TClonesArray * jetList = new TClonesArray("TParticle",1000);
2377 TClonesArray * bkgList = new TClonesArray("TParticle",1000);
2378
2379 TLorentzVector bkg(0,0,0,0);
2380 TLorentzVector lv (0,0,0,0);
2381
2382 Double_t ptjet = 0.0;
2383 Double_t ptbkg = 0.0;
2384 Int_t n0 = 0;
2385 Int_t n1 = 0;
2386 Bool_t b1 = kFALSE;
2387 Bool_t b0 = kFALSE;
2388
2389 TIter next(pl) ;
2390 TParticle * particle = 0 ;
2391
2392 while ( (particle = dynamic_cast<TParticle*>(next())) ) {
2393
2394 b0 = kFALSE;
2395 b1 = kFALSE;
2396
2397 //Particles in jet
2398 SetJet(particle, b0, fCone, etal, phil) ;
2399
2400 if(b0){
2401 new((*jetList)[n0++]) TParticle(*particle) ;
2402 particle->Momentum(lv);
2403 if(particle->Pt() > ptcut ){
2404 jet+=lv;
2405 ptjet+=particle->Pt();
2406 }
2407 }
2408
2409 //Background around (phi_gamma-pi, eta_leading)
2410 SetJet(particle, b1, fCone,etal, phig) ;
2411
2412 if(b1) {
2413 new((*bkgList)[n1++]) TParticle(*particle) ;
2414 if(particle->Pt() > ptcut ){
2415 bkg+=lv;
2416 ptbkg+=particle->Pt();
2417 }
2418 }
2419 }
2420
2421 ptjet = jet.Pt();
2422 ptbkg = bkg.Pt();
2423
e8daeb92 2424 if(strstr(fOptionGJ,"deb") || strstr(fOptionGJ,"deb all"))
1305bc01 2425 Info("MakeJet","Gamma pt %f, Jet pt %f, Bkg pt %f",ptg,ptjet,ptbkg);
2426
2427
2428 //Fill histograms
2429
2430 Double_t rat = ptjet/ptg ;
2431 Double_t ratbkg = ptbkg/ptg ;
2432
2433 dynamic_cast<TH2F*>
2434 (fListHistos->FindObject("Jet"+conf+"Ratio"))->Fill(ptg,rat);
2435 dynamic_cast<TH2F*>
2436 (fListHistos->FindObject("Jet"+conf+"Pt")) ->Fill(ptg,ptjet);
2437 dynamic_cast<TH2F*>
2438 (fListHistos->FindObject("Bkg"+conf+"Ratio"))->Fill(ptg,ratbkg);
2439 dynamic_cast<TH2F*>
2440 (fListHistos->FindObject("Bkg"+conf+"Pt")) ->Fill(ptg,ptbkg);
2441
2442
2443 if(IsJetSelected(ptg,ptjet,conf) || fSelect){
e8daeb92 2444 if(strstr(fOptionGJ,"deb") || strstr(fOptionGJ,"deb all"))
1305bc01 2445 Info("MakeJet","JetSelected");
2446 dynamic_cast<TH1F*>(fListHistos->FindObject("N"+conf+"Jet"))->
2447 Fill(ptg);
2448 dynamic_cast<TH2F*>(fListHistos->FindObject("N"+conf+"Leading"))
2449 ->Fill(ptg,ptl);
2450 FillJetHistos(jetList, ptg, conf, "Jet");
2451 FillJetHistos(bkgList, ptg, conf, "Bkg");
2452 }
2453
2454 jetList ->Delete();
2455 bkgList ->Delete();
2456}
2457
2458//____________________________________________________________________________
2459void AliPHOSGammaJet::MakeJetAnyConeOrPt(TClonesArray * pl, Double_t ptg,
2460 Double_t phig, Double_t ptl,
2461 Double_t phil, Double_t etal,
2462 TString conf){
e8daeb92 2463
2464 //Fill the jet with the particles around the leading particle with
2465 //R=fCone(i) and pt_th = fPtThres(i). Calculate the energy of the jet and
2466 //check if we select it. Fill jet i histograms
1305bc01 2467
2468 TClonesArray * jetList = new TClonesArray("TParticle",1000);
2469 TClonesArray * bkgList = new TClonesArray("TParticle",1000);
2470
2471 Double_t ptjet = 0.0;
2472 Double_t ptbkg = 0.0;
2473
2474 Int_t n1 = 0;
2475 Int_t n0 = 0;
2476 Bool_t b1 = kFALSE;
2477 Bool_t b0 = kFALSE;
2478
2479 //Create as many jets as cones and pt thresholds are defined
2480 Double_t maxcut = fJetRatioMaxCut;
2481 Double_t mincut = fJetRatioMinCut;
2482
2483 if(conf == "TPC"){
2484 maxcut = fJetTPCRatioMaxCut;
2485 mincut = fJetTPCRatioMinCut;
2486 }
2487
2488 Double_t ratjet = 0;
2489 Double_t ratbkg = 0;
2490
2491 for(Int_t icone = 0; icone<fNCone; icone++) {
2492 for(Int_t ipt = 0; ipt<fNPt;ipt++) {
2493
2494 TString cone = fNameCones[icone] ;
2495 TString ptcut = fNamePtThres[ipt] ;
2496
2497 TIter next(pl) ;
2498 TParticle * particle = 0 ;
2499
2500 ptjet = 0 ;
2501 ptbkg = 0 ;
2502
2503 while ( (particle = dynamic_cast<TParticle*>(next())) ) {
2504 b1 = kFALSE;
2505 b0 = kFALSE;
2506
2507 SetJet(particle, b0, fCones[icone],etal, phil) ;
2508 SetJet(particle, b1, fCones[icone],etal, phig) ;
2509
2510 if(b0){
2511 new((*jetList)[n0++]) TParticle(*particle) ;
2512 if(particle->Pt() > fPtThres[ipt] )
2513 ptjet+=particle->Pt();
2514 }
2515 if(b1) {
2516 new((*bkgList)[n1++]) TParticle(*particle) ;
2517 if(particle->Pt() > fPtThres[ipt] )
2518 ptbkg+=particle->Pt();
2519 }
2520
2521 }
2522
2523 //Fill histograms
2524 if(ptjet > 0.) {
2525
e8daeb92 2526 if(strstr(fOptionGJ,"deb")){
1305bc01 2527 Info("MakeJetAnyPt","cone %f, ptcut %f",fCones[icone],fPtThres[ipt]);
2528 Info("MakeJetAnyPt","pT: Gamma %f, Jet %f, Bkg %f",ptg,ptjet,ptbkg);
2529 }
2530
2531 ratjet = ptjet /ptg;
2532 ratbkg = ptbkg/ptg;
2533
2534 dynamic_cast<TH2F*>
2535 (fListHistos->FindObject("Jet"+conf+"RatioCone"+cone+"Pt"+ptcut))
2536 ->Fill(ptg,ratjet);
2537 dynamic_cast<TH2F*>
2538 (fListHistos->FindObject("Jet"+conf+"PtCone"+cone+"Pt"+ptcut))
2539 ->Fill(ptg,ptjet);
2540
2541 dynamic_cast<TH2F*>
2542 (fListHistos->FindObject("Bkg"+conf+"RatioCone"+cone+"Pt"+ptcut))
2543 ->Fill(ptg,ratbkg);
2544
2545 dynamic_cast<TH2F*>
2546 (fListHistos->FindObject("Bkg"+conf+"PtCone"+cone+"Pt"+ptcut))
2547 ->Fill(ptg,ptbkg);
2548
2549
2550 //Select Jet
2551 if((ratjet < maxcut) && (ratjet > mincut)){
2552
2553 dynamic_cast<TH1F*>
2554 (fListHistos->FindObject("N"+conf+"JetCone"+cone+"Pt"+ptcut))->
2555 Fill(ptg);
2556 dynamic_cast<TH2F*>
2557 (fListHistos->FindObject("N"+conf+"LeadingCone"+cone+"Pt"+ptcut))
2558 ->Fill(ptg,ptl);
2559
2560 FillJetHistosAnyConeOrPt
2561 (jetList,ptg,conf,"Jet",fNameCones[icone],fNamePtThres[ipt]);
2562 FillJetHistosAnyConeOrPt
2563 (bkgList,ptg,conf,"Bkg",fNameCones[icone],fNamePtThres[ipt]);
2564
2565 }
2566 } //ptjet > 0
2567 jetList ->Delete();
2568 bkgList ->Delete();
2569 }//for pt threshold
2570 }// for cone
2571}
2572
2573//____________________________________________________________________________
2574void AliPHOSGammaJet::MakePhoton(TLorentzVector & particle)
2575{
e8daeb92 2576 //Fast reconstruction for photons
1305bc01 2577 Double_t energy = particle.E() ;
2578 Double_t modenergy = MakeEnergy(energy) ;
2579 //Info("MakePhoton","Energy %f, Modif %f",energy,modenergy);
2580
2581 // get the detected direction
2582 TVector3 pos = particle.Vect();
2583 pos*=460./energy;
2584 TVector3 modpos = MakePosition(energy, pos) ;
2585 modpos *= modenergy / 460.;
2586
2587 Float_t modtheta = modpos.Theta();
2588 Float_t modphi = modpos.Phi();
2589
2590 // Set the modified 4-momentum of the reconstructed particle
2591 Float_t py = modenergy*TMath::Sin(modphi)*TMath::Sin(modtheta);
2592 Float_t px = modenergy*TMath::Cos(modphi)*TMath::Sin(modtheta);
2593 Float_t pz = modenergy*TMath::Cos(modtheta);
2594
2595 particle.SetPxPyPzE(px,py,pz,modenergy);
2596
2597}
2598
2599//____________________________________________________________________________
dffafdc2 2600TVector3 AliPHOSGammaJet::MakePosition(const Double_t energy, const TVector3 pos)
1305bc01 2601{
2602 // Smears the impact position according to the energy dependent position resolution
2603 // A gaussian position distribution is assumed
2604
2605 TVector3 newpos ;
2606
2607 Double_t sigma = SigmaP(energy) ;
2608 Double_t x = fRan.Gaus( pos.X(), sigma ) ;
2609 Double_t z = fRan.Gaus( pos.Z(), sigma ) ;
2610 Double_t y = pos.Y() ;
2611
2612 newpos.SetX(x) ;
2613 newpos.SetY(y) ;
2614 newpos.SetZ(z) ;
2615
2616 // Info("MakePosition","Theta dif %f",pos.Theta()-newpos.Theta());
2617// Info("MakePosition","Phi dif %f",pos.Phi()-newpos.Phi());
2618 return newpos ;
2619}
2620
77414c90 2621//____________________________________________________________________________
2622void AliPHOSGammaJet::Pi0Decay(Double_t mPi0, TLorentzVector &p0,
2623 TLorentzVector &p1, TLorentzVector &p2, Double_t &angle) {
2624 // Perform isotropic decay pi0 -> 2 photons
1305bc01 2625 // p0 is pi0 4-momentum (inut)
77414c90 2626 // p1 and p2 are photon 4-momenta (output)
2627 // cout<<"Boost vector"<<endl;
2628 TVector3 b = p0.BoostVector();
2629 //cout<<"Parameters"<<endl;
2630 //Double_t mPi0 = p0.M();
2631 Double_t phi = TMath::TwoPi() * gRandom->Rndm();
2632 Double_t cosThe = 2 * gRandom->Rndm() - 1;
2633 Double_t cosPhi = TMath::Cos(phi);
2634 Double_t sinPhi = TMath::Sin(phi);
2635 Double_t sinThe = TMath::Sqrt(1-cosThe*cosThe);
2636 Double_t ePi0 = mPi0/2.;
2637 //cout<<"ePi0 "<<ePi0<<endl;
2638 //cout<<"Components"<<endl;
2639 p1.SetPx(+ePi0*cosPhi*sinThe);
2640 p1.SetPy(+ePi0*sinPhi*sinThe);
2641 p1.SetPz(+ePi0*cosThe);
2642 p1.SetE(ePi0);
2643 //cout<<"p1: "<<p1.Px()<<" "<<p1.Py()<<" "<<p1.Pz()<<" "<<p1.E()<<endl;
2644 //cout<<"p1 Mass: "<<p1.Px()*p1.Px()+p1.Py()*p1.Py()+p1.Pz()*p1.Pz()-p1.E()*p1.E()<<endl;
2645 p2.SetPx(-ePi0*cosPhi*sinThe);
2646 p2.SetPy(-ePi0*sinPhi*sinThe);
2647 p2.SetPz(-ePi0*cosThe);
2648 p2.SetE(ePi0);
2649 //cout<<"p2: "<<p2.Px()<<" "<<p2.Py()<<" "<<p2.Pz()<<" "<<p2.E()<<endl;
2650 //cout<<"p2 Mass: "<<p2.Px()*p2.Px()+p2.Py()*p2.Py()+p2.Pz()*p2.Pz()-p2.E()*p2.E()<<endl;
2651 //cout<<"Boost "<<b.X()<<" "<<b.Y()<<" "<<b.Z()<<endl;
2652 p1.Boost(b);
2653 //cout<<"p1: "<<p1.Px()<<" "<<p1.Py()<<" "<<p1.Pz()<<" "<<p1.E()<<endl;
2654 p2.Boost(b);
2655 //cout<<"p2: "<<p2.Px()<<" "<<p2.Py()<<" "<<p2.Pz()<<" "<<p2.E()<<endl;
2656 //cout<<"angle"<<endl;
2657 angle = p1.Angle(p2.Vect());
2658 //cout<<angle<<endl;
2659}
1305bc01 2660
2661//____________________________________________________________________________
2662void AliPHOSGammaJet::Plot(TString what, Option_t * option) const
2663{
e8daeb92 2664 //Plot some relevant histograms of the analysis
1305bc01 2665 TH2F * h = dynamic_cast<TH2F*>(fOutputFile->Get(what));
2666 if(h){
2667 h->Draw();
2668 }
2669 else if (what == "all") {
2670 TCanvas * can = new TCanvas("GammaJet", "Gamma-Jet Study",10,40,1600,1200);
2671 can->cd() ;
2672 can->Range(0,0,22,20);
2673 TPaveLabel *pl1 = new TPaveLabel(1,18,20,19.5,"Titre","br");
2674 pl1->SetFillColor(18);
2675 pl1->SetTextFont(32);
2676 pl1->SetTextColor(49);
2677 pl1->Draw();
2678 Int_t index ;
2679 TPad * pad = 0;
2680 Float_t begx = -0.29, begy = 0., endx = 0., endy = 0.30 ;
2681 for (index = 0 ; index < fListHistos->GetEntries() ; index++) {
2682 TString name("pad") ;
2683 name += index ;
2684 begx += 0.30 ;
2685 endx += 0.30 ;
2686 if (begx >= 1.0 || endx >= 1.0) {
2687 begx = 0.01 ;
2688 endx = 0.30 ;
2689 begy += 0.30 ;
2690 endy += 0.30 ;
2691 }
2692 printf("%f %f %f %f \n", begx, begy, endx, endy) ;
2693 pad = new TPad(name,"This is a pad",begx,begy,endx,endy,33);
2694 pad->Draw();
2695 pad->cd() ;
2696 fListHistos->At(index)->Draw(option) ;
2697 pad->Modified() ;
2698 pad->Update() ;
2699 can->cd() ;
2700 }
2701
2702 }
2703 else{
2704 Info("Draw", "Histogram %s does not exist or unknown option", what.Data()) ;
2705 fOutputFile->ls();
2706 }
2707}
2708
77414c90 2709//____________________________________________________________________________
702ab87e 2710void AliPHOSGammaJet::Print(const Option_t * opt) const
1305bc01 2711{
e8daeb92 2712
2713 //Print some relevant parameters set for the analysis
1305bc01 2714 if(! opt)
2715 return;
2716
1305bc01 2717 Info("Print", "%s %s", GetName(), GetTitle() ) ;
2718
2719 printf("Eta cut : %f\n", fEtaCut) ;
2720 printf("D phi max cut : %f\n", fPhiMaxCut) ;
2721 printf("D phi min cut : %f\n", fPhiMinCut) ;
2722 printf("Leading Ratio max cut : %f\n", fRatioMaxCut) ;
2723 printf("Leading Ratio min cut : %f\n", fRatioMinCut) ;
2724 printf("Jet Ratio max cut : %f\n", fJetRatioMaxCut) ;
2725 printf("Jet Ratio min cut : %f\n", fJetRatioMinCut) ;
2726 printf("Jet TPC Ratio max cut : %f\n", fJetTPCRatioMaxCut) ;
2727 printf("Jet TPC Ratio min cut : %f\n", fJetTPCRatioMinCut) ;
2728 printf("Fast recons : %d\n", fOptFast);
2729 printf("Inv Mass max cut : %f\n", fInvMassMaxCut) ;
2730 printf("Inv Mass min cut : %f\n", fInvMassMinCut) ;
2731
2732}
2733
2734//___________________________________________________________________
2735void AliPHOSGammaJet::SetJet(TParticle * part, Bool_t & b, Float_t cone,
2736 Double_t eta, Double_t phi)
77414c90 2737{
e8daeb92 2738
2739 //Check if the particle is inside the cone defined by the leading particle
1305bc01 2740 b = kFALSE;
2741
2742 if(phi > TMath::TwoPi())
2743 phi-=TMath::TwoPi();
2744 if(phi < 0.)
2745 phi+=TMath::TwoPi();
2746
2747 Double_t rad = 10000 + cone;
2748
2749 if(TMath::Abs(part->Phi()-phi) <= (TMath::TwoPi() - cone))
2750 rad = TMath::Sqrt(TMath::Power(part->Eta()-eta,2)+
2751 TMath::Power(part->Phi()-phi,2));
2752 else{
2753 if(part->Phi()-phi > TMath::TwoPi() - cone)
2754 rad = TMath::Sqrt(TMath::Power(part->Eta()-eta,2)+
2755 TMath::Power((part->Phi()-TMath::TwoPi())-phi,2));
2756 if(part->Phi()-phi < -(TMath::TwoPi() - cone))
2757 rad = TMath::Sqrt(TMath::Power(part->Eta()-eta,2)+
2758 TMath::Power((part->Phi()+TMath::TwoPi())-phi,2));
2759 }
2760
2761 if(rad < cone )
2762 b = kTRUE;
77414c90 2763
77414c90 2764}
1305bc01 2765
77414c90 2766//____________________________________________________________________________
1305bc01 2767Double_t AliPHOSGammaJet::SigmaE(Double_t energy)
77414c90 2768{
1305bc01 2769 // Calculates the energy dependent energy resolution
2770
2771 Double_t rv = -1 ;
2772
2773 rv = TMath::Sqrt( TMath::Power(fResPara1/energy, 2)
2774 + TMath::Power(fResPara2/TMath::Sqrt(energy), 2)
2775 + TMath::Power(fResPara3, 2) ) ;
2776
2777 return rv * energy ;
77414c90 2778}
2779
2780//____________________________________________________________________________
1305bc01 2781Double_t AliPHOSGammaJet::SigmaP(Double_t energy)
77414c90 2782{
1305bc01 2783 // Calculates the energy dependent position resolution
2784
2785 Double_t sigma = TMath::Sqrt(TMath::Power(fPosParaA,2) +
2786 TMath::Power(fPosParaB,2) / energy) ;
2787
2788
2789 return sigma ; // in cm
77414c90 2790}
2791