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