changed Resolution, Material, and PhotonQA task to be able to run on the grid
[u/mrichter/AliRoot.git] / PWGGA / GammaConv / AliAnaConvIsolation.cxx
1 /**************************************************************************
2  * This file is property of and copyright by the ALICE HLT Project        *
3  * ALICE Experiment at CERN, All rights reserved.                         *
4  *                                                                        *
5  * Primary Author: Svein Lindal <slindal@fys.uio.no>                      *
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
16 /// @file   AliAnalysisTaskGammaJet.cxx
17 /// @author Svein Lindal
18 /// @brief  Class used to run isolation studies of conversion gamma/pions
19
20
21 #include "AliAnaConvIsolation.h"
22 #include "AliAODTrack.h"
23
24 #include "TObject.h"
25 #include "TClonesArray.h"
26 #include "TH2F.h"
27 #include "TList.h"
28 #include "AliAODConversionPhoton.h"
29
30
31 #include <iostream>
32
33 // Gamma - jet correlation analysis task
34 // Author: Svein Lindal
35
36
37 using namespace std;
38 ClassImp(AliAnaConvIsolation)
39
40 //________________________________________________________________________
41 AliAnaConvIsolation::AliAnaConvIsolation () : TObject(),
42   fIsoCurve(NULL),
43   fCurveFunction("0.1*x"),
44   fConeSize(0), 
45   fMinPt(0.1),
46   fMaxPtThreshold(0),
47   fSumPtThreshold(0),
48   fMaxPtFraction(0),
49   fSumPtFraction(0)
50   // fHistograms(NULL),
51   // fHistogramMaxPt(50)  
52 {
53   //Constructor
54   // for(Int_t i = 0; i < 2; i++){
55   //   fhMaxPtInCone[i] = NULL;
56   //   fhSumPtInCone[i] = NULL;
57   //   fhSumPtVsMaxPt[i] = NULL;
58   //   fhPtCandidates[i] = NULL;
59   //   fhTrackMult[i] = NULL;
60   // }
61   
62   
63 }
64
65 //________________________________________________________________________
66 AliAnaConvIsolation::AliAnaConvIsolation(Float_t coneSize, Float_t maxPtThreshold, Float_t sumPtThreshold, Float_t maxPtFraction, Float_t sumPtFraction) :
67   TObject(), 
68   fIsoCurve(NULL),
69   fCurveFunction("0.1*x"),
70   fConeSize(coneSize), 
71   fMinPt(0.1),
72   fMaxPtThreshold(maxPtThreshold),
73   fSumPtThreshold(sumPtThreshold),
74   fMaxPtFraction(maxPtFraction),
75   fSumPtFraction(sumPtFraction)
76   // fHistograms(NULL),
77   // fHistogramMaxPt(50)
78 {
79   //Constructor
80   // for(Int_t i = 0; i < 2; i++){
81   //   fhMaxPtInCone[i] = NULL;
82   //   fhSumPtInCone[i] = NULL;
83   //   fhSumPtVsMaxPt[i] = NULL;
84   //   fhPtCandidates[i] = NULL;
85   //   fhTrackMult[i] = NULL;
86   // }
87
88   fIsoCurve = new TF1("Isolation_curve", fCurveFunction.Data(), 0, 100);
89
90 }
91
92
93 //________________________________________________________________________________
94 AliAnaConvIsolation::~AliAnaConvIsolation() {
95   //Destructor
96
97 }
98
99
100 //________________________________________________________________________
101 // void AliAnaConvIsolation::CreateHistograms()
102 // {
103         //Create histograms
104   // if(!fHistograms) fHistograms = new TList();
105   // fHistograms->SetName(Form("Isolation_histo_cone_%f_maxPt_%f_sumPt_%f", fConeSize, fSumPtThreshold, fMaxPtThreshold));
106
107   //fHistograms->Add(fIsoCurve);
108
109   //cout << "Creatin isolation histograms conesize :" << fConeSize << endl;
110
111   //Create histograms add, to outputlis
112   // for(Int_t i = 0; i < 2; i++){
113   //   fhMaxPtInCone[i] = new TH2F(Form("fhMaxPtInCone_%s_%f", (i==0)? "nonIso" : "isolated", fConeSize), 
114   //                            Form("Max pt nonIso particle in cone %f vs candidate pt", fConeSize), 
115   //                            200, 0, fHistogramMaxPt, 200, 0, fHistogramMaxPt);
116   //   fHistograms->Add(fhMaxPtInCone[i]);
117     
118   //   fhSumPtInCone[i] = new TH2F(Form("fhSumPtInCone_%s_%f",  (i==0)? "nonIso" : "isolated", fConeSize), 
119   //                         Form("Sum pt in cone %f vs candidate pt %s", fConeSize,  (i==0)? "nonIsoay" : "isolated"), 
120   //                         200, 0, fHistogramMaxPt, 200, 0, fHistogramMaxPt);
121   //   fHistograms->Add(fhSumPtInCone[i]);
122
123   //   fhSumPtVsMaxPt[i] = new TH2F(Form("fhSumPtVsMaxPt_%s_%f",  (i==0)? "nonIso" : "isolated", fConeSize), 
124   //                             Form("fhSumPtVsMaxPt_%s_%f",  (i==0)? "nonIso" : "isolated", fConeSize), 
125   //                             200, 0, fHistogramMaxPt, 200, 0, fHistogramMaxPt);
126   //   fHistograms->Add(fhSumPtVsMaxPt[i]);
127   // }
128
129   // for(Int_t iIso = 0; iIso < 2; iIso++){
130   //   fhPtCandidates[iIso] = new TH1F(Form("fhPtCandidates_%s_%f", (iIso==0)? "nonIso" : "isolated", fConeSize),
131   //                                Form("Pt of %s candidates, cone size %f", (iIso==0)? "nonIsolated" : "isolated", fConeSize),
132   //                             200, 0 , fHistogramMaxPt);
133   //   fhPtCandidates[iIso]->GetXaxis()->SetTitle("P_{T} (GeV/c)");
134   //   fhPtCandidates[iIso]->GetYaxis()->SetTitle("dN/dP_{T} (c/GeV)");
135   //   fhPtCandidates[iIso]->SetMarkerStyle(kFullCircle);
136   //   fHistograms->Add(fhPtCandidates[iIso]);
137   // }
138
139   // for(int iDec = 0; iDec < 2; iDec++) {
140   //   for(int iIso = 0; iIso < 2; iIso++) {
141
142   //     fHistSumPt[iIso][iDec] = new TH1F(Form("fHistSumPt_%f_%s_%s", fConeSize, (iIso==0)?"nonIso":"isolated" , (iDec==0)?"noDec":"decay" ),  
143   //                                    Form("P_{T} distribution cone %f %s %s", fConeSize, (iIso==0)?"nonIso":"isolated" , (iDec==0)?"noDec":"decay" ), 
144   //                                 150, 0.1, 50);
145   //     fHistSumPt[iIso][iDec]->GetXaxis()->SetTitle("P_{T} (GeV/c)");
146   //     fHistSumPt[iIso][iDec]->GetYaxis()->SetTitle("dN/dP_{T} (c/GeV)");
147   //     fHistSumPt[iIso][iDec]->SetMarkerStyle(kFullCircle);
148   //     fHistograms->Add(fHistSumPt[iIso][iDec]);
149      
150       
151   //     fHistMaxPt[iIso][iDec] = new TH1F(Form("fHistMaxPt_%f_%s_%s", fConeSize, (iIso==0)?"nonIso":"isolated" , (iDec==0)?"noDec":"decay" ),  
152   //                                    Form("P_{T} distribution cone %f %s %s", fConeSize, (iIso==0)?"nonIso":"isolated" , (iDec==0)?"noDec":"decay" ), 
153   //                                 150, 0.1, 50);
154   //     fHistMaxPt[iIso][iDec]->GetXaxis()->SetTitle("P_{T} (GeV/c)");
155   //     fHistMaxPt[iIso][iDec]->GetYaxis()->SetTitle("dN/dP_{T} (c/GeV)");
156   //     fHistMaxPt[iIso][iDec]->SetMarkerStyle(kFullCircle);
157   //     fHistograms->Add(fHistMaxPt[iIso][iDec]);
158      
159  
160   //   }
161   // }
162
163
164   // for(Int_t iIso = 0; iIso < 2; iIso++){
165   //   fhTrackMult[iIso] = new TH1F(Form("fhTrackMult_%s_%f", (iIso==0)? "nonIso" : "isolated", fConeSize),
166   //                                Form("Pt of %s candidates, cone size %f", (iIso==0)? "nonIsolated" : "isolated", fConeSize),
167   //                             150, 0 , 150);
168   //   fhTrackMult[iIso]->GetXaxis()->SetTitle("n tracks in event");
169   //   fhTrackMult[iIso]->GetYaxis()->SetTitle("dN/dNTracks");
170   //   fhTrackMult[iIso]->SetMarkerStyle(kFullCircle);
171   //   fHistograms->Add(fhTrackMult[iIso]);
172   // }
173
174 //}
175
176 //________________________________________________________________________
177 Bool_t AliAnaConvIsolation::IsLeading(AliAODConversionParticle * particle, const TObjArray * tracks, const TObjArray * aodParticles) {
178   //Check if particle is highest pt particle in cone
179   for(Int_t ip = 0; ip < aodParticles->GetEntriesFast(); ip++) {
180         AliAODConversionParticle * oparticle = static_cast<AliAODConversionParticle*>(aodParticles->UncheckedAt(ip));
181         if(oparticle == particle) continue;
182         if(oparticle->Pt() > particle->Pt() && 
183            IsInCone(particle->Eta() - oparticle->Eta(), 
184                                 particle->Phi() - oparticle->Phi(), 
185                                 fConeSize) ) {
186           return kFALSE;
187         }
188   }
189
190  for(Int_t it = 0; it < tracks->GetEntriesFast(); it++) {
191         AliAODTrack * track = static_cast<AliAODTrack*>(tracks->UncheckedAt(it));
192         if(track->Pt() > particle->Pt()  && 
193            IsInCone(particle->Eta() - track->Eta(), particle->Phi() - track->Phi(), fConeSize) ) return kFALSE;
194   }
195
196   return kTRUE;
197 }
198
199 ///________________________________________________________________________
200 Int_t AliAnaConvIsolation::IsLeading(const AliAODConversionParticle * particle, const TObjArray * tracks, const Int_t tIDs[4]) const {
201
202   Bool_t leadingEvent = kTRUE;
203
204   //Is there a higher pt particle within cone ?
205   Float_t tPt = particle->Pt();
206   for(int it = 0; it < tracks->GetEntriesFast(); it++) {
207     AliVTrack * track = static_cast<AliVTrack*>(tracks->UncheckedAt(it));
208         Int_t tid = track->GetID();
209         if(tid == tIDs[0] || tid == tIDs[1] || tid == tIDs[2] || tid == tIDs[3]) continue;
210         if(track->Pt() > tPt) {
211           if(IsInCone(particle->Eta() - track->Eta(), particle->Phi() - track->Phi(), 0.8)) return 0;
212           leadingEvent = kFALSE;  
213         }
214   }     
215
216   if(leadingEvent) {
217         return 2;
218   }  else {
219         return 1;
220   }
221
222   return 0;
223 }
224
225 //_________________________________________________________________________
226 Bool_t AliAnaConvIsolation::IsIsolated(AliAODConversionPhoton * particle, const TClonesArray * const tracks, const Int_t nSpawn, const Int_t * const spawn, Bool_t &leading) {
227   //See header file for documentation
228
229   leading = kTRUE;
230
231   Float_t ptSum = 0.0;
232   Float_t ptMax = 0.0;
233
234   for(int it = 0; it < tracks->GetEntriesFast(); it++) {
235     AliAODTrack * track = dynamic_cast<AliAODTrack*>(tracks->At(it));
236
237     if (track) {
238
239       if(track->Pt() < fMinPt) continue;
240
241                         if(nSpawn && spawn) { ;}
242
243       ///Ignore tracks that are grandchildren of pion
244       // if ( particle->IsMySpawn(track->GetID(), nSpawn, spawn)) 
245       //        continue;
246       
247       
248       if ( IsInCone(particle->Eta() - track->Eta(), particle->Phi() - track->Phi(), fConeSize) ) {
249         ptSum += track->Pt();
250         if(track->Pt() > ptMax) ptMax = track->Pt();
251         if(track->Pt() > particle->Pt()) leading = kFALSE;
252       }
253     } else {
254       cout << "Bad track"<<endl;
255     }
256   }
257   
258
259   Bool_t isolated = EvaluateIsolationCriteria( ptSum, particle->Pt());
260
261   //FillHistograms(particle->Pt(), ptMax, ptSum, isolated, tracks->GetEntriesFast());
262
263   return isolated;
264
265 }
266
267
268 //_________________________________________________________________________
269 Bool_t AliAnaConvIsolation::IsIsolated(const AliAODConversionPhoton * const particle, const TClonesArray * const tracks, Bool_t &leading ) {
270   //See header file for documentation
271
272   leading = kTRUE;
273
274   Float_t ptSum = 0.0;
275   Float_t ptMax = 0.0;
276
277   for(int it = 0; it < tracks->GetEntriesFast(); it++) {
278   
279     AliAODTrack * track = dynamic_cast<AliAODTrack*>(tracks->At(it));
280    
281     if (track) {
282     
283       if(track->Pt() < fMinPt) continue;
284       
285       if ( (track->GetID() == particle->GetTrackLabel(0)) || track->GetID() == particle->GetTrackLabel(1) )  
286         continue;
287       
288       if ( IsInCone(particle->Eta() - track->Eta(), particle->Phi() - track->Phi(), fConeSize) ) {
289         ptSum += track->Pt();
290         if(track->Pt() > ptMax) ptMax = track->Pt();
291         if(track->Pt() > particle->Pt()) leading = kFALSE;
292       }
293     } else {
294       cout << "Bad track"<<endl;
295     }
296   }
297   
298
299   Bool_t isolated = EvaluateIsolationCriteria( ptSum, particle->Pt());
300
301   //FillHistograms(particle->Pt(), ptMax, ptSum, isolated, tracks->GetEntriesFast());
302
303   
304   return isolated;
305
306 }
307
308 // ///___________________________________________________________________________
309 // void AliAnaConvIsolation::FillHistograms(Float_t pt, Float_t ptMax, Float_t ptSum, Bool_t isolated, Int_t nTracks) {
310 //   //Fill histograms
311 //   fhMaxPtInCone[isolated]->Fill(pt, ptMax);
312 //   fhSumPtInCone[isolated]->Fill(pt, ptSum);
313 //   fhSumPtVsMaxPt[isolated]->Fill(ptMax, ptSum);
314 //   fhPtCandidates[isolated]->Fill(pt);
315 //   fhTrackMult[isolated]->Fill(nTracks);
316 // }
317
318
319 ///_____________________________________________________________________________
320 Bool_t AliAnaConvIsolation::EvaluateIsolationCriteria(Float_t ptSum, Float_t pt) const {
321   //Evaluate isolation criteria
322   return (ptSum < fIsoCurve->Eval(pt));
323
324 }