Coding violations corrected
[u/mrichter/AliRoot.git] / PWG4 / AliAnaGammaPhos.cxx
CommitLineData
ce9d0223 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
16/* $Id: */
17
18//_________________________________________________________________________
19// A basic analysis task to analyse photon detected by PHOS
d40a7452 20// A basic analysis task to analyse photon detected by PHOS
21// A basic analysis task to analyse photon detected by PHOS
22// A basic analysis task to analyse photon detected by PHOS
23// A basic analysis task to analyse photon detected by PHOS
ce9d0223 24//
25//*-- Yves Schutz
26//////////////////////////////////////////////////////////////////////////////
27
28#include <TCanvas.h>
29#include <TChain.h>
30#include <TFile.h>
31#include <TH1.h>
ce9d0223 32#include <TLegend.h>
33#include <TNtuple.h>
34#include <TROOT.h>
35#include <TVector3.h>
36
37#include "AliAnaGammaPhos.h"
38#include "AliAnalysisManager.h"
9ea45617 39#include "AliESDEvent.h"
40#include "AliESDCaloCluster.h"
ce9d0223 41#include "AliAODEvent.h"
42#include "AliAODHandler.h"
43#include "AliLog.h"
44
45//______________________________________________________________________________
46AliAnaGammaPhos::AliAnaGammaPhos() :
47 fChain(0x0),
48 fDebug(0),
49 fESD(0x0),
50 fAOD(0x0),
51 fAODPhotons(0x0),
52 fPhotonsInPhos(0),
53 fTreeA(0x0),
54 fPhotonId(1.0),
55 fOutputList(0x0),
3bb2c538 56 fhPHOSPos(0),
ce9d0223 57 fhPHOS(0),
58 fhPHOSEnergy(0),
59 fhPHOSDigits(0),
60 fhPHOSRecParticles(0),
61 fhPHOSPhotons(0),
62 fhPHOSInvariantMass(0),
63 fhPHOSDigitsEvent(0)
64{
65 //Default constructor
66}
67//______________________________________________________________________________
68AliAnaGammaPhos::AliAnaGammaPhos(const char *name) :
69 AliAnalysisTask(name,""),
70 fChain(0x0),
71 fDebug(0),
72 fESD(0x0),
73 fAOD(0x0),
74 fAODPhotons(0x0),
75 fPhotonsInPhos(0),
76 fTreeA(0x0),
77 fPhotonId(1.0),
78 fOutputList(0x0),
3bb2c538 79 fhPHOSPos(0),
ce9d0223 80 fhPHOS(0),
81 fhPHOSEnergy(0),
82 fhPHOSDigits(0),
83 fhPHOSRecParticles(0),
84 fhPHOSPhotons(0),
85 fhPHOSInvariantMass(0),
86 fhPHOSDigitsEvent(0)
87{
88 // Constructor.
89 // Input slot #0
90 DefineInput(0, TChain::Class());
91 // Output slots
92 DefineOutput(0, TTree::Class()) ;
93 DefineOutput(1, TList::Class()) ;
94}
95
d40a7452 96//____________________________________________________________________________
97AliAnaGammaPhos::AliAnaGammaPhos(const AliAnaGammaPhos& ap) :
98 AliAnalysisTask(ap.GetName(),""),
99 fChain(ap.fChain),
100 fDebug(ap.fDebug),
101 fESD(ap.fESD),
102 fAOD(ap.fAOD),
103 fAODPhotons(ap.fAODPhotons),
104 fPhotonsInPhos(ap.fPhotonsInPhos),
105 fTreeA(ap.fTreeA),
106 fPhotonId(ap.fPhotonId),
107 fOutputList(ap.fOutputList),
108 fhPHOSPos(ap.fhPHOSPos),
109 fhPHOS(ap.fhPHOS),
110 fhPHOSEnergy(ap.fhPHOSEnergy),
111 fhPHOSDigits(ap.fhPHOSDigits),
112 fhPHOSRecParticles(ap.fhPHOSRecParticles),
113 fhPHOSPhotons(ap.fhPHOSPhotons),
114 fhPHOSInvariantMass(ap.fhPHOSInvariantMass),
115 fhPHOSDigitsEvent(ap.fhPHOSDigitsEvent)
116{
117 // cpy ctor
118}
119
120//_____________________________________________________________________________
121AliAnaGammaPhos& AliAnaGammaPhos::operator = (const AliAnaGammaPhos& ap)
122{
123// assignment operator
124
125 this->~AliAnaGammaPhos();
126 new(this) AliAnaGammaPhos(ap);
127 return *this;
128}
129
ce9d0223 130//______________________________________________________________________________
131AliAnaGammaPhos::~AliAnaGammaPhos()
132{
133 // dtor
c0acc236 134 // fOutputList->Clear() ;
135 //delete fOutputList ;
ce9d0223 136}
137
138
139//______________________________________________________________________________
140void AliAnaGammaPhos::ConnectInputData(const Option_t*)
141{
142 // Initialisation of branch container and histograms
143
144 AliInfo(Form("*** Initialization of %s", GetName())) ;
145
146 // Get input data
147 fChain = dynamic_cast<TChain *>(GetInputData(0)) ;
148 if (!fChain) {
149 AliError(Form("Input 0 for %s not found\n", GetName()));
150 return ;
151 }
152
9ea45617 153 fESD = new AliESDEvent();
154 fESD->ReadFromTree(fChain);
155
ce9d0223 156}
157
158//________________________________________________________________________
159void AliAnaGammaPhos::CreateOutputObjects()
160{
161 // Create the outputs containers
162
163 OpenFile(0) ;
691685d6 164 AliAODHandler* handler = (AliAODHandler*) ((AliAnalysisManager::GetAnalysisManager())->GetOutputEventHandler());
c0acc236 165 fTreeA = handler->GetTree() ;
ce9d0223 166 fAOD = handler->GetAOD();
4b707925 167 fAODPhotons = fAOD->GetCaloClusters() ;
c0acc236 168
ce9d0223 169
170 OpenFile(1) ;
171
172 fhPHOSPos = new TNtuple("PHOSPos" , "Position in PHOS" , "x:y:z");
173 fhPHOS = new TNtuple("PHOS" , "PHOS" , "event:digits:clusters:photons");
bdcfac30 174 fhPHOSEnergy = new TH1D("PHOSEnergy" , "PHOSEnergy" , 100, 0., 100. ) ;
ce9d0223 175 fhPHOSDigits = new TH1I("PHOSDigitsCluster" , "PHOSDigits" , 20 , 0 , 20 ) ;
176 fhPHOSRecParticles = new TH1D("PHOSRecParticles" , "PHOSRecParticles" , 20 , 0., 20. ) ;
177 fhPHOSPhotons = new TH1I("PHOSPhotons" , "PHOSPhotons" , 20 , 0 , 20 ) ;
178 fhPHOSInvariantMass = new TH1D("PHOSInvariantMass" , "PHOSInvariantMass" , 400, 0., 400.) ;
179 fhPHOSDigitsEvent = new TH1I("PHOSDigitsEvent" , "PHOSDigitsEvent" , 30 , 0 , 30 ) ;
180
181 // create output container
182
183 fOutputList = new TList() ;
184 fOutputList->SetName(GetName()) ;
185
186 fOutputList->AddAt(fhPHOSPos, 0) ;
187 fOutputList->AddAt(fhPHOS, 1) ;
188 fOutputList->AddAt(fhPHOSEnergy, 2) ;
189 fOutputList->AddAt(fhPHOSDigits, 3) ;
190 fOutputList->AddAt(fhPHOSRecParticles, 4) ;
191 fOutputList->AddAt(fhPHOSPhotons, 5) ;
192 fOutputList->AddAt(fhPHOSInvariantMass, 6) ;
193 fOutputList->AddAt(fhPHOSDigitsEvent, 7) ;
194}
195
196//______________________________________________________________________________
197void AliAnaGammaPhos::Exec(Option_t *)
198{
199 // Processing of one event
200
201 Long64_t entry = fChain->GetReadEntry() ;
202
203 if (!fESD) {
204 AliError("fESD is not connected to the input!") ;
205 return ;
206 }
207
208 if ( !((entry-1)%100) )
209 AliInfo(Form("%s ----> Processing event # %lld", (dynamic_cast<TChain *>(fChain))->GetFile()->GetName(), entry)) ;
210
211 //************************ PHOS *************************************
212
213 Int_t firstPhosCluster = fESD->GetFirstPHOSCluster() ;
d40a7452 214 const Int_t kNumberOfPhosClusters = fESD->GetNumberOfPHOSClusters() ;
ce9d0223 215
d40a7452 216 TVector3 ** phosVector = new TVector3*[kNumberOfPhosClusters] ;
217 Float_t * phosPhotonsEnergy = new Float_t[kNumberOfPhosClusters] ;
ce9d0223 218 Int_t phosCluster ;
219 Int_t numberOfDigitsInPhos = 0 ;
220
221 fPhotonsInPhos = 0 ;
222 // loop over the PHOS Cluster
d40a7452 223 for(phosCluster = firstPhosCluster ; phosCluster < firstPhosCluster + kNumberOfPhosClusters ; phosCluster++) {
ce9d0223 224 AliESDCaloCluster * caloCluster = fESD->GetCaloCluster(phosCluster) ;
225 if (caloCluster) {
226 Float_t pos[3] ;
c0acc236 227 caloCluster->GetPosition( pos ) ;
228 fhPHOSEnergy->Fill( caloCluster->E() ) ;
ce9d0223 229 fhPHOSPos->Fill( pos[0], pos[1], pos[2] ) ;
230 fhPHOSDigits->Fill(entry, caloCluster->GetNumberOfDigits() ) ;
231 numberOfDigitsInPhos += caloCluster->GetNumberOfDigits() ;
4b707925 232 Double_t * pid = caloCluster->GetPid() ;
ce9d0223 233 if(pid[AliPID::kPhoton] > GetPhotonId() ) {
234 phosVector[fPhotonsInPhos] = new TVector3(pos[0],pos[1],pos[2]) ;
c0acc236 235 phosPhotonsEnergy[fPhotonsInPhos]=caloCluster->E() ;
ce9d0223 236 //new ((*fAODPhotons)[fPhotonsInPhos++;]) AliAODPhoton ( );
237 }
238 }
239 } //PHOS clusters
240
d40a7452 241 fhPHOSRecParticles->Fill(kNumberOfPhosClusters);
ce9d0223 242 fhPHOSPhotons->Fill(fPhotonsInPhos);
243 fhPHOSDigitsEvent->Fill(numberOfDigitsInPhos);
d40a7452 244 fhPHOS->Fill(entry, numberOfDigitsInPhos, kNumberOfPhosClusters, fPhotonsInPhos) ;
ce9d0223 245
246 // invariant Mass
247 if (fPhotonsInPhos > 1 ) {
248 Int_t phosPhoton1, phosPhoton2 ;
249 for(phosPhoton1 = 0 ; phosPhoton1 < fPhotonsInPhos ; phosPhoton1++) {
250 for(phosPhoton2 = phosPhoton1 + 1 ; phosPhoton2 < fPhotonsInPhos ; phosPhoton2++) {
251 Float_t tempMass = TMath::Sqrt( 2 * phosPhotonsEnergy[phosPhoton1] * phosPhotonsEnergy[phosPhoton2] *
252 ( 1 - TMath::Cos(phosVector[phosPhoton1]->Angle(*phosVector[phosPhoton2])) )
253 );
254 fhPHOSInvariantMass->Fill(tempMass*1000.);
255 }
256 }
257 }
32392f88 258
ce9d0223 259 PostData(0, fTreeA) ;
260 PostData(1, fOutputList);
261
262 delete [] phosVector ;
263 delete [] phosPhotonsEnergy ;
264
265}
266
267
268//______________________________________________________________________________
269void AliAnaGammaPhos::Init()
270{
271 // Intialisation of parameters
272 AliInfo("Doing initialisation") ;
273 SetPhotonId(0.9) ;
274}
275
276//______________________________________________________________________________
277void AliAnaGammaPhos::Terminate(Option_t *)
278{
279 // Processing when the event loop is ended
280
c0acc236 281// fOutputList = (TList*)GetOutputData(0);
282// fhPHOSPos = (TNtuple*)fOutputList->At(0);
283// fhPHOS = (TNtuple*)fOutputList->At(1);
284// fhPHOSEnergy = (TH1D*)fOutputList->At(2);
285// fhPHOSDigits = (TH1I*)fOutputList->At(3);
286// fhPHOSRecParticles = (TH1D*)fOutputList->At(4);
287// fhPHOSPhotons = (TH1I*)fOutputList->At(5);
288// fhPHOSInvariantMass = (TH1D*)fOutputList->At(6);
289// fhPHOSDigitsEvent = (TH1I*)fOutputList->At(7);
ce9d0223 290
32392f88 291
ce9d0223 292 Bool_t problem = kFALSE ;
293 AliInfo(Form(" *** %s Report:", GetName())) ;
294 printf(" PHOSEnergy Mean : %5.3f , RMS : %5.3f \n", fhPHOSEnergy->GetMean(), fhPHOSEnergy->GetRMS() ) ;
295 printf(" PHOSDigits Mean : %5.3f , RMS : %5.3f \n", fhPHOSDigits->GetMean(), fhPHOSDigits->GetRMS() ) ;
296 printf(" PHOSRecParticles Mean : %5.3f , RMS : %5.3f \n", fhPHOSRecParticles->GetMean(), fhPHOSRecParticles->GetRMS() ) ;
297 printf(" PHOSPhotons Mean : %5.3f , RMS : %5.3f \n", fhPHOSPhotons->GetMean(), fhPHOSPhotons->GetRMS() ) ;
298 printf(" PHOSInvariantMass Mean : %5.3f , RMS : %5.3f \n", fhPHOSInvariantMass->GetMean(), fhPHOSInvariantMass->GetRMS() ) ;
299 printf(" PHOSDigitsEvent Mean : %5.3f , RMS : %5.3f \n", fhPHOSDigitsEvent->GetMean(), fhPHOSDigitsEvent->GetRMS() ) ;
300
301 TCanvas * cPHOS = new TCanvas("cPHOS", "PHOS ESD Test", 400, 10, 600, 700) ;
302 cPHOS->Divide(3, 2);
303
304 cPHOS->cd(1) ;
305 if ( fhPHOSEnergy->GetMaximum() > 0. )
306 gPad->SetLogy();
307 fhPHOSEnergy->SetAxisRange(0, 25.);
308 fhPHOSEnergy->SetLineColor(2);
309 fhPHOSEnergy->Draw();
310
311 cPHOS->cd(2) ;
312 fhPHOSDigits->SetAxisRange(0,25.);
313 fhPHOSDigits->SetLineColor(2);
314 fhPHOSDigits->Draw();
315
316 cPHOS->cd(3) ;
317 if ( fhPHOSRecParticles->GetMaximum() > 0. )
318 gPad->SetLogy();
319 fhPHOSRecParticles->SetAxisRange(0, 25.);
320 fhPHOSRecParticles->SetLineColor(2);
321 fhPHOSRecParticles->Draw();
322
323 cPHOS->cd(4) ;
324 if ( fhPHOSPhotons->GetMaximum() > 0. )
325 gPad->SetLogy();
326 fhPHOSPhotons->SetAxisRange(0,25.);
327 fhPHOSPhotons->SetLineColor(2);
328 fhPHOSPhotons->Draw();
329
330 cPHOS->cd(5) ;
331 fhPHOSInvariantMass->SetLineColor(2);
332 fhPHOSInvariantMass->Draw();
333
334 cPHOS->cd(6) ;
335 if ( fhPHOSDigitsEvent->GetMaximum() > 0. )
336 gPad->SetLogy();
337 fhPHOSDigitsEvent->SetAxisRange(0,40.);
338 fhPHOSDigitsEvent->SetLineColor(2);
339 fhPHOSDigitsEvent->Draw();
340
341 cPHOS->Print("PHOS.eps");
342
343 char line[1024] ;
344 sprintf(line, ".!tar -zcf %s.tar.gz *.eps", GetName()) ;
345 gROOT->ProcessLine(line);
346 sprintf(line, ".!rm -fR *.eps");
347 gROOT->ProcessLine(line);
348
349 AliInfo(Form("!!! All the eps files are in %s.tar.gz !!!", GetName())) ;
350
351 char * report ;
352 if(problem)
353 report="Problems found, please check!!!";
354 else
355 report="OK";
356
357 AliInfo(Form("*** %s Summary Report: %s \n",GetName(), report)) ;
358}