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