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