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