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