New Gamma package
[u/mrichter/AliRoot.git] / PWG4 / AliAnaGammaJet.cxx
CommitLineData
f9cea31c 1/**************************************************************************
2 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
3 * *
4 * Author: The ALICE Off-line Project. *
5 * Contributors are mentioned in the code where appropriate. *
6 * *
7 * Permission to use, copy, modify and distribute this software and its *
8 * documentation strictly for non-commercial purposes is hereby granted *
9 * without fee, provided that the above copyright notice appears in all *
10 * copies and that both the copyright notice and this permission notice *
11 * appear in the supporting documentation. The authors make no claims *
12 * about the suitability of this software for any purpose. It is *
13 * provided "as is" without express or implied warranty. *
14 **************************************************************************/
15/* $Id$ */
16
17/* History of cvs commits:
18 *
19 * $Log$
20 *
21 */
22
23//_________________________________________________________________________
24// Class for the analysis of gamma-jet correlations
25// Basically it seaches for a prompt photon in the Calorimeters acceptance,
26// if so we construct a jet around the highest pt particle in the opposite
27// side in azimuth, inside the Central Tracking System (ITS+TPC) and
28// EMCAL acceptances (only when PHOS detects the gamma). First the leading
29// particle and then the jet have to fullfill several conditions
30// (energy, direction ..) to be accepted. Then the fragmentation function
31// of this jet is constructed
32// Class created from old AliPHOSGammaJet
33// (see AliRoot versions previous Release 4-09)
34//
35//*-- Author: Gustavo Conesa (LNF-INFN)
36//////////////////////////////////////////////////////////////////////////////
37
38
39// --- ROOT system ---
40
41#include <TFile.h>
42#include <TParticle.h>
43#include <TH2.h>
44
45#include "AliAnaGammaJet.h"
46#include "AliESD.h"
47#include "AliESDtrack.h"
48#include "AliESDCaloCluster.h"
49#include "Riostream.h"
50#include "AliLog.h"
51
52ClassImp(AliAnaGammaJet)
53
54//____________________________________________________________________________
55AliAnaGammaJet::AliAnaGammaJet(const char *name) :
56 AliAnaGammaDirect(name),
57 fSeveralConeAndPtCuts(0),
58 fPbPb(kFALSE),
59 fJetsOnlyInCTS(0),
60 fEtaEMCALCut(0.),fPhiMaxCut(0.),
61 fPhiMinCut(0.),
62 fInvMassMaxCut(0.), fInvMassMinCut(0.),
63 fJetCTSRatioMaxCut(0.),
64 fJetCTSRatioMinCut(0.), fJetRatioMaxCut(0.),
65 fJetRatioMinCut(0.), fNCone(0),
66 fNPt(0), fCone(0), fPtThreshold(0),
67 fPtJetSelectionCut(0.0),
68 fAngleMaxParam(), fSelect(0)
69{
70
71 // ctor
72 fAngleMaxParam.Set(4) ;
73 fAngleMaxParam.Reset(0.);
74
75 for(Int_t i = 0; i<10; i++){
76 fCones[i] = 0.0 ;
77 fNameCones[i] = "" ;
78 fPtThres[i] = 0.0 ;
79 fNamePtThres[i] = "" ;
80 if( i < 6 ){
81 fJetXMin1[i] = 0.0 ;
82 fJetXMin2[i] = 0.0 ;
83 fJetXMax1[i] = 0.0 ;
84 fJetXMax2[i] = 0.0 ;
85 fBkgMean[i] = 0.0 ;
86 fBkgRMS[i] = 0.0 ;
87 if( i < 2 ){
88 fJetE1[i] = 0.0 ;
89 fJetE2[i] = 0.0 ;
90 fJetSigma1[i] = 0.0 ;
91 fJetSigma2[i] = 0.0 ;
92 fPhiEMCALCut[i] = 0.0 ;
93 }
94 }
95 }
96
97 TList * list = gDirectory->GetListOfKeys() ;
98 TIter next(list) ;
99 TH2F * h = 0 ;
100 Int_t index ;
101 for (index = 0 ; index < list->GetSize()-1 ; index++) {
102 //-1 to avoid GammaJet Task
103 h = dynamic_cast<TH2F*>(gDirectory->Get(list->At(index)->GetName())) ;
104 fOutputContainer->Add(h) ;
105 }
106
107 // Input slot #0 works with an Ntuple
108 DefineInput(0, TChain::Class());
109 // Output slot #0 writes into a TH1 container
110 DefineOutput(0, TObjArray::Class()) ;
111
112}
113
114
115//____________________________________________________________________________
116AliAnaGammaJet::AliAnaGammaJet(const AliAnaGammaJet & gj) :
117 AliAnaGammaDirect(gj),
118 fSeveralConeAndPtCuts(gj.fSeveralConeAndPtCuts),
119 fPbPb(gj.fPbPb), fJetsOnlyInCTS(gj.fJetsOnlyInCTS),
120 fEtaEMCALCut(gj.fEtaEMCALCut),
121 fPhiMaxCut(gj.fPhiMaxCut), fPhiMinCut(gj.fPhiMinCut),
122 fInvMassMaxCut(gj.fInvMassMaxCut), fInvMassMinCut(gj.fInvMassMinCut),
123 fRatioMinCut(gj.fRatioMinCut),
124 fJetCTSRatioMaxCut(gj.fJetCTSRatioMaxCut),
125 fJetCTSRatioMinCut(gj.fJetCTSRatioMinCut), fJetRatioMaxCut(gj.fJetRatioMaxCut),
126 fJetRatioMinCut(gj.fJetRatioMinCut), fNCone(gj.fNCone),
127 fNPt(gj.fNPt), fCone(gj.fCone), fPtThreshold(gj.fPtThreshold),
128 fPtJetSelectionCut(gj.fPtJetSelectionCut),
129 fOutputContainer(0), fAngleMaxParam(gj.fAngleMaxParam),
130 fSelect(gj.fSelect)
131{
132 // cpy ctor
133 SetName (gj.GetName()) ;
134 SetTitle(gj.GetTitle()) ;
135
136 for(Int_t i = 0; i<10; i++){
137 fCones[i] = gj.fCones[i] ;
138 fNameCones[i] = gj.fNameCones[i] ;
139 fPtThres[i] = gj.fPtThres[i] ;
140 fNamePtThres[i] = gj.fNamePtThres[i] ;
141 if( i < 6 ){
142 fJetXMin1[i] = gj.fJetXMin1[i] ;
143 fJetXMin2[i] = gj.fJetXMin2[i] ;
144 fJetXMax1[i] = gj.fJetXMax1[i] ;
145 fJetXMax2[i] = gj.fJetXMax2[i] ;
146 fBkgMean[i] = gj.fBkgMean[i] ;
147 fBkgRMS[i] = gj.fBkgRMS[i] ;
148 if( i < 2 ){
149 fJetE1[i] = gj.fJetE1[i] ;
150 fJetE2[i] = gj.fJetE2[i] ;
151 fJetSigma1[i] = gj.fJetSigma1[i] ;
152 fJetSigma2[i] = gj.fJetSigma2[i] ;
153 fPhiEMCALCut[i] = gj.fPhiEMCALCut[i] ;
154 }
155 }
156 }
157}
158
159//____________________________________________________________________________
160AliAnaGammaJet::~AliAnaGammaJet()
161{
162 // Remove all pointers
163 fOutputContainer->Clear() ;
164 delete fOutputContainer ;
165
166 delete fhChargeRatio ;
167 delete fhPi0Ratio ;
168 delete fhDeltaPhiCharge ;
169 delete fhDeltaPhiPi0 ;
170 delete fhDeltaEtaCharge ;
171 delete fhDeltaEtaPi0 ;
172 delete fhAnglePair ;
173 delete fhAnglePairAccepted ;
174 delete fhAnglePairNoCut ;
175 delete fhAnglePairLeadingCut ;
176 delete fhAnglePairAngleCut ;
177 delete fhAnglePairAllCut ;
178 delete fhAnglePairLeading ;
179 delete fhInvMassPairNoCut ;
180 delete fhInvMassPairLeadingCut ;
181 delete fhInvMassPairAngleCut ;
182 delete fhInvMassPairAllCut ;
183 delete fhInvMassPairLeading ;
184 delete fhNBkg ;
185 delete fhNLeading ;
186 delete fhNJet ;
187 delete fhJetRatio ;
188 delete fhJetPt ;
189 delete fhBkgRatio ;
190 delete fhBkgPt ;
191 delete fhJetFragment ;
192 delete fhBkgFragment ;
193 delete fhJetPtDist ;
194 delete fhBkgPtDist ;
195
196 delete [] fhJetRatios;
197 delete [] fhJetPts;
198 delete [] fhBkgRatios;
199 delete [] fhBkgPts;
200
201 delete [] fhNLeadings;
202 delete [] fhNJets;
203 delete [] fhNBkgs;
204
205 delete [] fhJetFragments;
206 delete [] fhBkgFragments;
207 delete [] fhJetPtDists;
208 delete [] fhBkgPtDists;
209
210}
211
212//____________________________________________________________________________
213Double_t AliAnaGammaJet::CalculateJetRatioLimit(const Double_t ptg,
214 const Double_t *par,
215 const Double_t *x) {
216
217 //Parametrized cut for the energy of the jet.
218
219 Double_t epp = par[0] + par[1] * ptg ;
220 Double_t spp = par[2] + par[3] * ptg ;
221 Double_t f = x[0] + x[1] * ptg ;
222 Double_t epb = epp + par[4] ;
223 Double_t spb = TMath::Sqrt(spp*spp+ par[5]*par[5]) ;
224 Double_t rat = (epb - spb * f) / ptg ;
225 //Info("CalculateLimit","epp %f, spp %f, f %f", epp, spp, f);
226 //Info("CalculateLimit","epb %f, spb %f, rat %f", epb, spb, rat);
227 return rat ;
228}
229
230
231//____________________________________________________________________________
232void AliAnaGammaJet::Exec(Option_t *)
233{
234
235 // Processing of one event
236
237 //Get ESDs
238 Long64_t entry = GetChain()->GetReadEntry() ;
239
240 if (!GetESD()) {
241 AliError("fESD is not connected to the input!") ;
242 return ;
243 }
244
245 if (GetPrintInfo())
246 AliInfo(Form("%s ----> Processing event # %lld", (dynamic_cast<TChain *>(GetChain()))->GetFile()->GetName(), entry)) ;
247
248 //CreateTLists with arrays of TParticles. Filled with particles only relevant for the analysis.
249
250 TClonesArray * particleList = new TClonesArray("TParticle",1000); // All particles refitted in CTS and detected in EMCAL (jet)
251 TClonesArray * plCTS = new TClonesArray("TParticle",1000); // All particles refitted in Central Tracking System (ITS+TPC)
252 TClonesArray * plNe = new TClonesArray("TParticle",1000); // All particles measured in Jet Calorimeter
253 TClonesArray * plCalo = new TClonesArray("TParticle",1000); // All particles measured in Prompt Gamma calorimeter
254
255
256 TParticle *pGamma = new TParticle(); //It will contain the kinematics of the found prompt gamma
257 TParticle *pLeading = new TParticle(); //It will contain the kinematics of the found leading particle
258
259 Bool_t iIsInPHOSorEMCAL = kFALSE ; //To check if Gamma was in any calorimeter
260
261 //Fill lists with photons, neutral particles and charged particles
262 //look for the highest energy photon in the event inside fCalorimeter
263 AliDebug(2, "Fill particle lists, get prompt gamma");
264
265 //Fill particle lists
266
267 if(GetCalorimeter() == "PHOS")
268 CreateParticleList(particleList, plCTS,plNe,plCalo);
269 else if(GetCalorimeter() == "EMCAL")
270 CreateParticleList(particleList, plCTS,plCalo,plNe);
271 else
272 AliError("No calorimeter selected");
273
274 //Search highest energy prompt gamma in calorimeter
275 GetPromptGamma(plCalo, plCTS, pGamma, iIsInPHOSorEMCAL) ;
276
277 AliDebug(1, Form("Is Gamma in %s? %d",GetCalorimeter().Data(),iIsInPHOSorEMCAL));
278 AliDebug(3,Form("Charged list entries %d, Neutral list entries %d, %s list entries %d",
279 plCTS->GetEntries(),plNe->GetEntries(), GetCalorimeter().Data(), plCalo->GetEntries()));
280
281 //If there is any prompt photon in fCalorimeter,
282 //search jet leading particle
283
284 if(iIsInPHOSorEMCAL){
285 if (GetPrintInfo())
286 AliInfo(Form("Prompt Gamma: pt %f, phi %f, eta %f", pGamma->Pt(),pGamma->Phi(),pGamma->Eta())) ;
287
288 AliDebug(2, "Get Leading Particles, Make jet");
289
290 //Search leading particles in CTS and EMCAL
291 if(GetLeadingParticle(plCTS, plNe, pGamma, pLeading)){
292 if (GetPrintInfo())
293 AliInfo(Form("Leading: pt %f, phi %f, eta %f", pLeading->Pt(),pLeading->Phi(),pLeading->Eta())) ;
294
295 //Search Jet
296 if(!fSeveralConeAndPtCuts)
297 MakeJet(particleList,pGamma,pLeading,"");
298 else{
299 for(Int_t icone = 0; icone<fNCone; icone++) {
300 for(Int_t ipt = 0; ipt<fNPt;ipt++) {
301 TString lastname ="Cone"+ fNameCones[icone]+"Pt"+ fNamePtThres[ipt];
302 MakeJet(particleList, pGamma, pLeading,lastname);
303 }//icone
304 }//ipt
305 }//fSeveralConeAndPtCuts
306 }//Leading
307 }//Gamma in Calo
308
309 AliDebug(2, "End of analysis, delete pointers");
310
311 particleList->Delete() ;
312 plCTS->Delete() ;
313 plNe->Delete() ;
314 plCalo->Delete() ;
315 pLeading->Delete();
316 pGamma->Delete();
317
318 delete plNe ;
319 delete plCalo ;
320 delete plCTS ;
321 delete particleList ;
322 // delete pLeading;
323 // delete pGamma;
324
325 PostData(0, fOutputContainer);
326}
327
328//____________________________________________________________________________
329void AliAnaGammaJet::FillJetHistos(TClonesArray * pl, Double_t ptg, Double_t ptl, TString type, TString lastname)
330{
331 //Fill histograms wth jet fragmentation
332 //and number of selected jets and leading particles
333 //and the background multiplicity
334 TParticle * particle = 0 ;
335 Int_t ipr = 0;
336 Float_t charge = 0;
337
338 TIter next(pl) ;
339 while ( (particle = dynamic_cast<TParticle*>(next())) ) {
340 ipr++ ;
341 Double_t pt = particle->Pt();
342
343 charge = TDatabasePDG::Instance()
344 ->GetParticle(particle->GetPdgCode())->Charge();
345 if(charge != 0){//Only jet Charged particles
346 dynamic_cast<TH2F*>
347 (fOutputContainer->FindObject(type+"Fragment"+lastname))
348 ->Fill(ptg,pt/ptg);
349 dynamic_cast<TH2F*>
350 (fOutputContainer->FindObject(type+"PtDist"+lastname))
351 ->Fill(ptg,pt);
352 }//charged
353
354 }//while
355
356 if(type == "Bkg")
357 dynamic_cast<TH1F*>
358 (fOutputContainer->FindObject("NBkg"+lastname))
359 ->Fill(ipr);
360 else{
361 dynamic_cast<TH1F*>
362 (fOutputContainer->FindObject("NJet"+lastname))->
363 Fill(ptg);
364 dynamic_cast<TH2F*>
365 (fOutputContainer->FindObject("NLeading"+lastname))
366 ->Fill(ptg,ptl);
367 }
368
369}
370
371//____________________________________________________________________________
372void AliAnaGammaJet::GetLeadingCharge(TClonesArray * pl, TParticle * pGamma, TParticle * pLeading) const
373{
374 //Search for the charged particle with highest with
375 //Phi=Phi_gamma-Pi and pT=0.1E_gamma
376 Double_t pt = -100.;
377 Double_t phi = -100.;
378
379 for(Int_t ipr = 0;ipr < pl->GetEntries() ; ipr ++ ){
380
381 TParticle * particle = dynamic_cast<TParticle *>(pl->At(ipr)) ;
382
383 Double_t ptl = particle->Pt();
384 Double_t rat = ptl/pGamma->Pt() ;
385 Double_t phil = particle->Phi() ;
386 Double_t phig = pGamma->Phi();
387
388 //Selection within angular and energy limits
389 if(((phig-phil)> fPhiMinCut) && ((phig-phil)<fPhiMaxCut) &&
390 (rat > fRatioMinCut) && (rat < fRatioMaxCut) && (ptl > pt)) {
391 phi = phil ;
392 pt = ptl ;
393 pLeading->SetMomentum(particle->Px(),particle->Py(),particle->Pz(),particle->Energy());
394 AliDebug(4,Form("Charge in CTS: pt %f eta %f phi %f pt/Eg %f \n", pt, particle->Eta(), phi,rat)) ;
395 }
396 }
397
398 AliDebug(3,Form("Leading in CTS: pt %f eta %f phi %f pt/Eg %f \n", pt, pLeading->Eta(), phi,pt/pGamma->Pt())) ;
399
400}
401
402
403//____________________________________________________________________________
404void AliAnaGammaJet::GetLeadingPi0(TClonesArray * pl, TParticle * pGamma, TParticle * pLeading)
405{
406
407 //Search for the neutral pion with highest with
408 //Phi=Phi_gamma-Pi and pT=0.1E_gamma
409 Double_t pt = -100.;
410 Double_t phi = -100.;
411 Double_t ptl = -100.;
412 Double_t rat = -100.;
413 Double_t phil = -100. ;
414 Double_t ptg = pGamma->Pt();
415 Double_t phig = pGamma->Phi();
416
417 TIter next(pl);
418 TParticle * particle = 0;
419
420 Int_t iPrimary = -1;
421 TLorentzVector gammai,gammaj;
422 Double_t angle = 0., e = 0., invmass = 0.;
423 Double_t anglef = 0., ef = 0., invmassf = 0.;
424 Int_t ksPdg = 0;
425 Int_t jPrimary=-1;
426
427 while ( (particle = (TParticle*)next()) ) {
428 iPrimary++;
429
430 ksPdg = particle->GetPdgCode();
431 ptl = particle->Pt();
432 if(ksPdg == 111){ //2 gamma overlapped, found with PID
433 rat = ptl/ptg ;
434 phil = particle->Phi() ;
435 //Selection within angular and energy limits
436 if((ptl> pt)&& (rat > fRatioMinCut) && (rat < fRatioMaxCut) &&
437 ((phig-phil)>fPhiMinCut)&&((phig-phil)<fPhiMaxCut)){
438 phi = phil ;
439 pt = ptl ;
440 pLeading->SetMomentum(particle->Px(),particle->Py(),particle->Pz(),particle->Energy());
441 AliDebug(4,Form("Pi0 candidate: pt %f eta %f phi %f pt/Eg %f \n", pLeading->Pt(), pLeading->Eta(), pLeading->Phi(), pLeading->Pt()/pGamma->Pt())) ;
442 }// cuts
443 }// pdg = 111
444 else if(ksPdg == 22){//1 gamma
445 //Search the photon companion in case it comes from a Pi0 decay
446 //Apply several cuts to select the good pair
447 particle->Momentum(gammai);
448 jPrimary=-1;
449 TIter next2(pl);
450 while ( (particle = (TParticle*)next2()) ) {
451 jPrimary++;
452 if(jPrimary>iPrimary){
453 ksPdg = particle->GetPdgCode();
454
455 if(ksPdg == 22){
456 particle->Momentum(gammaj);
457 //Info("GetLeadingPi0","Egammai %f, Egammaj %f",
458 //gammai.Pt(),gammaj.Pt());
459 ptl = (gammai+gammaj).Pt();
460 phil = (gammai+gammaj).Phi();
461 if(phil < 0)
462 phil+=TMath::TwoPi();
463 rat = ptl/ptg ;
464 invmass = (gammai+gammaj).M();
465 angle = gammaj.Angle(gammai.Vect());
466 //Info("GetLeadingPi0","Angle %f", angle);
467 e = (gammai+gammaj).E();
468 //Fill histograms with no cuts applied.
469 fhAnglePairNoCut->Fill(e,angle);
470 fhInvMassPairNoCut->Fill(ptg,invmass);
471 //First cut on the energy and azimuth of the pair
472 if((rat > fRatioMinCut) && (rat < fRatioMaxCut) &&
473 ((phig-phil)>fPhiMinCut)&&((phig-phil)<fPhiMaxCut)){
474
475 fhAnglePairLeadingCut->Fill(e,angle);
476 fhInvMassPairLeadingCut->Fill(ptg,invmass);
477 //Second cut on the aperture of the pair
478 if(IsAngleInWindow(angle,e)){
479 fhAnglePairAngleCut->Fill(e,angle);
480 fhInvMassPairAngleCut->Fill(ptg,invmass);
481
482 //Info("GetLeadingPi0","InvMass %f", invmass);
483 //Third cut on the invariant mass of the pair
484 if((invmass>fInvMassMinCut) && (invmass<fInvMassMaxCut)){
485 fhInvMassPairAllCut->Fill(ptg,invmass);
486 fhAnglePairAllCut->Fill(e,angle);
487 if(ptl > pt ){
488 pt = ptl;
489 phi = phil ;
490 ef = e ;
491 anglef = angle ;
492 invmassf = invmass ;
493 pLeading->SetMomentum(particle->Px(),particle->Py(),particle->Pz(),particle->Energy());
494 AliDebug(4,Form("Pi0 candidate: pt %f eta %f phi %f pt/Eg %f \n", pLeading->Pt(), pLeading->Eta(), pLeading->Phi(), pLeading->Pt()/pGamma->Pt())) ;
495 }
496 }//cuts
497 }//(invmass>0.125) && (invmass<0.145)
498 }//gammaj.Angle(gammai.Vect())<0.04
499 }//if pdg = 22
500 }//iprimary<jprimary
501 }//while
502 }// if pdg = 22
503 // }
504 }//while
505
506 if(ef > 0.){//Final pi0 found, highest pair energy, fill histograms
507 fhInvMassPairLeading->Fill(ptg,invmassf);
508 fhAnglePairLeading->Fill(ef,anglef);
509 }
510
511 AliDebug(3,Form("Leading EMCAL: pt %f eta %f phi %f pt/Eg %f \n", pLeading->Pt(), pLeading->Eta(), pLeading->Phi(), pLeading->Pt()/pGamma->Pt())) ;
512}
513
514//____________________________________________________________________________
515Bool_t AliAnaGammaJet::GetLeadingParticle(TClonesArray * plCTS, TClonesArray * plNe,
516 TParticle * pGamma, TParticle * pLeading)
517{
518 //Search Charged or Neutral leading particle, select the highest one.
519
520 TParticle * pLeadingCh = new TParticle();
521 TParticle * pLeadingPi0 = new TParticle();
522
523 Double_t ptg = pGamma->Pt();
524 Double_t phig = pGamma->Phi();
525 Double_t etag = pGamma->Eta();
526
527 if(GetCalorimeter() == "PHOS" && !fJetsOnlyInCTS)
528 {
529 AliDebug(3, "GetLeadingPi0");
530 GetLeadingPi0 (plNe, pGamma, pLeadingPi0) ;
531 AliDebug(3, "GetLeadingCharge");
532 GetLeadingCharge(plCTS, pGamma, pLeadingCh) ;
533
534 Double_t ptch = pLeadingCh->Pt();
535 Double_t phich = pLeadingCh->Phi();
536 Double_t etach = pLeadingCh->Eta();
537 Double_t ptpi = pLeadingPi0->Pt();
538 Double_t phipi = pLeadingPi0->Phi();
539 Double_t etapi = pLeadingPi0->Eta();
540
541 //Is leading cone inside EMCAL acceptance?
542
543 Bool_t insidech = kFALSE ;
544 if((phich - fCone) > fPhiEMCALCut[0] &&
545 (phich + fCone) < fPhiEMCALCut[1] &&
546 (etach-fCone) < fEtaEMCALCut )
547 insidech = kTRUE ;
548
549 Bool_t insidepi = kFALSE ;
550 if((phipi - fCone) > fPhiEMCALCut[0] &&
551 (phipi + fCone) < fPhiEMCALCut[1] &&
552 (etapi-fCone) < fEtaEMCALCut )
553 insidepi = kTRUE ;
554
555 AliDebug(2,Form("Leading: charged pt %f, pi0 pt %f",ptch,ptpi)) ;
556
557 if (ptch > 0 || ptpi > 0)
558 {
559 if(insidech && (ptch > ptpi))
560 {
561 if (GetPrintInfo())
562 AliInfo("Leading found in CTS");
563 pLeading->SetMomentum(pLeadingCh->Px(),pLeadingCh->Py(),pLeadingCh->Pz(),pLeadingCh->Energy());
564 AliDebug(3,Form("Final leading found in CTS, pt %f, phi %f, eta %f",ptch,phich,etach)) ;
565 fhChargeRatio->Fill(ptg,ptch/pGamma->Pt());
566 fhDeltaPhiCharge->Fill(ptg,pGamma->Phi()-phich);
567 fhDeltaEtaCharge->Fill(ptg,pGamma->Eta()-etach);
568 return 1;
569 }
570
571 else if((ptpi > ptch) && insidepi)
572 {
573 if (GetPrintInfo())
574 AliInfo("Leading found in EMCAL");
575 pLeading->SetMomentum(pLeadingPi0->Px(),pLeadingPi0->Py(),pLeadingPi0->Pz(),pLeadingPi0->Energy());
576 AliDebug(3,Form("Final leading found in EMCAL, pt %f, phi %f, eta %f",ptpi,phipi,etapi)) ;
577 fhPi0Ratio ->Fill(ptg,ptpi/ptg);
578 fhDeltaPhiPi0->Fill(ptg,phig-phipi);
579 fhDeltaEtaPi0->Fill(ptg,etag-etapi);
580 return 1;
581 }
582
583 else{
584 AliDebug(3,"NO LEADING PARTICLE FOUND");}
585 return 0;
586 }
587 else{
588 AliDebug(3,"NO LEADING PARTICLE FOUND");
589 return 0;
590 }
591 }
592
593 else
594 {
595 //No calorimeter present for Leading particle detection
596 GetLeadingCharge(plCTS, pGamma, pLeading) ;
597 Double_t ptch = pLeading->Pt();
598 Double_t phich = pLeading->Phi();
599 Double_t etach = pLeading->Eta();
600 if(ptch > 0){
601 fhChargeRatio->Fill(ptg,ptch/ptg);
602 fhDeltaPhiCharge->Fill(ptg,phig-phich);
603 fhDeltaEtaCharge->Fill(ptg,etag-etach);
604 AliDebug(3,Form("Leading found : pt %f, phi %f, eta %f",ptch,phich,etach)) ;
605 return 1;
606 }
607 else
608 {
609 AliDebug(3,"NO LEADING PARTICLE FOUND");
610 return 0;
611 }
612 }
613}
614
615 //____________________________________________________________________________
616void AliAnaGammaJet::Init(const Option_t * )
617{
618// // Initialisation of branch container
619 AliAnaGammaDirect::Init();
620
621 //Initialize the parameters of the analysis.
622 //fCalorimeter="PHOS";
623 fAngleMaxParam.Set(4) ;
624 fAngleMaxParam.AddAt(0.4,0);//={0.4,-0.25,0.025,-2e-4};
625 fAngleMaxParam.AddAt(-0.25,1) ;
626 fAngleMaxParam.AddAt(0.025,2) ;
627 fAngleMaxParam.AddAt(-2e-4,3) ;
628 fSeveralConeAndPtCuts = kFALSE ;
629 //fPrintInfo = kTRUE;
630 fPbPb = kFALSE ;
631 fInvMassMaxCut = 0.16 ;
632 fInvMassMinCut = 0.11 ;
633 //fJetsOnlyInCTS = kTRUE ;
634 fEtaEMCALCut = 0.7 ;
635 fPhiEMCALCut[0] = 80. *TMath::Pi()/180.;
636 fPhiEMCALCut[1] = 190.*TMath::Pi()/180.;
637 fPhiMaxCut = 3.4 ;
638 fPhiMinCut = 2.9 ;
639
640 //Jet selection parameters
641 //Fixed cut (old)
642 fRatioMaxCut = 1.0 ;
643 fRatioMinCut = 0.1 ;
644 fJetRatioMaxCut = 1.2 ;
645 fJetRatioMinCut = 0.3 ;
646 fJetsOnlyInCTS = kFALSE ;
647 fJetCTSRatioMaxCut = 1.2 ;
648 fJetCTSRatioMinCut = 0.3 ;
649 fSelect = 0 ;
650
651 //Cut depending on gamma energy
652
653 fPtJetSelectionCut = 20.; //For Low pt jets+BKG, another limits applied
654 //Reconstructed jet energy dependence parameters
655 //e_jet = a1+e_gamma b2.
656 //Index 0-> Pt>2 GeV r = 0.3; Index 1-> Pt>0.5 GeV r = 0.3
657 fJetE1[0] = -5.75; fJetE1[1] = -4.1;
658 fJetE2[0] = 1.005; fJetE2[1] = 1.05;
659
660 //Reconstructed sigma of jet energy dependence parameters
661 //s_jet = a1+e_gamma b2.
662 //Index 0-> Pt>2 GeV r = 0.3; Index 1-> Pt>0.5 GeV r = 0.3
663 fJetSigma1[0] = 2.65; fJetSigma1[1] = 2.75;
664 fJetSigma2[0] = 0.0018; fJetSigma2[1] = 0.033;
665
666 //Background mean energy and RMS
667 //Index 0-> No BKG; Index 1-> BKG > 2 GeV;
668 //Index 2-> (low pt jets)BKG > 0.5 GeV;
669 //Index > 2, same for CTS conf
670 fBkgMean[0] = 0.; fBkgMean[1] = 8.8 ; fBkgMean[2] = 69.5;
671 fBkgMean[3] = 0.; fBkgMean[4] = 6.4; fBkgMean[5] = 48.6;
672 fBkgRMS[0] = 0.; fBkgRMS[1] = 7.5; fBkgRMS[2] = 22.0;
673 fBkgRMS[3] = 0.; fBkgRMS[4] = 5.4; fBkgRMS[5] = 13.2;
674
675 //Factor x of min/max = E -+ x * sigma. Obtained after selecting the
676 //limits for monoenergetic jets.
677 //Index 0-> No BKG; Index 1-> BKG > 2 GeV;
678 //Index 2-> (low pt jets) BKG > 0.5 GeV;
679 //Index > 2, same for CTS conf
680
681 fJetXMin1[0] =-0.69 ; fJetXMin1[1] = 0.39 ; fJetXMin1[2] =-0.88 ;
682 fJetXMin1[3] =-2.0 ; fJetXMin1[4] =-0.442 ; fJetXMin1[5] =-1.1 ;
683 fJetXMin2[0] = 0.066; fJetXMin2[1] = 0.038; fJetXMin2[2] = 0.034;
684 fJetXMin2[3] = 0.25 ; fJetXMin2[4] = 0.113; fJetXMin2[5] = 0.077 ;
685 fJetXMax1[0] =-3.8 ; fJetXMax1[1] =-0.76 ; fJetXMax1[2] =-3.6 ;
686 fJetXMax1[3] =-2.7 ; fJetXMax1[4] =-1.21 ; fJetXMax1[5] =-3.7 ;
687 fJetXMax2[0] =-0.012; fJetXMax2[1] =-0.022; fJetXMax2[2] = 0.016;
688 fJetXMax2[3] =-0.024; fJetXMax2[4] =-0.008; fJetXMax2[5] = 0.027;
689
690
691 //Different cones and pt thresholds to construct the jet
692
693 fCone = 0.3 ;
694 fPtThreshold = 0. ;
695 fNCone = 3 ;
696 fNPt = 3 ;
697 fCones[1] = 0.2 ; fNameCones[1] = "02" ;
698 fPtThres[0] = 0.0 ; fNamePtThres[0] = "00" ;
699 fCones[0] = 0.3 ; fNameCones[0] = "03" ;
700 fPtThres[1] = 0.5 ; fNamePtThres[1] = "05" ;
701 fCones[2] = 0.4 ; fNameCones[2] = "04" ;
702 fPtThres[2] = 1.0 ; fNamePtThres[2] = "10" ;
703 //Fill particle lists when PID is ok
704 // fEMCALPID = kFALSE;
705 // fPHOSPID = kFALSE;
706
707 //Initialization of histograms
708
709 MakeHistos() ;
710
711}
712
713//__________________________________________________________________________-
714Bool_t AliAnaGammaJet::IsAngleInWindow(const Float_t angle,const Float_t e) {
715 //Check if the opening angle of the candidate pairs is inside
716 //our selection windowd
717
718 Bool_t result = kFALSE;
719 Double_t mpi0 = 0.1349766;
720 Double_t max = fAngleMaxParam.At(0)*TMath::Exp(fAngleMaxParam.At(1)*e)
721 +fAngleMaxParam.At(2)+fAngleMaxParam.At(3)*e;
722 Double_t arg = (e*e-2*mpi0*mpi0)/(e*e);
723 Double_t min = 100. ;
724 if(arg>0.)
725 min = TMath::ACos(arg);
726
727 if((angle<max)&&(angle>=min))
728 result = kTRUE;
729
730 return result;
731}
732
733//__________________________________________________________________________-
734Bool_t AliAnaGammaJet::IsJetSelected(const Double_t ptg, const Double_t ptj){
735 //Check if the energy of the reconstructed jet is within an energy window
736
737 Double_t par[6];
738 Double_t xmax[2];
739 Double_t xmin[2];
740
741 Int_t iCTS = 0;
742 if(fJetsOnlyInCTS)
743 iCTS = 3 ;
744
745 if(!fPbPb){
746 //Phythia alone, jets with pt_th > 0.2, r = 0.3
747 par[0] = fJetE1[0]; par[1] = fJetE2[0];
748 //Energy of the jet peak
749 //e_jet = fJetE1[0]+fJetE2[0]*e_gamma, simulation fit
750 par[2] = fJetSigma1[0]; par[3] = fJetSigma2[0];
751 //Sigma of the jet peak
752 //sigma_jet = fJetSigma1[0]+fJetSigma2[0]*e_gamma, simulation fit
753 par[4] = fBkgMean[0 + iCTS]; par[5] = fBkgRMS[0 + iCTS];
754 //Parameters reserved for PbPb bkg.
755 xmax[0] = fJetXMax1[0 + iCTS]; xmax[1] = fJetXMax2[0 + iCTS];
756 xmin[0] = fJetXMin1[0 + iCTS]; xmin[1] = fJetXMin2[0 + iCTS];
757 //Factor that multiplies sigma to obtain the best limits,
758 //by observation, of mono jet ratios (ptjet/ptg)
759 //X_jet = fJetX1[0]+fJetX2[0]*e_gamma
760
761 }
762 else{
763 if(ptg > fPtJetSelectionCut){
764 //Phythia +PbPb with pt_th > 2 GeV/c, r = 0.3
765 par[0] = fJetE1[0]; par[1] = fJetE2[0];
766 //Energy of the jet peak, same as in pp
767 //e_jet = fJetE1[0]+fJetE2[0]*e_gamma, simulation fit
768 par[2] = fJetSigma1[0]; par[3] = fJetSigma2[0];
769 //Sigma of the jet peak, same as in pp
770 //sigma_jet = fJetSigma1[0]+fJetSigma2[0]*e_gamma, simulation fit
771 par[4] = fBkgMean[1 + iCTS]; par[5] = fBkgRMS[1 + iCTS];
772 //Mean value and RMS of PbPb Bkg
773 xmax[0] = fJetXMax1[1 + iCTS]; xmax[1] = fJetXMax2[1 + iCTS];
774 xmin[0] = fJetXMin1[1 + iCTS]; xmin[1] = fJetXMin2[1 + iCTS];
775 //Factor that multiplies sigma to obtain the best limits,
776 //by observation, of mono jet ratios (ptjet/ptg) mixed with PbPb Bkg,
777 //pt_th > 2 GeV, r = 0.3
778 //X_jet = fJetX1[0]+fJetX2[0]*e_gamma
779
780 }
781 else{
782 //Phythia + PbPb with pt_th > 0.5 GeV/c, r = 0.3
783 par[0] = fJetE1[1]; par[1] = fJetE2[1];
784 //Energy of the jet peak, pt_th > 2 GeV/c, r = 0.3
785 //e_jet = fJetE1[0]+fJetE2[0]*e_gamma, simulation fit
786 par[2] = fJetSigma1[1]; par[3] = fJetSigma2[1];
787 //Sigma of the jet peak, pt_th > 2 GeV/c, r = 0.3
788 //sigma_jet = fJetSigma1[0]+fJetSigma2[0]*e_gamma, simulation fit
789 par[4] = fBkgMean[2 + iCTS]; par[5] = fBkgRMS[2 + iCTS];
790 //Mean value and RMS of PbPb Bkg in a 0.3 cone, pt > 2 GeV.
791 xmax[0] = fJetXMax1[2 + iCTS]; xmax[1] = fJetXMax2[2 + iCTS];
792 xmin[0] = fJetXMin1[2 + iCTS]; xmin[1] = fJetXMin2[2 + iCTS];
793 //Factor that multiplies sigma to obtain the best limits,
794 //by observation, of mono jet ratios (ptjet/ptg) mixed with PbPb Bkg,
795 //pt_th > 2 GeV, r = 0.3
796 //X_jet = fJetX1[0]+fJetX2[0]*e_gamma
797
798 }//If low pt jet in bkg
799 }//if Bkg
800
801 //Calculate minimum and maximum limits of the jet ratio.
802 Double_t min = CalculateJetRatioLimit(ptg, par, xmin);
803 Double_t max = CalculateJetRatioLimit(ptg, par, xmax);
804
805 AliDebug(3,Form("Jet selection? : Limits min %f, max %f, pt_jet %f, pt_gamma %f, pt_jet / pt_gamma %f",min,max,ptj,ptg,ptj/ptg));
806
807 if(( min < ptj/ptg ) && ( max > ptj/ptg))
808 return kTRUE;
809 else
810 return kFALSE;
811
812}
813
814//____________________________________________________________________________
815void AliAnaGammaJet::MakeHistos()
816{
817 // Create histograms to be saved in output file and
818 // stores them in fOutputContainer
819
820 fOutputContainer = new TObjArray(10000) ;
821
822 TObjArray * outputContainer =GetOutputContainer();
823 for(Int_t i = 0; i < outputContainer->GetEntries(); i++ )
824 fOutputContainer->Add(outputContainer->At(i)) ;
825
826 //
827 fhChargeRatio = new TH2F
828 ("ChargeRatio","p_{T leading charge} /p_{T #gamma} vs p_{T #gamma}",
829 120,0,120,120,0,1);
830 fhChargeRatio->SetYTitle("p_{T lead charge} /p_{T #gamma}");
831 fhChargeRatio->SetXTitle("p_{T #gamma} (GeV/c)");
832 fOutputContainer->Add(fhChargeRatio) ;
833
834 fhDeltaPhiCharge = new TH2F
835 ("DeltaPhiCharge","#phi_{#gamma} - #phi_{charge} vs p_{T #gamma}",
836 200,0,120,200,0,6.4);
837 fhDeltaPhiCharge->SetYTitle("#Delta #phi");
838 fhDeltaPhiCharge->SetXTitle("p_{T #gamma} (GeV/c)");
839 fOutputContainer->Add(fhDeltaPhiCharge) ;
840
841 fhDeltaEtaCharge = new TH2F
842 ("DeltaEtaCharge","#eta_{#gamma} - #eta_{charge} vs p_{T #gamma}",
843 200,0,120,200,-2,2);
844 fhDeltaEtaCharge->SetYTitle("#Delta #eta");
845 fhDeltaEtaCharge->SetXTitle("p_{T #gamma} (GeV/c)");
846 fOutputContainer->Add(fhDeltaEtaCharge) ;
847
848 //
849 if(!fJetsOnlyInCTS){
850 fhPi0Ratio = new TH2F
851 ("Pi0Ratio","p_{T leading #pi^{0}} /p_{T #gamma} vs p_{T #gamma}",
852 120,0,120,120,0,1);
853 fhPi0Ratio->SetYTitle("p_{T lead #pi^{0}} /p_{T #gamma}");
854 fhPi0Ratio->SetXTitle("p_{T #gamma} (GeV/c)");
855 fOutputContainer->Add(fhPi0Ratio) ;
856
857 fhDeltaPhiPi0 = new TH2F
858 ("DeltaPhiPi0","#phi_{#gamma} - #phi_{ #pi^{0}} vs p_{T #gamma}",
859 200,0,120,200,0,6.4);
860 fhDeltaPhiPi0->SetYTitle("#Delta #phi");
861 fhDeltaPhiPi0->SetXTitle("p_{T #gamma} (GeV/c)");
862 fOutputContainer->Add(fhDeltaPhiPi0) ;
863
864 fhDeltaEtaPi0 = new TH2F
865 ("DeltaEtaPi0","#eta_{#gamma} - #eta_{ #pi^{0}} vs p_{T #gamma}",
866 200,0,120,200,-2,2);
867 fhDeltaEtaPi0->SetYTitle("#Delta #eta");
868 fhDeltaEtaPi0->SetXTitle("p_{T #gamma} (GeV/c)");
869 fOutputContainer->Add(fhDeltaEtaPi0) ;
870
871 //
872 fhAnglePair = new TH2F
873 ("AnglePair",
874 "Angle between #pi^{0} #gamma pair vs p_{T #pi^{0}}",
875 200,0,50,200,0,0.2);
876 fhAnglePair->SetYTitle("Angle (rad)");
877 fhAnglePair->SetXTitle("E_{ #pi^{0}} (GeV/c)");
878 fOutputContainer->Add(fhAnglePair) ;
879
880 fhAnglePairAccepted = new TH2F
881 ("AnglePairAccepted",
882 "Angle between #pi^{0} #gamma pair vs p_{T #pi^{0}}, both #gamma in eta<0.7, inside window",
883 200,0,50,200,0,0.2);
884 fhAnglePairAccepted->SetYTitle("Angle (rad)");
885 fhAnglePairAccepted->SetXTitle("E_{ #pi^{0}} (GeV/c)");
886 fOutputContainer->Add(fhAnglePairAccepted) ;
887
888 fhAnglePairNoCut = new TH2F
889 ("AnglePairNoCut",
890 "Angle between all #gamma pair vs p_{T #pi^{0}}",200,0,50,200,0,0.2);
891 fhAnglePairNoCut->SetYTitle("Angle (rad)");
892 fhAnglePairNoCut->SetXTitle("E_{ #pi^{0}} (GeV/c)");
893 fOutputContainer->Add(fhAnglePairNoCut) ;
894
895 fhAnglePairLeadingCut = new TH2F
896 ("AnglePairLeadingCut",
897 "Angle between all #gamma pair that have a good phi and pt vs p_{T #pi^{0}}",
898 200,0,50,200,0,0.2);
899 fhAnglePairLeadingCut->SetYTitle("Angle (rad)");
900 fhAnglePairLeadingCut->SetXTitle("E_{ #pi^{0}} (GeV/c)");
901 fOutputContainer->Add(fhAnglePairLeadingCut) ;
902
903 fhAnglePairAngleCut = new TH2F
904 ("AnglePairAngleCut",
905 "Angle between all #gamma pair (angle + leading cut) vs p_{T #pi^{0}}"
906 ,200,0,50,200,0,0.2);
907 fhAnglePairAngleCut->SetYTitle("Angle (rad)");
908 fhAnglePairAngleCut->SetXTitle("E_{ #pi^{0}} (GeV/c)");
909 fOutputContainer->Add(fhAnglePairAngleCut) ;
910
911 fhAnglePairAllCut = new TH2F
912 ("AnglePairAllCut",
913 "Angle between all #gamma pair (angle + inv mass cut+leading) vs p_{T #pi^{0}}"
914 ,200,0,50,200,0,0.2);
915 fhAnglePairAllCut->SetYTitle("Angle (rad)");
916 fhAnglePairAllCut->SetXTitle("E_{ #pi^{0}} (GeV/c)");
917 fOutputContainer->Add(fhAnglePairAllCut) ;
918
919 fhAnglePairLeading = new TH2F
920 ("AnglePairLeading",
921 "Angle between all #gamma pair finally selected vs p_{T #pi^{0}}",
922 200,0,50,200,0,0.2);
923 fhAnglePairLeading->SetYTitle("Angle (rad)");
924 fhAnglePairLeading->SetXTitle("E_{ #pi^{0}} (GeV/c)");
925 fOutputContainer->Add(fhAnglePairLeading) ;
926
927 //
928 fhInvMassPairNoCut = new TH2F
929 ("InvMassPairNoCut","Invariant Mass of all #gamma pair vs p_{T #gamma}",
930 120,0,120,360,0,0.5);
931 fhInvMassPairNoCut->SetYTitle("Invariant Mass (GeV/c^{2})");
932 fhInvMassPairNoCut->SetXTitle("p_{T #gamma} (GeV/c)");
933 fOutputContainer->Add(fhInvMassPairNoCut) ;
934
935 fhInvMassPairLeadingCut = new TH2F
936 ("InvMassPairLeadingCut",
937 "Invariant Mass of #gamma pair (leading cuts) vs p_{T #gamma}",
938 120,0,120,360,0,0.5);
939 fhInvMassPairLeadingCut->SetYTitle("Invariant Mass (GeV/c^{2})");
940 fhInvMassPairLeadingCut->SetXTitle("p_{T #gamma} (GeV/c)");
941 fOutputContainer->Add(fhInvMassPairLeadingCut) ;
942
943 fhInvMassPairAngleCut = new TH2F
944 ("InvMassPairAngleCut",
945 "Invariant Mass of #gamma pair (angle cut) vs p_{T #gamma}",
946 120,0,120,360,0,0.5);
947 fhInvMassPairAngleCut->SetYTitle("Invariant Mass (GeV/c^{2})");
948 fhInvMassPairAngleCut->SetXTitle("p_{T #gamma} (GeV/c)");
949 fOutputContainer->Add(fhInvMassPairAngleCut) ;
950
951 fhInvMassPairAllCut = new TH2F
952 ("InvMassPairAllCut",
953 "Invariant Mass of #gamma pair (angle+invmass cut+leading) vs p_{T #gamma}",
954 120,0,120,360,0,0.5);
955 fhInvMassPairAllCut->SetYTitle("Invariant Mass (GeV/c^{2})");
956 fhInvMassPairAllCut->SetXTitle("p_{T #gamma} (GeV/c)");
957 fOutputContainer->Add(fhInvMassPairAllCut) ;
958
959 fhInvMassPairLeading = new TH2F
960 ("InvMassPairLeading",
961 "Invariant Mass of #gamma pair selected vs p_{T #gamma}",
962 120,0,120,360,0,0.5);
963 fhInvMassPairLeading->SetYTitle("Invariant Mass (GeV/c^{2})");
964 fhInvMassPairLeading->SetXTitle("p_{T #gamma} (GeV/c)");
965 fOutputContainer->Add(fhInvMassPairLeading) ;
966 }
967
968
969 if(!fSeveralConeAndPtCuts){
970
971 //Count
972 fhNBkg = new TH1F("NBkg","bkg multiplicity",9000,0,9000);
973 fhNBkg->SetYTitle("counts");
974 fhNBkg->SetXTitle("N");
975 fOutputContainer->Add(fhNBkg) ;
976
977 fhNLeading = new TH2F
978 ("NLeading","Accepted Jet Leading", 240,0,120,240,0,120);
979 fhNLeading->SetYTitle("p_{T charge} (GeV/c)");
980 fhNLeading->SetXTitle("p_{T #gamma}(GeV/c)");
981 fOutputContainer->Add(fhNLeading) ;
982
983 fhNJet = new TH1F("NJet","Accepted jets",240,0,120);
984 fhNJet->SetYTitle("N");
985 fhNJet->SetXTitle("p_{T #gamma}(GeV/c)");
986 fOutputContainer->Add(fhNJet) ;
987
988 //Ratios and Pt dist of reconstructed (not selected) jets
989 //Jet
990 fhJetRatio = new TH2F
991 ("JetRatio","p_{T jet lead}/p_{T #gamma} vs p_{T #gamma}",
992 240,0,120,200,0,10);
993 fhJetRatio->SetYTitle("p_{T jet lead #pi^{0}}/p_{T #gamma}");
994 fhJetRatio->SetXTitle("p_{T #gamma} (GeV/c)");
995 fOutputContainer->Add(fhJetRatio) ;
996
997 fhJetPt = new TH2F
998 ("JetPt", "p_{T jet lead} vs p_{T #gamma}",240,0,120,400,0,200);
999 fhJetPt->SetYTitle("p_{T jet}");
1000 fhJetPt->SetXTitle("p_{T #gamma} (GeV/c)");
1001 fOutputContainer->Add(fhJetPt) ;
1002
1003 //Bkg
1004
1005 fhBkgRatio = new TH2F
1006 ("BkgRatio","p_{T bkg lead}/p_{T #gamma} vs p_{T #gamma}",
1007 240,0,120,200,0,10);
1008 fhBkgRatio->SetYTitle("p_{T bkg lead charge}/p_{T #gamma}");
1009 fhBkgRatio->SetXTitle("p_{T #gamma} (GeV/c)");
1010 fOutputContainer->Add(fhBkgRatio) ;
1011
1012 fhBkgPt = new TH2F
1013 ("BkgPt","p_{T jet lead} vs p_{T #gamma}",240,0,120,400,0,200);
1014 fhBkgPt->SetYTitle("p_{T jet lead charge}/p_{T #gamma}");
1015 fhBkgPt->SetXTitle("p_{T #gamma} (GeV/c)");
1016 fOutputContainer->Add(fhBkgPt) ;
1017
1018 //Jet Distributions
1019
1020 fhJetFragment =
1021 new TH2F("JetFragment","x = p_{T i charged}/p_{T #gamma}",
1022 240,0.,120.,1000,0.,1.2);
1023 fhJetFragment->SetYTitle("x_{T}");
1024 fhJetFragment->SetXTitle("p_{T #gamma}");
1025 fOutputContainer->Add(fhJetFragment) ;
1026
1027 fhBkgFragment = new TH2F
1028 ("BkgFragment","x = p_{T i charged}/p_{T #gamma}",
1029 240,0.,120.,1000,0.,1.2);
1030 fhBkgFragment->SetYTitle("x_{T}");
1031 fhBkgFragment->SetXTitle("p_{T #gamma}");
1032 fOutputContainer->Add(fhBkgFragment) ;
1033
1034 fhJetPtDist =
1035 new TH2F("JetPtDist","x = p_{T i charged}",240,0.,120.,400,0.,200.);
1036 fhJetPtDist->SetXTitle("p_{T #gamma} (GeV/c)");
1037 fOutputContainer->Add(fhJetPtDist) ;
1038
1039 fhBkgPtDist = new TH2F
1040 ("BkgPtDist","x = p_{T i charged}",240,0.,120.,400,0.,200.);
1041 fhBkgPtDist->SetXTitle("p_{T #gamma} (GeV/c)");
1042 fOutputContainer->Add(fhBkgPtDist) ;
1043
1044 }
1045 else{
1046 //If we want to study the jet for different cones and pt
1047
1048 for(Int_t icone = 0; icone<fNCone; icone++){
1049 for(Int_t ipt = 0; ipt<fNPt;ipt++){
1050
1051 //Jet
1052
1053 fhJetRatios[icone][ipt] = new TH2F
1054 ("JetRatioCone"+fNameCones[icone]+"Pt"+fNamePtThres[ipt],
1055 "p_{T jet lead #pi^{0}}/p_{T #gamma} vs p_{T #gamma}, cone ="
1056 +fNameCones[icone]+", pt>" +fNamePtThres[ipt]+" GeV/c",
1057 240,0,120,200,0,10);
1058 fhJetRatios[icone][ipt]->
1059 SetYTitle("p_{T jet lead #pi^{0}}/p_{T #gamma}");
1060 fhJetRatios[icone][ipt]->SetXTitle("p_{T #gamma} (GeV/c)");
1061 fOutputContainer->Add(fhJetRatios[icone][ipt]) ;
1062
1063
1064 fhJetPts[icone][ipt] = new TH2F
1065 ("JetPtCone"+fNameCones[icone]+"Pt"+fNamePtThres[ipt],
1066 "p_{T jet lead #pi^{0}}/p_{T #gamma} vs p_{T #gamma}, cone ="
1067 +fNameCones[icone]+", pt>" +fNamePtThres[ipt]+" GeV/c",
1068 240,0,120,400,0,200);
1069 fhJetPts[icone][ipt]->
1070 SetYTitle("p_{T jet lead #pi^{0}}/p_{T #gamma}");
1071 fhJetPts[icone][ipt]->SetXTitle("p_{T #gamma} (GeV/c)");
1072 fOutputContainer->Add(fhJetPts[icone][ipt]) ;
1073
1074 //Bkg
1075 fhBkgRatios[icone][ipt] = new TH2F
1076 ("BkgRatioCone"+fNameCones[icone]+"Pt"+fNamePtThres[ipt],
1077 "p_{T bkg lead #pi^{0}}/p_{T #gamma} vs p_{T #gamma}, cone ="
1078 +fNameCones[icone]+", pt>" +fNamePtThres[ipt]+" GeV/c",
1079 240,0,120,200,0,10);
1080 fhBkgRatios[icone][ipt]->
1081 SetYTitle("p_{T bkg lead #pi^{0}}/p_{T #gamma}");
1082 fhBkgRatios[icone][ipt]->SetXTitle("p_{T #gamma} (GeV/c)");
1083 fOutputContainer->Add(fhBkgRatios[icone][ipt]) ;
1084
1085 fhBkgPts[icone][ipt] = new TH2F
1086 ("BkgPtCone"+fNameCones[icone]+"Pt"+fNamePtThres[ipt],
1087 "p_{T jet lead #pi^{0}}/p_{T #gamma} vs p_{T #gamma}, cone ="
1088 +fNameCones[icone]+", pt>" +fNamePtThres[ipt]+" GeV/c",
1089 240,0,120,400,0,200);
1090 fhBkgPts[icone][ipt]->
1091 SetYTitle("p_{T jet lead #pi^{0}}/p_{T #gamma}");
1092 fhBkgPts[icone][ipt]->SetXTitle("p_{T #gamma} (GeV/c)");
1093 fOutputContainer->Add(fhBkgPts[icone][ipt]) ;
1094
1095 //Counts
1096 fhNBkgs[icone][ipt] = new TH1F
1097 ("NBkgCone"+fNameCones[icone]+"Pt"+fNamePtThres[ipt],
1098 "bkg multiplicity cone ="+fNameCones[icone]+", pt>"
1099 +fNamePtThres[ipt]+" GeV/c",9000,0,9000);
1100 fhNBkgs[icone][ipt]->SetYTitle("counts");
1101 fhNBkgs[icone][ipt]->SetXTitle("N");
1102 fOutputContainer->Add(fhNBkgs[icone][ipt]) ;
1103
1104 fhNLeadings[icone][ipt] = new TH2F
1105 ("NLeadingCone"+fNameCones[icone]+"Pt"+fNamePtThres[ipt],
1106 "p_{T #gamma} vs p_{T #pi^{0}} cone ="+fNameCones[icone]+", pt>"
1107 +fNamePtThres[ipt]+" GeV/c",120,0,120,120,0,120);
1108 fhNLeadings[icone][ipt]->SetYTitle("p_{T #pi^{0}}(GeV/c)");
1109 fhNLeadings[icone][ipt]->SetXTitle("p_{T #gamma}(GeV/c)");
1110 fOutputContainer->Add(fhNLeadings[icone][ipt]) ;
1111
1112 fhNJets[icone][ipt] = new TH1F
1113 ("NJetCone"+fNameCones[icone]+"Pt"+fNamePtThres[ipt],
1114 "Number of neutral jets, cone ="+fNameCones[icone]+", pt>"
1115 +fNamePtThres[ipt]+" GeV/c",120,0,120);
1116 fhNJets[icone][ipt]->SetYTitle("N");
1117 fhNJets[icone][ipt]->SetXTitle("p_{T #gamma}(GeV/c)");
1118 fOutputContainer->Add(fhNJets[icone][ipt]) ;
1119
1120 //Fragmentation Function
1121 fhJetFragments[icone][ipt] = new TH2F
1122 ("JetFragmentCone"+fNameCones[icone]+"Pt"+fNamePtThres[ipt],
1123 "x_{T} = p_{T i}/p_{T #gamma}, cone ="+fNameCones[icone]+", pt>"
1124 +fNamePtThres[ipt]+" GeV/c",120,0.,120.,240,0.,1.2);
1125 fhJetFragments[icone][ipt]->SetYTitle("x_{T}");
1126 fhJetFragments[icone][ipt]->SetXTitle("p_{T #gamma}");
1127 fOutputContainer->Add(fhJetFragments[icone][ipt]) ;
1128
1129 fhBkgFragments[icone][ipt] = new TH2F
1130 ("BkgFragmentCone"+fNameCones[icone]+"Pt"+fNamePtThres[ipt],
1131 "x_{T} = p_{T i}/p_{T #gamma}, cone ="+fNameCones[icone]+", pt>"
1132 +fNamePtThres[ipt]+" GeV/c",120,0.,120.,240,0.,1.2);
1133 fhBkgFragments[icone][ipt]->SetYTitle("x_{T}");
1134 fhBkgFragments[icone][ipt]->SetXTitle("p_{T #gamma}");
1135 fOutputContainer->Add(fhBkgFragments[icone][ipt]) ;
1136
1137 //Jet particle distribution
1138
1139 fhJetPtDists[icone][ipt] = new TH2F
1140 ("JetPtDistCone"+fNameCones[icone]+"Pt"+fNamePtThres[ipt],
1141 "p_{T i}, cone ="+fNameCones[icone]+", pt>" +fNamePtThres[ipt]+
1142 " GeV/c",120,0.,120.,120,0.,120.);
1143 fhJetPtDists[icone][ipt]->SetXTitle("p_{T #gamma} (GeV/c)");
1144 fOutputContainer->Add(fhJetPtDists[icone][ipt]) ;
1145
1146 fhBkgPtDists[icone][ipt] = new TH2F
1147 ("BkgPtDistCone"+fNameCones[icone]+"Pt"+fNamePtThres[ipt],
1148 "p_{T i}, cone ="+fNameCones[icone]+", pt>" +fNamePtThres[ipt]+
1149 " GeV/c",120,0.,120.,120,0.,120.);
1150 fhBkgPtDists[icone][ipt]->SetXTitle("p_{T #gamma} (GeV/c)");
1151 fOutputContainer->Add(fhBkgPtDists[icone][ipt]) ;
1152
1153 }//ipt
1154 } //icone
1155 }//If we want to study any cone or pt threshold
1156}
1157
1158//____________________________________________________________________________
1159void AliAnaGammaJet::MakeJet(TClonesArray * pl, TParticle * pGamma, TParticle* pLeading,TString lastname)
1160{
1161 //Fill the jet with the particles around the leading particle with
1162 //R=fCone and pt_th = fPtThres. Caclulate the energy of the jet and
1163 //check if we select it. Fill jet histograms
1164
1165 TClonesArray * jetList = new TClonesArray("TParticle",1000);
1166 TClonesArray * bkgList = new TClonesArray("TParticle",1000);
1167
1168 TLorentzVector jet (0,0,0,0);
1169 TLorentzVector bkg(0,0,0,0);
1170 TLorentzVector lv (0,0,0,0);
1171
1172 Double_t ptjet = 0.0;
1173 Double_t ptbkg = 0.0;
1174 Int_t n0 = 0;
1175 Int_t n1 = 0;
1176 Bool_t b1 = kFALSE;
1177 Bool_t b0 = kFALSE;
1178
1179 Double_t ptg = pGamma->Pt();
1180 Double_t phig = pGamma->Phi();
1181 Double_t ptl = pLeading->Pt();
1182 Double_t phil = pLeading->Phi();
1183 Double_t etal = pLeading->Eta();
1184
1185 Float_t ptcut = 0. ;
1186 if(fPbPb){
1187 if(ptg > fPtJetSelectionCut) ptcut = 2. ;
1188 else ptcut = 0.5;
1189 }
1190
1191 TIter next(pl) ;
1192 TParticle * particle = 0 ;
1193 while ( (particle = dynamic_cast<TParticle*>(next())) ) {
1194
1195 b0 = kFALSE;
1196 b1 = kFALSE;
1197
1198 //Particles in jet
1199 SetJet(particle, b0, fCone, etal, phil) ;
1200
1201 if(b0){
1202 new((*jetList)[n0++]) TParticle(*particle) ;
1203 particle->Momentum(lv);
1204 if(particle->Pt() > ptcut ){
1205 jet+=lv;
1206 ptjet+=particle->Pt();
1207 }
1208 }
1209
1210 //Background around (phi_gamma-pi, eta_leading)
1211 SetJet(particle, b1, fCone,etal, phig) ;
1212
1213 if(b1) {
1214 new((*bkgList)[n1++]) TParticle(*particle) ;
1215 particle->Momentum(lv);
1216 if(particle->Pt() > ptcut ){
1217 bkg+=lv;
1218 ptbkg+=particle->Pt();
1219 }
1220 }
1221 }
1222
1223 ptjet = jet.Pt();
1224 ptbkg = bkg.Pt();
1225
1226 if(ptjet > 0.) {
1227
1228 AliDebug(2,Form("Gamma pt %f, Jet pt %f, Bkg pt %f",ptg,ptjet,ptbkg));
1229
1230 //Fill histograms
1231
1232 Double_t ratjet = ptjet/ptg ;
1233 Double_t ratbkg = ptbkg/ptg ;
1234
1235 dynamic_cast<TH2F*>
1236 (fOutputContainer->FindObject("JetRatio"+lastname))
1237 ->Fill(ptg,ratjet);
1238 dynamic_cast<TH2F*>
1239 (fOutputContainer->FindObject("JetPt"+lastname))
1240 ->Fill(ptg,ptjet);
1241
1242 dynamic_cast<TH2F*>
1243 (fOutputContainer->FindObject("BkgRatio"+lastname))
1244 ->Fill(ptg,ratbkg);
1245
1246 dynamic_cast<TH2F*>
1247 (fOutputContainer->FindObject("BkgPt"+lastname))
1248 ->Fill(ptg,ptbkg);
1249
1250
1251 //Jet selection
1252 Bool_t kSelect = kFALSE;
1253 if(fSelect == 0)
1254 kSelect = kTRUE; //Accept all jets, no restriction
1255 else if(fSelect == 1){
1256 //Selection with parametrized cuts
1257 if(IsJetSelected(ptg,ptjet)) kSelect = kTRUE;
1258 }
1259 else if(fSelect == 2){
1260 //Simple selection
1261 if(!fJetsOnlyInCTS){
1262 if((ratjet < fJetRatioMaxCut) && (ratjet > fJetRatioMinCut )) kSelect = kTRUE;
1263 }
1264 else{
1265 if((ratjet < fJetCTSRatioMaxCut) && (ratjet > fJetCTSRatioMinCut )) kSelect = kTRUE;
1266 }
1267 }
1268 else
1269 AliError("Jet selection option larger than 2, DONT SELECT JETS");
1270
1271
1272 if(kSelect){
1273 if (GetPrintInfo())
1274 AliInfo(Form("Jet Selected: pt %f ", ptjet)) ;
1275
1276 FillJetHistos(jetList, ptg, ptl,"Jet",lastname);
1277 FillJetHistos(bkgList, ptg, ptl, "Bkg",lastname);
1278 }
1279 } //ptjet > 0
1280
1281 jetList ->Delete();
1282 bkgList ->Delete();
1283
1284}
1285
1286//____________________________________________________________________________
1287void AliAnaGammaJet::Print(const Option_t * opt) const
1288{
1289
1290 //Print some relevant parameters set for the analysis
1291 if(! opt)
1292 return;
1293
1294 Info("Print", "%s %s", GetName(), GetTitle() ) ;
1295 if(!fJetsOnlyInCTS)
1296 printf("EMCAL Acceptance cut : phi min %f, phi max %f, eta %f \n", fPhiEMCALCut[0], fPhiEMCALCut[1], fEtaEMCALCut) ;
1297
1298 printf("Phi_Leading < %f\n", fPhiMaxCut) ;
1299 printf("Phi_Leading > %f\n", fPhiMinCut) ;
1300 printf("pT Leading / pT Gamma < %f\n", fRatioMaxCut) ;
1301 printf("pT Leading / pT Gamma > %f\n", fRatioMinCut) ;
1302 printf("M_pair < %f\n", fInvMassMaxCut) ;
1303 printf("M_pair > %f\n", fInvMassMinCut) ;
1304 if(fSelect == 2){
1305 printf("pT Jet / pT Gamma < %f\n", fJetRatioMaxCut) ;
1306 printf("pT Jet / pT Gamma > %f\n", fJetRatioMinCut) ;
1307 printf("pT Jet (Only CTS)/ pT Gamma < %f\n", fJetCTSRatioMaxCut) ;
1308 printf("pT Jet (Only CTS)/ pT Gamma > %f\n", fJetCTSRatioMinCut) ;
1309 }
1310
1311
1312}
1313
1314//___________________________________________________________________
1315void AliAnaGammaJet::SetJet(TParticle * part, Bool_t & b, Float_t cone,
1316 Double_t eta, Double_t phi)
1317{
1318
1319 //Check if the particle is inside the cone defined by the leading particle
1320 b = kFALSE;
1321
1322 if(phi > TMath::TwoPi())
1323 phi-=TMath::TwoPi();
1324 if(phi < 0.)
1325 phi+=TMath::TwoPi();
1326
1327 Double_t rad = 10000 + cone;
1328
1329 if(TMath::Abs(part->Phi()-phi) <= (TMath::TwoPi() - cone))
1330 rad = TMath::Sqrt(TMath::Power(part->Eta()-eta,2)+
1331 TMath::Power(part->Phi()-phi,2));
1332 else{
1333 if(part->Phi()-phi > TMath::TwoPi() - cone)
1334 rad = TMath::Sqrt(TMath::Power(part->Eta()-eta,2)+
1335 TMath::Power((part->Phi()-TMath::TwoPi())-phi,2));
1336 if(part->Phi()-phi < -(TMath::TwoPi() - cone))
1337 rad = TMath::Sqrt(TMath::Power(part->Eta()-eta,2)+
1338 TMath::Power((part->Phi()+TMath::TwoPi())-phi,2));
1339 }
1340
1341 if(rad < cone )
1342 b = kTRUE;
1343
1344}
1345
1346
1347void AliAnaGammaJet::Terminate(Option_t *)
1348{
1349 // The Terminate() function is the last function to be called during
1350 // a query. It always runs on the client, it can be used to present
1351 // the results graphically or save the results to file.
1352
1353
1354}