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