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