1 /**************************************************************************
2 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
4 * Author: The ALICE Off-line Project. *
5 * Contributors are mentioned in the code where appropriate. *
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 **************************************************************************/
18 //_________________________________________________________________________
19 // A basic analysis task to analyse photon detected by PHOS
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
26 //////////////////////////////////////////////////////////////////////////////
37 #include "AliAnaGammaPhos.h"
38 #include "AliAnalysisManager.h"
39 #include "AliESDEvent.h"
40 #include "AliESDCaloCluster.h"
41 #include "AliAODEvent.h"
42 #include "AliAODHandler.h"
45 //______________________________________________________________________________
46 AliAnaGammaPhos::AliAnaGammaPhos() :
60 fhPHOSRecParticles(0),
62 fhPHOSInvariantMass(0),
67 //______________________________________________________________________________
68 AliAnaGammaPhos::AliAnaGammaPhos(const char *name) :
69 AliAnalysisTask(name,""),
83 fhPHOSRecParticles(0),
85 fhPHOSInvariantMass(0),
90 DefineInput(0, TChain::Class());
92 DefineOutput(0, TTree::Class()) ;
93 DefineOutput(1, TList::Class()) ;
96 //____________________________________________________________________________
97 AliAnaGammaPhos::AliAnaGammaPhos(const AliAnaGammaPhos& ap) :
98 AliAnalysisTask(ap.GetName(),""),
103 fAODPhotons(ap.fAODPhotons),
104 fPhotonsInPhos(ap.fPhotonsInPhos),
106 fPhotonId(ap.fPhotonId),
107 fOutputList(ap.fOutputList),
108 fhPHOSPos(ap.fhPHOSPos),
110 fhPHOSEnergy(ap.fhPHOSEnergy),
111 fhPHOSDigits(ap.fhPHOSDigits),
112 fhPHOSRecParticles(ap.fhPHOSRecParticles),
113 fhPHOSPhotons(ap.fhPHOSPhotons),
114 fhPHOSInvariantMass(ap.fhPHOSInvariantMass),
115 fhPHOSDigitsEvent(ap.fhPHOSDigitsEvent)
120 //_____________________________________________________________________________
121 AliAnaGammaPhos& AliAnaGammaPhos::operator = (const AliAnaGammaPhos& ap)
123 // assignment operator
125 this->~AliAnaGammaPhos();
126 new(this) AliAnaGammaPhos(ap);
130 //______________________________________________________________________________
131 AliAnaGammaPhos::~AliAnaGammaPhos()
134 // fOutputList->Clear() ;
135 //delete fOutputList ;
139 //______________________________________________________________________________
140 void AliAnaGammaPhos::ConnectInputData(const Option_t*)
142 // Initialisation of branch container and histograms
144 AliInfo(Form("*** Initialization of %s", GetName())) ;
147 fChain = dynamic_cast<TChain *>(GetInputData(0)) ;
149 AliError(Form("Input 0 for %s not found\n", GetName()));
153 fESD = new AliESDEvent();
154 fESD->ReadFromTree(fChain);
158 //________________________________________________________________________
159 void AliAnaGammaPhos::CreateOutputObjects()
161 // Create the outputs containers
164 AliAODHandler* handler = (AliAODHandler*) ((AliAnalysisManager::GetAnalysisManager())->GetOutputEventHandler());
165 fTreeA = handler->GetTree() ;
166 fAOD = handler->GetAOD();
167 fAODPhotons = fAOD->GetCaloClusters() ;
172 fhPHOSPos = new TNtuple("PHOSPos" , "Position in PHOS" , "x:y:z");
173 fhPHOS = new TNtuple("PHOS" , "PHOS" , "event:digits:clusters:photons");
174 fhPHOSEnergy = new TH1D("PHOSEnergy" , "PHOSEnergy" , 100, 0., 100. ) ;
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 ) ;
181 // create output container
183 fOutputList = new TList() ;
184 fOutputList->SetName(GetName()) ;
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) ;
196 //______________________________________________________________________________
197 void AliAnaGammaPhos::Exec(Option_t *)
199 // Processing of one event
201 Long64_t entry = fChain->GetReadEntry() ;
204 AliError("fESD is not connected to the input!") ;
208 // if ( !((entry-1)%100) )
209 AliInfo(Form("%s ----> Processing event # %lld", (dynamic_cast<TChain *>(fChain))->GetFile()->GetName(), entry)) ;
211 //************************ PHOS *************************************
212 TRefArray * caloClustersArr = new TRefArray();
213 fESD->GetPHOSClusters(caloClustersArr);
215 const Int_t kNumberOfPhosClusters = caloClustersArr->GetEntries() ;
217 TVector3 ** phosVector = new TVector3*[kNumberOfPhosClusters] ;
218 Float_t * phosPhotonsEnergy = new Float_t[kNumberOfPhosClusters] ;
220 Int_t numberOfDigitsInPhos = 0 ;
223 // loop over the PHOS Cluster
224 for(phosCluster = 0 ; phosCluster < kNumberOfPhosClusters ; phosCluster++) {
225 AliESDCaloCluster * caloCluster = (AliESDCaloCluster *) caloClustersArr->At(phosCluster) ;
226 //AliESDCaloCluster * caloCluster = fESD->GetCaloCluster(phosCluster) ;
229 caloCluster->GetPosition( pos ) ;
230 fhPHOSEnergy->Fill( caloCluster->E() ) ;
231 fhPHOSPos->Fill( pos[0], pos[1], pos[2] ) ;
232 fhPHOSDigits->Fill(entry, caloCluster->GetNCells() ) ;
233 numberOfDigitsInPhos += caloCluster->GetNCells() ;
234 Double_t * pid = caloCluster->GetPid() ;
235 if(pid[AliPID::kPhoton] > GetPhotonId() ) {
236 phosVector[fPhotonsInPhos] = new TVector3(pos[0],pos[1],pos[2]) ;
237 phosPhotonsEnergy[fPhotonsInPhos]=caloCluster->E() ;
238 //new ((*fAODPhotons)[fPhotonsInPhos++;]) AliAODPhoton ( );
243 fhPHOSRecParticles->Fill(kNumberOfPhosClusters);
244 fhPHOSPhotons->Fill(fPhotonsInPhos);
245 fhPHOSDigitsEvent->Fill(numberOfDigitsInPhos);
246 fhPHOS->Fill(entry, numberOfDigitsInPhos, kNumberOfPhosClusters, fPhotonsInPhos) ;
249 if (fPhotonsInPhos > 1 ) {
250 Int_t phosPhoton1, phosPhoton2 ;
251 for(phosPhoton1 = 0 ; phosPhoton1 < fPhotonsInPhos ; phosPhoton1++) {
252 for(phosPhoton2 = phosPhoton1 + 1 ; phosPhoton2 < fPhotonsInPhos ; phosPhoton2++) {
253 Float_t tempMass = TMath::Sqrt( 2 * phosPhotonsEnergy[phosPhoton1] * phosPhotonsEnergy[phosPhoton2] *
254 ( 1 - TMath::Cos(phosVector[phosPhoton1]->Angle(*phosVector[phosPhoton2])) )
256 fhPHOSInvariantMass->Fill(tempMass*1000.);
261 PostData(0, fTreeA) ;
262 PostData(1, fOutputList);
264 delete [] phosVector ;
265 delete [] phosPhotonsEnergy ;
270 //______________________________________________________________________________
271 void AliAnaGammaPhos::Init()
273 // Intialisation of parameters
274 AliInfo("Doing initialisation") ;
278 //______________________________________________________________________________
279 void AliAnaGammaPhos::Terminate(Option_t *)
281 // Processing when the event loop is ended
283 // fOutputList = (TList*)GetOutputData(0);
284 // fhPHOSPos = (TNtuple*)fOutputList->At(0);
285 // fhPHOS = (TNtuple*)fOutputList->At(1);
286 // fhPHOSEnergy = (TH1D*)fOutputList->At(2);
287 // fhPHOSDigits = (TH1I*)fOutputList->At(3);
288 // fhPHOSRecParticles = (TH1D*)fOutputList->At(4);
289 // fhPHOSPhotons = (TH1I*)fOutputList->At(5);
290 // fhPHOSInvariantMass = (TH1D*)fOutputList->At(6);
291 // fhPHOSDigitsEvent = (TH1I*)fOutputList->At(7);
294 Bool_t problem = kFALSE ;
295 AliInfo(Form(" *** %s Report:", GetName())) ;
296 printf(" PHOSEnergy Mean : %5.3f , RMS : %5.3f \n", fhPHOSEnergy->GetMean(), fhPHOSEnergy->GetRMS() ) ;
297 printf(" PHOSDigits Mean : %5.3f , RMS : %5.3f \n", fhPHOSDigits->GetMean(), fhPHOSDigits->GetRMS() ) ;
298 printf(" PHOSRecParticles Mean : %5.3f , RMS : %5.3f \n", fhPHOSRecParticles->GetMean(), fhPHOSRecParticles->GetRMS() ) ;
299 printf(" PHOSPhotons Mean : %5.3f , RMS : %5.3f \n", fhPHOSPhotons->GetMean(), fhPHOSPhotons->GetRMS() ) ;
300 printf(" PHOSInvariantMass Mean : %5.3f , RMS : %5.3f \n", fhPHOSInvariantMass->GetMean(), fhPHOSInvariantMass->GetRMS() ) ;
301 printf(" PHOSDigitsEvent Mean : %5.3f , RMS : %5.3f \n", fhPHOSDigitsEvent->GetMean(), fhPHOSDigitsEvent->GetRMS() ) ;
303 TCanvas * cPHOS = new TCanvas("cPHOS", "PHOS ESD Test", 400, 10, 600, 700) ;
307 if ( fhPHOSEnergy->GetMaximum() > 0. )
309 fhPHOSEnergy->SetAxisRange(0, 25.);
310 fhPHOSEnergy->SetLineColor(2);
311 fhPHOSEnergy->Draw();
314 fhPHOSDigits->SetAxisRange(0,25.);
315 fhPHOSDigits->SetLineColor(2);
316 fhPHOSDigits->Draw();
319 if ( fhPHOSRecParticles->GetMaximum() > 0. )
321 fhPHOSRecParticles->SetAxisRange(0, 25.);
322 fhPHOSRecParticles->SetLineColor(2);
323 fhPHOSRecParticles->Draw();
326 if ( fhPHOSPhotons->GetMaximum() > 0. )
328 fhPHOSPhotons->SetAxisRange(0,25.);
329 fhPHOSPhotons->SetLineColor(2);
330 fhPHOSPhotons->Draw();
333 fhPHOSInvariantMass->SetLineColor(2);
334 fhPHOSInvariantMass->Draw();
337 if ( fhPHOSDigitsEvent->GetMaximum() > 0. )
339 fhPHOSDigitsEvent->SetAxisRange(0,40.);
340 fhPHOSDigitsEvent->SetLineColor(2);
341 fhPHOSDigitsEvent->Draw();
343 cPHOS->Print("PHOS.eps");
346 sprintf(line, ".!tar -zcf %s.tar.gz *.eps", GetName()) ;
347 gROOT->ProcessLine(line);
348 sprintf(line, ".!rm -fR *.eps");
349 gROOT->ProcessLine(line);
351 AliInfo(Form("!!! All the eps files are in %s.tar.gz !!!", GetName())) ;
355 report="Problems found, please check!!!";
359 AliInfo(Form("*** %s Summary Report: %s \n",GetName(), report)) ;