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