AliPhysics  63e47e1 (63e47e1)
AliPhotonIsolation.cxx
Go to the documentation of this file.
1 
2 /**************************************************************************
3  * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
4  **
5  * Authors: Helena Horstmann, Daniel Mühlheim*
6  * Version 1.0*
7  **
8  * Permission to use, copy, modify and distribute this software and its *
9  * documentation strictly for non-commercial purposes is hereby granted *
10  * without fee, provided that the above copyright notice appears in all *
11  * copies and that both the copyright notice and this permission notice *
12  * appear in the supporting documentation. The authors make no claims*
13  * about the suitability of this software for any purpose. It is*
14  * provided "as is" without express or implied warranty.*
15  **************************************************************************/
16 
18 //---------------------------------------------
19 // Basic Isolation Class
20 //---------------------------------------------
22 
23 
24 #include "AliAnalysisManager.h"
25 #include "AliAODEvent.h"
26 #include "AliESDEvent.h"
27 #include "AliMCEventHandler.h"
28 #include "AliMCEvent.h"
29 #include "AliMCParticle.h"
30 #include "AliAODMCParticle.h"
31 #include "AliPhotonIsolation.h"
32 #include "AliV0ReaderV1.h"
33 #include "AliVParticle.h"
34 #include "TChain.h"
35 #include "TH1F.h"
36 #include "TAxis.h"
37 
38 #include <vector>
39 #include <map>
40 #include <utility>
41 
42 
43 class iostream;
44 
45 using namespace std;
46 
47 
48 ClassImp(AliPhotonIsolation)
49 
50 //_______________________________________________________________________
51 AliPhotonIsolation::AliPhotonIsolation(const char *name, Int_t photonType) : AliAnalysisTaskSE(name),
52 
53  fPhotonType(photonType),
54  fV0ReaderName(""),
55  fCorrTaskSetting(""),
56  fMapClustertoPtR1(),
57  fMapClustertoPtR2(),
58  fMapClustertoPtR3(),
59  fMapClustertoPtR4(),
60  fListHistos(NULL),
61  fHistTest(NULL),
62  fHistClusterEnergy(NULL),
63  fHistIso(NULL)
64 {
65  // Default constructor
66  DefineInput(0, TChain::Class());
67 }
68 
69 
70 //________________________________________________________________________
72  // default deconstructor
73  fMapClustertoPtR1.clear();
74  fMapClustertoPtR2.clear();
75  fMapClustertoPtR3.clear();
76  fMapClustertoPtR4.clear();
77 
78  if(fHistIso) delete fHistIso;
79  if(fHistTest) delete fHistTest;
80  if(fHistClusterEnergy) delete fHistClusterEnergy;
81  if(fListHistos != NULL){
82  delete fListHistos;
83  fListHistos = NULL; //#################################################################### Why?
84  }
85 }
86 
87 //________________________________________________________________________
89  fMapClustertoPtR1.clear();
90  fMapClustertoPtR2.clear();
91  fMapClustertoPtR3.clear();
92  fMapClustertoPtR4.clear();
93 }
94 
95 //________________________________________________________________________
96 //@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@ OUTPUT
98  if(fListHistos != NULL){
99  delete fListHistos;
100  fListHistos = NULL;
101  }
102  if(fListHistos == NULL){
103  fListHistos = new TList();
104  fListHistos->SetOwner(kTRUE);
105  fListHistos->SetName(Form("PhotonIsolation_%i",fPhotonType));
106  }
107 
108  //Create user output objects
109  fHistIso = new TH1F("fHistM02","M02 Distribution",200,0,50);
110  fHistTest = new TH1F("fHistClus","Number of clusters per Event",20,0,20);
111  fHistClusterEnergy = new TH1F("fHistEnergy","Pt Spectrum",280,0,140);
112  fListHistos->Add(fHistTest);
113  fListHistos->Add(fHistIso);
114  fListHistos->Add(fHistClusterEnergy);
115 
116 }
117 
118 //________________________________________________________________________
120  // Initialize function to be called once before analysis
121  fMapClustertoPtR1.clear();
122  fMapClustertoPtR2.clear();
123  fMapClustertoPtR3.clear();
124  fMapClustertoPtR4.clear();
125 }
126 
127 //________________________________________________________________________
128 //@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@ USEREXEC
130 
131  // main method of AliPhotonIsolation, first initialize and then process event
132 
133  Initialize();
134 
135  ProcessEvent(fInputEvent);
136 
137  //DebugIsolation();
138 
139  return;
140 }
141 
142 //________________________________________________________________________
143 //@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@ PROCESS EVENT
144 void AliPhotonIsolation::ProcessEvent(AliVEvent *event){
145  Bool_t debug = 0; //Debugging on/off
146  //###########################################################GET NUMBER OF CLUSTERS IN EVENT
147  Int_t nClus = 0;
148  TClonesArray * arrClusters = NULL;
149  if(!fCorrTaskSetting.CompareTo("")){
150  nClus = event->GetNumberOfCaloClusters();
151  } else {
152  arrClusters = dynamic_cast<TClonesArray*>(event->FindListObject(Form("%sClustersBranch",fCorrTaskSetting.Data())));
153  nClus = arrClusters->GetEntries();
154  }
155 
156  //#####################################################CREATE CLASS TO GET EVENT INFORMATION
157  AliESDEvent *esdev = dynamic_cast<AliESDEvent*>(event);
158  AliAODEvent *aodev = 0;
159  if (!esdev) {
160  aodev = dynamic_cast<AliAODEvent*>(event);
161  if (!aodev) {
162  AliError("Task needs AOD or ESD event, returning");
163  return;
164  }
165  }
166 
167  if(debug){
168  cout << "___________________________________________________This is the start of an event___________________________________" << endl;
169  }
170 
171  fHistTest->Fill(nClus);
172 
173  //#########################################################################LOOP OVER CLUSTERS
174  for(Int_t iclus=0;iclus < nClus;iclus++){
175 
176  if(debug){
177  cout << "________________________________̣̣This is a new cluster_________________________________" << endl;
178  }
179 
180  Float_t ptsum1 = 0., ptsum2 = 0., ptsum3 = 0., ptsum4 = 0.;
181  //#####################################################INITIALISE AND CREATE CLUSTER OBJECT
182  AliVCluster* cluster = NULL;
183  if(arrClusters){
184  if(esdev){
185  if(arrClusters)
186  cluster = new AliESDCaloCluster(*(AliESDCaloCluster*)arrClusters->At(iclus));
187  } else if(aodev){
188  if(arrClusters)
189  cluster = new AliAODCaloCluster(*(AliAODCaloCluster*)arrClusters->At(iclus));
190  }
191  }
192  else
193  cluster = event->GetCaloCluster(iclus); //create ALiVCluster object
194  if (!cluster){ continue;}
195  //####################################################DETERMINE ETA, PHI AND E_T FOR THE CLUSTER
196  Float_t eta_clus = 0., phi_clus = 0., theta_clus = 0; //define variables for cluster properties
197  Float_t clspos[3] = {0.,0.,0.}; //initialise position vector
198 
199  cluster->GetPosition(clspos); //get position of cluster
200  TVector3 clspos_v3(clspos); //create 3Vector to gain access to geometrical functions
201  eta_clus = clspos_v3.Eta(); //determine Pseudo rapidity
202  phi_clus = clspos_v3.Phi(); //determine Phi
203  theta_clus = clspos_v3.Theta(); //determine Theta in radians
204 
205  Float_t ET_clus = (cluster->E())*TMath::Sin(theta_clus); //Determine cluster energy (E_T=E*sin(theta), because m=0 or negligible)
206 
207  Float_t m02 = cluster->GetM02();
208  // cout << "M02: " << m02 << endl;
209 
210  fHistIso->Fill(m02);
211 
212  if(debug){
213  cout << "Clusterenergie: " << ET_clus << endl;
214  }
215 
216  fHistClusterEnergy->Fill(ET_clus);
217 
218  //#########################################################################LOOP OVER TRACKS
219  for (Int_t itr=0;itr<event->GetNumberOfTracks();itr++){
220  if(debug){
221  cout << "______________________________This is a new track_______________________" << endl;
222  }
223 
224  //#####################################################INITIALISE AND CREATE TRACK OBJECT
225  AliVTrack *inTrack = 0x0; //initialise general track variable
226  if(esdev){
227  inTrack = esdev->GetTrack(itr); //get track
228  if(!inTrack) continue; //test track
229  // AliESDtrack *esdt = dynamic_cast<AliESDtrack*>(inTrack); //create class with track information
230  } else if(aodev) { //same for AODs
231  inTrack = dynamic_cast<AliVTrack*>(aodev->GetTrack(itr));
232  if(!inTrack) continue;
233  // AliAODTrack *aodt = dynamic_cast<AliAODTrack*>(inTrack);
234  }
235 
236  //####################################################DETERMINE ETA AND PHI FOR THE TRACK
237  Float_t eta_track = 0., phi_track = 0.; //define variables for track properties
238  eta_track = inTrack->Eta();
239  phi_track = inTrack->Phi();
240 
241  //#####################################CHECK THE DISTANCE/ANGLE BETWEEN TRACK AND CLUSTER
242  Float_t R2 = (TMath::Power(eta_clus-eta_track,2)+TMath::Power(phi_clus-phi_track,2)); //calculate the squared radius
243 /*
244  if(debug){
245  cout << "PhiTrack: " << phi_track << endl;
246  cout << "EtaTrack: " << eta_track << endl;
247  cout << "PhiClus: " << phi_clus << endl;
248  cout << "EtaClus: " << eta_clus << endl;
249  cout << "R: " << R2 << endl;
250  }
251 */
252 /*
253  if(debug){
254  cout << << endl;
255  }
256 */
257 
258  Float_t pT = inTrack->Pt();//get pT of track
259 
260  //add up pT values of tracks in cone
261  if(R2<0.01){//compare to R=0.1
262  ptsum1 += pT;
263  }else if(R2<0.04){//compare to R=0.2
264  ptsum2 += pT;
265  }else if(R2<0.09){//compare to R=0.3
266  ptsum3 += pT;
267  }else if(R2<0.16){//compare to R=0.4
268  ptsum4 += pT;
269  }
270 
271  if(debug && R2<0.16){
272  cout << "R2: " << R2 << endl;
273  }
274 
275  } //end of for loop tracks
276 
277  ptsum2+=ptsum1; //add pT of innerst cone to second cone
278  ptsum3+=ptsum2; //add complete second to third
279  ptsum4+=ptsum3; //add complete third to fourth
280 
281  if(debug){
282  cout << "ptsum1: " << ptsum1 << endl;
283  cout << "ptsum2: " << ptsum2 << endl;
284  cout << "ptsum3: " << ptsum3 << endl;
285  cout << "ptsum4: " << ptsum4 << endl;
286  }
287 
288  //map Ptsum to clusterID, map all the ptsum that are bigger than 0
289  if(ptsum1>0){
290  fMapClustertoPtR1.insert(make_pair(cluster->GetID(),ptsum1));
291  fMapClustertoPtR2.insert(make_pair(cluster->GetID(),ptsum2));
292  fMapClustertoPtR3.insert(make_pair(cluster->GetID(),ptsum3));
293  fMapClustertoPtR4.insert(make_pair(cluster->GetID(),ptsum4));
294  }else if(ptsum2>0){
295  fMapClustertoPtR2.insert(make_pair(cluster->GetID(),ptsum2));
296  fMapClustertoPtR3.insert(make_pair(cluster->GetID(),ptsum3));
297  fMapClustertoPtR4.insert(make_pair(cluster->GetID(),ptsum4));
298  }else if(ptsum3>0){
299  fMapClustertoPtR3.insert(make_pair(cluster->GetID(),ptsum3));
300  fMapClustertoPtR4.insert(make_pair(cluster->GetID(),ptsum4));
301  }else if(ptsum4>0){
302  fMapClustertoPtR4.insert(make_pair(cluster->GetID(),ptsum4));
303  }
304 
305  if(debug){
306  cout << "Cluster ID: " << iclus << endl;
307  cout << "ptsumR4: " << ptsum4 << endl;
308  }
309 
310  }//end for loop clusters
311 
312  if(debug){
313  for (map<Int_t,Float_t>::iterator it= fMapClustertoPtR4.begin(); it != fMapClustertoPtR4.end(); ++it) {
314  cout << it->first << "\t" << it->second << endl;
315  }
316  }
317 
318  return;
319 }
320 
321 
322 //Determine Isolation dependent on ClusterID, Isolation Cone and Isolation Energy
323 
325 
326  Bool_t isolated = kFALSE;
327  if (R>0.09 && R<0.11){ //check Radius
328  if(fMapClustertoPtR1[clusterID]<isoPt) isolated = kTRUE;
329  //cout << "MapPt" << fMapClustertoPtR1[clusterID] << endl;
330  }else if(R>0.19 && R<0.21){ //check Radius
331  if(fMapClustertoPtR2[clusterID]<isoPt) isolated = kTRUE;
332  }else if(R>0.29 && R<0.31){ //check Radius
333  if(fMapClustertoPtR3[clusterID]<isoPt) isolated = kTRUE;
334  }else if(R>0.39 && R<0.41){ //check Radius
335  if(fMapClustertoPtR4[clusterID]<isoPt) isolated = kTRUE;
336  }else{
337  AliFatal("AliPhotonIsolation: GetIsolation - Radius of the Isolation Cone has to be set to R=0.1, R=0.2, R=0.3 or R=0.4");
338  }
339 
340  // cout << "_______________________________________________________________________Here I am to check out isolated Photons!" << endl;
341 
342  return isolated;
343 }
344 
345 
347 
348  return;
349 }
350 
351 
352 //
virtual void Terminate(Option_t *)
Bool_t GetIsolation(Int_t clusterID, Float_t R, Float_t isoPt)
void ProcessEvent(AliVEvent *event)
int Int_t
Definition: External.C:63
float Float_t
Definition: External.C:68
virtual void UserExec(Option_t *option)
const char Option_t
Definition: External.C:48
bool Bool_t
Definition: External.C:53