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