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