]> git.uio.no Git - u/mrichter/AliRoot.git/blame - PHOS/AliPHOSAnalyze.cxx
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
position of the particle initiating the hit in PHOS
[u/mrichter/AliRoot.git] / PHOS / AliPHOSAnalyze.cxx
CommitLineData
6ad0bfa0 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
b2a60966 16/* $Id$ */
17
6ad0bfa0 18//_________________________________________________________________________
5f20d3fb 19// Algorythm class to analyze PHOSv1 events:
b2a60966 20// Construct histograms and displays them.
21// Use the macro EditorBar.C for best access to the functionnalities
a3dfe79c 22//*--
2aad621e 23//*-- Author: Y. Schutz (SUBATECH) & Gines Martinez (SUBATECH)
6ad0bfa0 24//////////////////////////////////////////////////////////////////////////////
25
26// --- ROOT system ---
27
92862013 28#include "TFile.h"
29#include "TH1.h"
30#include "TPad.h"
6ad0bfa0 31#include "TH2.h"
32#include "TParticle.h"
33#include "TClonesArray.h"
34#include "TTree.h"
35#include "TMath.h"
36#include "TCanvas.h"
2bed9e3e 37#include "TStyle.h"
6ad0bfa0 38
39// --- Standard library ---
40
de9ec31b 41#include <iostream.h>
42#include <stdio.h>
92862013 43
6ad0bfa0 44// --- AliRoot header files ---
45
46#include "AliRun.h"
47#include "AliPHOSAnalyze.h"
48#include "AliPHOSClusterizerv1.h"
49#include "AliPHOSTrackSegmentMakerv1.h"
26d4b141 50#include "AliPHOSPIDv1.h"
6ad0bfa0 51#include "AliPHOSReconstructioner.h"
52#include "AliPHOSDigit.h"
53#include "AliPHOSTrackSegment.h"
54#include "AliPHOSRecParticle.h"
83974468 55#include "AliPHOSIndexToObject.h"
ed4205d8 56#include "AliPHOSHit.h"
2bed9e3e 57#include "AliPHOSCPVHit.h"
867ede0e 58#include "AliPHOSCpvRecPoint.h"
6ad0bfa0 59
60ClassImp(AliPHOSAnalyze)
61
6ad0bfa0 62//____________________________________________________________________________
63 AliPHOSAnalyze::AliPHOSAnalyze()
64{
b2a60966 65 // default ctor (useless)
6ad0bfa0 66
67 fRootFile = 0 ;
68}
69
70//____________________________________________________________________________
71AliPHOSAnalyze::AliPHOSAnalyze(Text_t * name)
72{
83974468 73 // ctor: analyze events from root file "name"
6ad0bfa0 74
92862013 75 Bool_t ok = OpenRootFile(name) ;
76 if ( !ok ) {
6ad0bfa0 77 cout << " AliPHOSAnalyze > Error opening " << name << endl ;
78 }
79 else {
eecb6765 80 //========== Get AliRun object from file
81 gAlice = (AliRun*) fRootFile->Get("gAlice") ;
fc879520 82
eecb6765 83 //=========== Get the PHOS object and associated geometry from the file
788dcca4 84 fPHOS = (AliPHOSv1 *)gAlice->GetDetector("PHOS") ;
eecb6765 85 fGeom = AliPHOSGeometry::GetInstance( fPHOS->GetGeometry()->GetName(), fPHOS->GetGeometry()->GetTitle() );
fc879520 86
eecb6765 87 //========== Initializes the Index to Object converter
88 fObjGetter = AliPHOSIndexToObject::GetInstance(fPHOS) ;
89 //========== Current event number
90 fEvt = -999 ;
fc879520 91
6ad0bfa0 92 }
2bed9e3e 93 fDebugLevel = 0;
69183710 94 fClu = 0 ;
95 fPID = 0 ;
96 fTrs = 0 ;
97 fRec = 0 ;
46b146ca 98 ResetHistograms() ;
6ad0bfa0 99}
100
88714635 101//____________________________________________________________________________
102AliPHOSAnalyze::AliPHOSAnalyze(const AliPHOSAnalyze & ana)
103{
104 // copy ctor
105 ( (AliPHOSAnalyze &)ana ).Copy(*this) ;
106}
107
108//____________________________________________________________________________
191c1c41 109void AliPHOSAnalyze::Copy(TObject & obj)
88714635 110{
111 // copy an analysis into an other one
112 TObject::Copy(obj) ;
113 // I do nothing more because the copy is silly but the Code checkers requires one
114}
115
6ad0bfa0 116//____________________________________________________________________________
117AliPHOSAnalyze::~AliPHOSAnalyze()
118{
119 // dtor
120
a3dfe79c 121 if(fRootFile->IsOpen()) fRootFile->Close() ;
122 if(fRootFile) {delete fRootFile ; fRootFile=0 ;}
123 if(fPHOS) {delete fPHOS ; fPHOS =0 ;}
124 if(fClu) {delete fClu ; fClu =0 ;}
125 if(fPID) {delete fPID ; fPID =0 ;}
126 if(fRec) {delete fRec ; fRec =0 ;}
127 if(fTrs) {delete fTrs ; fTrs =0 ;}
6ad0bfa0 128
129}
6ad0bfa0 130//____________________________________________________________________________
ed4205d8 131void AliPHOSAnalyze::DrawRecon(Int_t Nevent,Int_t Nmod){
132 //Draws pimary particles and reconstructed
133 //digits, RecPoints, RecPartices etc
134 //for event Nevent in the module Nmod.
135
136 TH2F * digitOccupancy = new TH2F("digitOccupancy","EMC digits", 64,-71.,71.,64,-71.,71.);
137 TH2F * emcOccupancy = new TH2F("emcOccupancy","EMC RecPoints",64,-71.,71.,64,-71.,71.);
138 TH2F * ppsdUp = new TH2F("ppsdUp","PPSD Up digits", 128,-71.,71.,128,-71.,71.) ;
139 TH2F * ppsdUpCl = new TH2F("ppsdUpCl","PPSD Up RecPoints",128,-71.,71.,128,-71.,71.) ;
140 TH2F * ppsdLow = new TH2F("ppsdLow","PPSD Low digits", 128,-71.,71.,128,-71.,71.) ;
141 TH2F * ppsdLowCl = new TH2F("ppsdLowCl","PPSD Low RecPoints",128,-71.,71.,128,-71.,71.) ;
142 TH2F * nbar = new TH2F("nbar","Primary nbar", 64,-71.,71.,64,-71.,71.);
143 TH2F * phot = new TH2F("phot","Primary Photon", 64,-71.,71.,64,-71.,71.);
144 TH2F * charg = new TH2F("charg","Primary charged",64,-71.,71.,64,-71.,71.);
145 TH2F * recPhot = new TH2F("recPhot","RecParticles with primary Photon",64,-71.,71.,64,-71.,71.);
146 TH2F * recNbar = new TH2F("recNbar","RecParticles with primary Nbar", 64,-71.,71.,64,-71.,71.);
147
46b146ca 148 //========== Create the Clusterizer
149 fClu = new AliPHOSClusterizerv1() ;
83974468 150
ed4205d8 151 gAlice->GetEvent(Nevent);
152
ed4205d8 153 TParticle * primary ;
154 Int_t iPrimary ;
2ab0c725 155 for ( iPrimary = 0 ; iPrimary < gAlice->GetNtrack() ; iPrimary++)
ed4205d8 156 {
2ab0c725 157 primary = gAlice->Particle(iPrimary) ;
ed4205d8 158 Int_t primaryType = primary->GetPdgCode() ;
159 if( (primaryType == 211)||(primaryType == -211)||(primaryType == 2212)||(primaryType == -2212) ) {
160 Int_t moduleNumber ;
161 Double_t primX, primZ ;
162 fGeom->ImpactOnEmc(primary->Theta(), primary->Phi(), moduleNumber, primX, primZ) ;
163 if(moduleNumber==Nmod)
164 charg->Fill(primZ,primX,primary->Energy()) ;
165 }
166 if( primaryType == 22 ) {
167 Int_t moduleNumber ;
168 Double_t primX, primZ ;
169 fGeom->ImpactOnEmc(primary->Theta(), primary->Phi(), moduleNumber, primX, primZ) ;
170 if(moduleNumber==Nmod)
171 phot->Fill(primZ,primX,primary->Energy()) ;
172 }
173 else{
174 if( primaryType == -2112 ) {
175 Int_t moduleNumber ;
176 Double_t primX, primZ ;
177 fGeom->ImpactOnEmc(primary->Theta(), primary->Phi(), moduleNumber, primX, primZ) ;
178 if(moduleNumber==Nmod)
179 nbar->Fill(primZ,primX,primary->Energy()) ;
180 }
181 }
182 }
69183710 183
ed4205d8 184 fPHOS->SetTreeAddress() ;
92862013 185
ed4205d8 186 gAlice->TreeD()->GetEvent(0) ;
187 gAlice->TreeR()->GetEvent(0) ;
188
189 TObjArray ** emcRecPoints = fPHOS->EmcRecPoints() ;
190 TObjArray ** ppsdRecPoints = fPHOS->PpsdRecPoints() ;
191 TClonesArray ** recParticleList = fPHOS->RecParticles() ;
192
193 Int_t iDigit ;
194 AliPHOSDigit * digit ;
195
196 for(iDigit = 0; iDigit < fPHOS->Digits()->GetEntries(); iDigit++)
92862013 197 {
ed4205d8 198 digit = (AliPHOSDigit *) fPHOS->Digits()->At(iDigit) ;
199 Int_t relid[4];
200 fGeom->AbsToRelNumbering(digit->GetId(), relid) ;
201 Float_t x,z ;
202 fGeom->RelPosInModule(relid,x,z) ;
203 Float_t e = fClu->Calibrate(digit->GetAmp()) ;
204 if(relid[0]==Nmod){
205 if(relid[1]==0) //EMC
206 digitOccupancy->Fill(x,z,e) ;
207 if((relid[1]>0)&&(relid[1]<17))
208 ppsdUp->Fill(x,z,e) ;
209 if(relid[1]>16)
210 ppsdLow->Fill(x,z,e) ;
211 }
212 }
213
214 Int_t irecp ;
215 TVector3 pos ;
216
217 for(irecp = 0; irecp < (*emcRecPoints)->GetEntries() ; irecp ++){
218 AliPHOSEmcRecPoint * emc= (AliPHOSEmcRecPoint*)(*emcRecPoints)->At(irecp) ;
219 if(emc->GetPHOSMod()==Nmod){
220 emc->GetLocalPosition(pos) ;
221 emcOccupancy->Fill(pos.X(),pos.Z(),emc->GetEnergy());
222 }
223 }
224
225 for(irecp = 0; irecp < (*ppsdRecPoints)->GetEntries() ; irecp ++){
226 AliPHOSPpsdRecPoint * ppsd= (AliPHOSPpsdRecPoint *)(*ppsdRecPoints)->At(irecp) ;
227 if(ppsd->GetPHOSMod()==Nmod){
228 ppsd->GetLocalPosition(pos) ;
229 if(ppsd->GetUp())
230 ppsdUpCl->Fill(pos.X(),pos.Z(),ppsd->GetEnergy());
231 else
232 ppsdLowCl->Fill(pos.X(),pos.Z(),ppsd->GetEnergy());
233 }
234 }
235
236 AliPHOSRecParticle * recParticle ;
237 Int_t iRecParticle ;
238 for(iRecParticle = 0; iRecParticle < (*recParticleList)->GetEntries() ;iRecParticle++ )
239 {
240 recParticle = (AliPHOSRecParticle *) (*recParticleList)->At(iRecParticle) ;
241
242 Int_t moduleNumberRec ;
243 Double_t recX, recZ ;
244 fGeom->ImpactOnEmc(recParticle->Theta(), recParticle->Phi(), moduleNumberRec, recX, recZ) ;
245 if(moduleNumberRec == Nmod){
246
247 Double_t minDistance = 5. ;
248 Int_t closestPrimary = -1 ;
249
250 Int_t numberofprimaries ;
251 Int_t * listofprimaries = recParticle->GetPrimaries(numberofprimaries) ;
252 Int_t index ;
253 TParticle * primary ;
254 Double_t distance = minDistance ;
11e586af 255
ed4205d8 256 for ( index = 0 ; index < numberofprimaries ; index++){
2ab0c725 257 primary = gAlice->Particle(listofprimaries[index]) ;
ed4205d8 258 Int_t moduleNumber ;
259 Double_t primX, primZ ;
260 fGeom->ImpactOnEmc(primary->Theta(), primary->Phi(), moduleNumber, primX, primZ) ;
261 if(moduleNumberRec == moduleNumber)
262 distance = TMath::Sqrt((recX-primX)*(recX-primX)+(recZ-primZ)*(recZ-primZ) ) ;
263 if(minDistance > distance)
46b146ca 264 {
ed4205d8 265 minDistance = distance ;
266 closestPrimary = listofprimaries[index] ;
46b146ca 267 }
ed4205d8 268 }
269
270 if(closestPrimary >=0 ){
271
2ab0c725 272 Int_t primaryType = gAlice->Particle(closestPrimary)->GetPdgCode() ;
ed4205d8 273
274 if(primaryType==22)
275 recPhot->Fill(recZ,recX,recParticle->Energy()) ;
276 else
277 if(primaryType==-2112)
278 recNbar->Fill(recZ,recX,recParticle->Energy()) ;
279 }
280 }
281 }
282
283
284 digitOccupancy->Draw("box") ;
285 emcOccupancy->SetLineColor(2) ;
286 emcOccupancy->Draw("boxsame") ;
287 ppsdUp->SetLineColor(3) ;
288 ppsdUp->Draw("boxsame") ;
289 ppsdLow->SetLineColor(4) ;
290 ppsdLow->Draw("boxsame") ;
291 phot->SetLineColor(8) ;
292 phot->Draw("boxsame") ;
293 nbar->SetLineColor(6) ;
294 nbar->Draw("boxsame") ;
295
296}
134ce69a 297//____________________________________________________________________________
ed4205d8 298 void AliPHOSAnalyze::Reconstruct(Int_t nevents,Int_t firstEvent )
c3b9b3f9 299{
134ce69a 300
ed4205d8 301 // Performs reconstruction of EMC and CPV (GPS2, IHEP or MIXT)
a3dfe79c 302 // for events from FirstEvent to Nevents
2bed9e3e 303
304 Int_t ievent ;
ed4205d8 305 for ( ievent=firstEvent; ievent<nevents; ievent++) {
306 if (ievent==firstEvent) {
2bed9e3e 307 cout << "Analyze > Starting Reconstructing " << endl ;
308 //========== Create the Clusterizer
309 fClu = new AliPHOSClusterizerv1() ;
fc879520 310
2bed9e3e 311 //========== Creates the track segment maker
312 fTrs = new AliPHOSTrackSegmentMakerv1() ;
fad3e5b9 313 // fTrs->UnsetUnfoldFlag() ;
2bed9e3e 314
fad3e5b9 315 //========== Creates the particle identifier
316 fPID = new AliPHOSPIDv1() ;
317 fPID->SetShowerProfileCuts(0.3, 1.8, 0.3, 1.8 ) ;
fc879520 318
2bed9e3e 319 //========== Creates the Reconstructioner
320 fRec = new AliPHOSReconstructioner(fClu, fTrs, fPID) ;
321 if (fDebugLevel != 0) fRec -> SetDebugReconstruction(kTRUE);
322 }
fad3e5b9 323
2bed9e3e 324 if (fDebugLevel != 0 ||
325 (ievent+1) % (Int_t)TMath::Power( 10, (Int_t)TMath::Log10(ievent+1) ) == 0)
326 cout << "======= Analyze ======> Event " << ievent+1 << endl ;
327
037cc66d 328 fPHOS->Enable() ;
fad3e5b9 329
037cc66d 330 gAlice->Hits2Digits() ;
fad3e5b9 331
2bed9e3e 332 //=========== Do the reconstruction
ed4205d8 333 fPHOS->Reconstruction(fRec);
fad3e5b9 334
ed4205d8 335 }
336
2bed9e3e 337 if(fClu) {delete fClu ; fClu =0 ;}
338 if(fPID) {delete fPID ; fPID =0 ;}
339 if(fRec) {delete fRec ; fRec =0 ;}
340 if(fTrs) {delete fTrs ; fTrs =0 ;}
341
342}
343
344//-------------------------------------------------------------------------------------
8c0cd6e9 345void AliPHOSAnalyze::ReadAndPrintCPV(Int_t EvFirst, Int_t EvLast)
2bed9e3e 346{
fad3e5b9 347// //
348// // Read and print generated and reconstructed hits in CPV
349// // for events from EvFirst to Nevent.
350// // If only EvFirst is defined, print only this one event.
351// // Author: Yuri Kharlov
352// // 12 October 2000
353// //
354
355// if (EvFirst!=0 && EvLast==0) EvLast=EvFirst;
356// for ( Int_t ievent=EvFirst; ievent<=EvLast; ievent++) {
8c0cd6e9 357
fad3e5b9 358// //========== Event Number>
359// cout << endl << "==== ReadAndPrintCPV ====> Event is " << ievent+1 << endl ;
2bed9e3e 360
fad3e5b9 361// //=========== Connects the various Tree's for evt
362// Int_t ntracks = gAlice->GetEvent(ievent);
2bed9e3e 363
fad3e5b9 364// //========== Creating branches ===================================
365// AliPHOSRecPoint::RecPointsList ** emcRecPoints = fPHOS->EmcRecPoints() ;
366// gAlice->TreeR()->SetBranchAddress( "PHOSEmcRP" , emcRecPoints ) ;
2bed9e3e 367
fad3e5b9 368// AliPHOSRecPoint::RecPointsList ** cpvRecPoints = fPHOS->PpsdRecPoints() ;
369// gAlice->TreeR()->SetBranchAddress( "PHOSPpsdRP", cpvRecPoints ) ;
2bed9e3e 370
fad3e5b9 371// // Read and print CPV hits
a3dfe79c 372
fad3e5b9 373// AliPHOSCPVModule cpvModule;
374// TClonesArray *cpvHits;
375// Int_t nCPVhits;
376// AliPHOSCPVHit *cpvHit;
377// TLorentzVector p;
378// Float_t xgen, zgen;
379// Int_t ipart;
380// Int_t nGenHits = 0;
381// for (Int_t itrack=0; itrack<ntracks; itrack++) {
382// //=========== Get the Hits Tree for the Primary track itrack
383// gAlice->ResetHits();
384// gAlice->TreeH()->GetEvent(itrack);
385// Int_t iModule = 0 ;
386// for (iModule=0; iModule < fGeom->GetNCPVModules(); iModule++) {
387// cpvModule = fPHOS->GetCPVModule(iModule);
388// cpvHits = cpvModule.Hits();
389// nCPVhits = cpvHits->GetEntriesFast();
390// for (Int_t ihit=0; ihit<nCPVhits; ihit++) {
391// nGenHits++;
392// cpvHit = (AliPHOSCPVHit*)cpvHits->UncheckedAt(ihit);
393// p = cpvHit->GetMomentum();
394// xgen = cpvHit->X();
395// zgen = cpvHit->Y();
396// ipart = cpvHit->GetIpart();
397// printf("CPV hit in module %d: ",iModule+1);
398// printf(" p = (%f, %f, %f, %f) GeV,\n",
399// p.Px(),p.Py(),p.Pz(),p.Energy());
400// printf(" (X,Z) = (%8.4f, %8.4f) cm, ipart = %d\n",
401// xgen,zgen,ipart);
402// }
403// }
404// }
405
406// // Read and print CPV reconstructed points
407
408// //=========== Gets the Reconstruction TTree
409// gAlice->TreeR()->GetEvent(0) ;
410// printf("Recpoints: %d\n",(*fPHOS->CpvRecPoints())->GetEntries());
411// TIter nextRP(*fPHOS->CpvRecPoints() ) ;
412// AliPHOSCpvRecPoint *cpvRecPoint ;
413// Int_t nRecPoints = 0;
414// while( ( cpvRecPoint = (AliPHOSCpvRecPoint *)nextRP() ) ) {
415// nRecPoints++;
416// TVector3 locpos;
417// cpvRecPoint->GetLocalPosition(locpos);
418// Int_t phosModule = cpvRecPoint->GetPHOSMod();
419// printf("CPV recpoint in module %d: (X,Z) = (%f,%f) cm\n",
420// phosModule,locpos.X(),locpos.Z());
421// }
422// printf("This event has %d generated hits and %d reconstructed points\n",
423// nGenHits,nRecPoints);
424// }
2bed9e3e 425}
134ce69a 426
2bed9e3e 427//____________________________________________________________________________
428void AliPHOSAnalyze::AnalyzeCPV(Int_t Nevents)
429{
037cc66d 430// //
431// // Analyzes CPV characteristics
432// // Author: Yuri Kharlov
433// // 9 October 2000
434// //
435
436// // Book histograms
437
438// TH1F *hDx = new TH1F("hDx" ,"CPV x-resolution@reconstruction",100,-5. , 5.);
439// TH1F *hDz = new TH1F("hDz" ,"CPV z-resolution@reconstruction",100,-5. , 5.);
440// TH1F *hDr = new TH1F("hDr" ,"CPV r-resolution@reconstruction",100, 0. , 5.);
441// TH1S *hNrp = new TH1S("hNrp" ,"CPV rec.point multiplicity", 21,-0.5,20.5);
442// TH1S *hNrpX = new TH1S("hNrpX","CPV rec.point Phi-length" , 21,-0.5,20.5);
443// TH1S *hNrpZ = new TH1S("hNrpZ","CPV rec.point Z-length" , 21,-0.5,20.5);
444
445// cout << "Start CPV Analysis"<< endl ;
446// for ( Int_t ievent=0; ievent<Nevents; ievent++) {
2bed9e3e 447
037cc66d 448// //========== Event Number>
449// // if ( (ievent+1) % (Int_t)TMath::Power( 10, (Int_t)TMath::Log10(ievent+1) ) == 0)
450// cout << endl << "==== AnalyzeCPV ====> Event is " << ievent+1 << endl ;
2bed9e3e 451
037cc66d 452// //=========== Connects the various Tree's for evt
453// Int_t ntracks = gAlice->GetEvent(ievent);
2bed9e3e 454
037cc66d 455// //========== Creating branches ===================================
456// AliPHOSRecPoint::RecPointsList ** emcRecPoints = fPHOS->EmcRecPoints() ;
457// gAlice->TreeR()->SetBranchAddress( "PHOSEmcRP" , emcRecPoints ) ;
2bed9e3e 458
037cc66d 459// AliPHOSRecPoint::RecPointsList ** cpvRecPoints = fPHOS->PpsdRecPoints() ;
460// gAlice->TreeR()->SetBranchAddress( "PHOSPpsdRP", cpvRecPoints ) ;
8c0cd6e9 461
037cc66d 462// // Create and fill arrays of hits for each CPV module
8c0cd6e9 463
037cc66d 464// Int_t nOfModules = fGeom->GetNModules();
465// TClonesArray **hitsPerModule = new TClonesArray *[nOfModules];
466// Int_t iModule = 0;
467// for (iModule=0; iModule < nOfModules; iModule++)
468// hitsPerModule[iModule] = new TClonesArray("AliPHOSCPVHit",100);
8c0cd6e9 469
037cc66d 470// AliPHOSCPVModule cpvModule;
471// TClonesArray *cpvHits;
472// Int_t nCPVhits;
473// AliPHOSCPVHit *cpvHit;
474// TLorentzVector p;
475// Float_t xzgen[2];
476// Int_t ipart;
134ce69a 477
037cc66d 478// // First go through all primary tracks and fill the arrays
479// // of hits per each CPV module
2bed9e3e 480
037cc66d 481// for (Int_t itrack=0; itrack<ntracks; itrack++) {
482// // Get the Hits Tree for the Primary track itrack
483// gAlice->ResetHits();
484// gAlice->TreeH()->GetEvent(itrack);
485// for (Int_t iModule=0; iModule < nOfModules; iModule++) {
486// cpvModule = fPHOS->GetCPVModule(iModule);
487// cpvHits = cpvModule.Hits();
488// nCPVhits = cpvHits->GetEntriesFast();
489// for (Int_t ihit=0; ihit<nCPVhits; ihit++) {
490// cpvHit = (AliPHOSCPVHit*)cpvHits->UncheckedAt(ihit);
491// p = cpvHit->GetMomentum();
492// xzgen[0] = cpvHit->X();
493// xzgen[1] = cpvHit->Y();
494// ipart = cpvHit->GetIpart();
495// TClonesArray &lhits = *(TClonesArray *)hitsPerModule[iModule];
496// new(lhits[hitsPerModule[iModule]->GetEntriesFast()]) AliPHOSCPVHit(*cpvHit);
497// }
498// cpvModule.Clear();
499// }
500// }
501// for (iModule=0; iModule < nOfModules; iModule++) {
502// Int_t nsum = hitsPerModule[iModule]->GetEntriesFast();
503// printf("Module %d has %d hits\n",iModule,nsum);
504// }
505
506// // Then go through reconstructed points and for each find
507// // the closeset hit
508// // The distance from the rec.point to the closest hit
509// // gives the coordinate resolution of the CPV
510
511// // Get the Reconstruction Tree
512// gAlice->TreeR()->GetEvent(0) ;
513// TIter nextRP(*fPHOS->PpsdRecPoints() ) ;
514// AliPHOSCpvRecPoint *cpvRecPoint ;
515// Float_t xgen, zgen;
516// while( ( cpvRecPoint = (AliPHOSCpvRecPoint *)nextRP() ) ) {
517// TVector3 locpos;
518// cpvRecPoint->GetLocalPosition(locpos);
519// Int_t phosModule = cpvRecPoint->GetPHOSMod();
520// Int_t rpMult = cpvRecPoint->GetDigitsMultiplicity();
521// Int_t rpMultX, rpMultZ;
522// cpvRecPoint->GetClusterLengths(rpMultX,rpMultZ);
523// Float_t xrec = locpos.X();
524// Float_t zrec = locpos.Z();
525// Float_t dxmin = 1.e+10;
526// Float_t dzmin = 1.e+10;
527// Float_t r2min = 1.e+10;
528// Float_t r2;
529
530// cpvHits = hitsPerModule[phosModule-1];
531// Int_t nCPVhits = cpvHits->GetEntriesFast();
532// for (Int_t ihit=0; ihit<nCPVhits; ihit++) {
533// cpvHit = (AliPHOSCPVHit*)cpvHits->UncheckedAt(ihit);
534// xgen = cpvHit->X();
535// zgen = cpvHit->Y();
536// r2 = TMath::Power((xgen-xrec),2) + TMath::Power((zgen-zrec),2);
537// if ( r2 < r2min ) {
538// r2min = r2;
539// dxmin = xgen - xrec;
540// dzmin = zgen - zrec;
541// }
542// }
543// hDx ->Fill(dxmin);
544// hDz ->Fill(dzmin);
545// hDr ->Fill(TMath::Sqrt(r2min));
546// hNrp ->Fill(rpMult);
547// hNrpX->Fill(rpMultX);
548// hNrpZ->Fill(rpMultZ);
549// }
550// delete [] hitsPerModule;
551// }
552// // Save histograms
553
554// Text_t outputname[80] ;
555// sprintf(outputname,"%s.analyzed",fRootFile->GetName());
556// TFile output(outputname,"RECREATE");
557// output.cd();
558
559// hDx ->Write() ;
560// hDz ->Write() ;
561// hDr ->Write() ;
562// hNrp ->Write() ;
563// hNrpX->Write() ;
564// hNrpZ->Write() ;
565
566// // Plot histograms
567
568// TCanvas *cpvCanvas = new TCanvas("CPV","CPV analysis",20,20,800,400);
569// gStyle->SetOptStat(111111);
570// gStyle->SetOptFit(1);
571// gStyle->SetOptDate(1);
572// cpvCanvas->Divide(3,2);
573
574// cpvCanvas->cd(1);
575// gPad->SetFillColor(10);
576// hNrp->SetFillColor(16);
577// hNrp->Draw();
578
579// cpvCanvas->cd(2);
580// gPad->SetFillColor(10);
581// hNrpX->SetFillColor(16);
582// hNrpX->Draw();
583
584// cpvCanvas->cd(3);
585// gPad->SetFillColor(10);
586// hNrpZ->SetFillColor(16);
587// hNrpZ->Draw();
588
589// cpvCanvas->cd(4);
590// gPad->SetFillColor(10);
591// hDx->SetFillColor(16);
592// hDx->Fit("gaus");
593// hDx->Draw();
594
595// cpvCanvas->cd(5);
596// gPad->SetFillColor(10);
597// hDz->SetFillColor(16);
598// hDz->Fit("gaus");
599// hDz->Draw();
600
601// cpvCanvas->cd(6);
602// gPad->SetFillColor(10);
603// hDr->SetFillColor(16);
604// hDr->Draw();
605
606// cpvCanvas->Print("CPV.ps");
c3b9b3f9 607
c3b9b3f9 608}
2bed9e3e 609
c3b9b3f9 610//____________________________________________________________________________
69183710 611 void AliPHOSAnalyze::InvariantMass(Int_t Nevents )
c3b9b3f9 612{
69183710 613 // Calculates Real and Mixed invariant mass distributions
ed4205d8 614
2f04ed65 615 const Int_t knMixedEvents = 4 ; //# of events used for calculation of 'mixed' distribution
616 Int_t mixedLoops = (Int_t )TMath::Ceil(Nevents/knMixedEvents) ;
69183710 617
618 //========== Booking Histograms
619 TH2D * hRealEM = new TH2D("hRealEM", "Real for EM particles", 250,0.,1.,40,0.,4.) ;
620 TH2D * hRealPhot = new TH2D("hRealPhot", "Real for kPhoton particles", 250,0.,1.,40,0.,4.) ;
621 TH2D * hMixedEM = new TH2D("hMixedEM", "Mixed for EM particles", 250,0.,1.,40,0.,4.) ;
622 TH2D * hMixedPhot= new TH2D("hMixedPhot","Mixed for kPhoton particles",250,0.,1.,40,0.,4.) ;
623
624 Int_t ievent;
ed4205d8 625 Int_t eventInMixedLoop ;
69183710 626
2f04ed65 627 Int_t nRecParticles[4];//knMixedEvents] ;
69183710 628
2f04ed65 629 AliPHOSRecParticle::RecParticlesList * allRecParticleList = new TClonesArray("AliPHOSRecParticle", knMixedEvents*1000) ;
69183710 630
ed4205d8 631 for(eventInMixedLoop = 0; eventInMixedLoop < mixedLoops; eventInMixedLoop++ ){
69183710 632 Int_t iRecPhot = 0 ;
c3b9b3f9 633
2f04ed65 634 for ( ievent=0; ievent < knMixedEvents; ievent++){
69183710 635
2f04ed65 636 Int_t absEventNumber = eventInMixedLoop*knMixedEvents + ievent ;
69183710 637
638 //=========== Connects the various Tree's for evt
ed4205d8 639 gAlice->GetEvent(absEventNumber);
640
69183710 641 //========== Creating branches ===================================
ed4205d8 642 fPHOS->SetTreeAddress() ;
69183710 643
ed4205d8 644 gAlice->TreeD()->GetEvent(0) ;
69183710 645 gAlice->TreeR()->GetEvent(0) ;
646
ed4205d8 647 TClonesArray ** recParticleList = fPHOS->RecParticles() ;
648
649
650 AliPHOSRecParticle * recParticle ;
69183710 651 Int_t iRecParticle ;
ed4205d8 652 for(iRecParticle = 0; iRecParticle < (*recParticleList)->GetEntries() ;iRecParticle++ )
69183710 653 {
ed4205d8 654 recParticle = (AliPHOSRecParticle *) (*recParticleList)->At(iRecParticle) ;
655 if((recParticle->GetType() == AliPHOSFastRecParticle::kGAMMA)||
656 (recParticle->GetType() == AliPHOSFastRecParticle::kNEUTRALEM)){
657 new( (*allRecParticleList)[iRecPhot] ) AliPHOSRecParticle(*recParticle) ;
69183710 658 iRecPhot++;
659 }
660 }
661
ed4205d8 662 nRecParticles[ievent] = iRecPhot-1 ;
c3b9b3f9 663 }
69183710 664
69183710 665 //Now calculate invariant mass:
666 Int_t irp1,irp2 ;
ed4205d8 667 Int_t nCurEvent = 0 ;
c3b9b3f9 668
ed4205d8 669 for(irp1 = 0; irp1 < allRecParticleList->GetEntries()-1; irp1++){
670 AliPHOSRecParticle * rp1 = (AliPHOSRecParticle *)allRecParticleList->At(irp1) ;
c3b9b3f9 671
ed4205d8 672 for(irp2 = irp1+1; irp2 < allRecParticleList->GetEntries(); irp2++){
673 AliPHOSRecParticle * rp2 = (AliPHOSRecParticle *)allRecParticleList->At(irp2) ;
69183710 674
ed4205d8 675 Double_t invMass ;
676 invMass = (rp1->Energy()+rp2->Energy())*(rp1->Energy()+rp2->Energy())-
69183710 677 (rp1->Px()+rp2->Px())*(rp1->Px()+rp2->Px())-
678 (rp1->Py()+rp2->Py())*(rp1->Py()+rp2->Py())-
679 (rp1->Pz()+rp2->Pz())*(rp1->Pz()+rp2->Pz()) ;
680
ed4205d8 681 if(invMass> 0)
682 invMass = TMath::Sqrt(invMass);
69183710 683
ed4205d8 684 Double_t pt ;
685 pt = TMath::Sqrt((rp1->Px()+rp2->Px() )*( rp1->Px()+rp2->Px() ) +(rp1->Py()+rp2->Py())*(rp1->Py()+rp2->Py()));
69183710 686
ed4205d8 687 if(irp1 > nRecParticles[nCurEvent])
688 nCurEvent++;
69183710 689
ed4205d8 690 if(irp2 <= nRecParticles[nCurEvent]){ //'Real' event
691 hRealEM->Fill(invMass,pt);
69183710 692 if((rp1->GetType() == AliPHOSFastRecParticle::kGAMMA)&&(rp2->GetType() == AliPHOSFastRecParticle::kGAMMA))
ed4205d8 693 hRealPhot->Fill(invMass,pt);
69183710 694 }
695 else{
ed4205d8 696 hMixedEM->Fill(invMass,pt);
69183710 697 if((rp1->GetType() == AliPHOSFastRecParticle::kGAMMA)&&(rp2->GetType() == AliPHOSFastRecParticle::kGAMMA))
ed4205d8 698 hMixedPhot->Fill(invMass,pt);
69183710 699 } //real-mixed
700
701 } //loop over second rp
702 }//loop over first rp
ed4205d8 703 allRecParticleList->Delete() ;
69183710 704 } //Loop over events
705
ed4205d8 706 delete allRecParticleList ;
69183710 707
708 //writing output
709 TFile output("invmass.root","RECREATE");
710 output.cd();
711
712 hRealEM->Write() ;
713 hRealPhot->Write() ;
714 hMixedEM->Write() ;
715 hMixedPhot->Write() ;
716
717 output.Write();
718 output.Close();
c3b9b3f9 719
720}
721
ed4205d8 722//____________________________________________________________________________
333d1c21 723 void AliPHOSAnalyze::ReadAndPrintEMC(Int_t EvFirst, Int_t EvLast)
ed4205d8 724{
037cc66d 725// //
726// // Read and print generated and reconstructed hits in EMC
727// // for events from EvFirst to Nevent.
728// // If only EvFirst is defined, print only this one event.
729// // Author: Yuri Kharlov
730// // 24 November 2000
731// //
732
733// if (EvFirst!=0 && EvLast==0) EvLast=EvFirst;
734// Int_t ievent;
735// for (ievent=EvFirst; ievent<=EvLast; ievent++) {
ed4205d8 736
037cc66d 737// //========== Event Number>
738// cout << endl << "==== ReadAndPrintEMC ====> Event is " << ievent+1 << endl ;
ed4205d8 739
037cc66d 740// //=========== Connects the various Tree's for evt
741// Int_t ntracks = gAlice->GetEvent(ievent);
742// fPHOS->SetTreeAddress() ;
ed4205d8 743
037cc66d 744// gAlice->TreeD()->GetEvent(0) ;
745// gAlice->TreeR()->GetEvent(0) ;
ed4205d8 746
037cc66d 747// // Loop over reconstructed particles
ed4205d8 748
037cc66d 749// TClonesArray ** recParticleList = fPHOS->RecParticles() ;
750// AliPHOSRecParticle * recParticle ;
751// Int_t iRecParticle ;
752// Int_t *primList;
753// Int_t nPrimary;
754// for(iRecParticle = 0; iRecParticle < (*recParticleList)->GetEntries() ;iRecParticle++ ) {
755// recParticle = (AliPHOSRecParticle *) (*recParticleList)->At(iRecParticle) ;
756// Float_t recE = recParticle->Energy();
757// primList = recParticle->GetPrimaries(nPrimary);
758// Int_t moduleNumberRec ;
759// Double_t recX, recZ ;
760// fGeom->ImpactOnEmc(recParticle->Theta(), recParticle->Phi(), moduleNumberRec, recX, recZ) ;
761// printf("Rec point: module %d, (X,Z) = (%8.4f,%8.4f) cm, E = %.3f GeV, primary = %d\n",
762// moduleNumberRec,recX,recZ,recE,*primList);
763// }
ed4205d8 764
037cc66d 765// // Read and print EMC hits from EMCn branches
ed4205d8 766
037cc66d 767// AliPHOSCPVModule emcModule;
768// TClonesArray *emcHits;
769// Int_t nEMChits;
770// AliPHOSCPVHit *emcHit;
771// TLorentzVector p;
772// Float_t xgen, zgen;
773// Int_t ipart, primary;
774// Int_t nGenHits = 0;
775// for (Int_t itrack=0; itrack<ntracks; itrack++) {
776// //=========== Get the Hits Tree for the Primary track itrack
777// gAlice->ResetHits();
778// gAlice->TreeH()->GetEvent(itrack);
779// Int_t iModule = 0 ;
780// for (iModule=0; iModule < fGeom->GetNModules(); iModule++) {
781// emcModule = fPHOS->GetEMCModule(iModule);
782// emcHits = emcModule.Hits();
783// nEMChits = emcHits->GetEntriesFast();
784// for (Int_t ihit=0; ihit<nEMChits; ihit++) {
785// nGenHits++;
786// emcHit = (AliPHOSCPVHit*)emcHits->UncheckedAt(ihit);
787// p = emcHit->GetMomentum();
788// xgen = emcHit->X();
789// zgen = emcHit->Y();
790// ipart = emcHit->GetIpart();
791// primary= emcHit->GetTrack();
792// printf("EMC hit A: module %d, ",iModule+1);
793// printf(" p = (%f .4, %f .4, %f .4, %f .4) GeV,\n",
794// p.Px(),p.Py(),p.Pz(),p.Energy());
795// printf(" (X,Z) = (%8.4f, %8.4f) cm, ipart = %d, primary = %d\n",
796// xgen,zgen,ipart,primary);
797// }
798// }
799// }
ed4205d8 800
037cc66d 801// // // Read and print EMC hits from PHOS branch
802
803// // for (Int_t itrack=0; itrack<ntracks; itrack++) {
804// // //=========== Get the Hits Tree for the Primary track itrack
805// // gAlice->ResetHits();
806// // gAlice->TreeH()->GetEvent(itrack);
807// // TClonesArray *hits = fPHOS->Hits();
808// // AliPHOSHit *hit ;
809// // Int_t ihit;
810// // for ( ihit = 0 ; ihit < hits->GetEntries() ; ihit++ ) {
811// // hit = (AliPHOSHit*)hits->At(ihit) ;
812// // Float_t hitXYZ[3];
813// // hitXYZ[0] = hit->X();
814// // hitXYZ[1] = hit->Y();
815// // hitXYZ[2] = hit->Z();
816// // ipart = hit->GetPid();
817// // primary = hit->GetPrimary();
818// // Int_t absId = hit->GetId();
819// // Int_t relId[4];
820// // fGeom->AbsToRelNumbering(absId, relId) ;
821// // Int_t module = relId[0];
822// // if (relId[1]==0 && !(hitXYZ[0]==0 && hitXYZ[2]==0))
823// // printf("EMC hit B: module %d, (X,Z) = (%8.4f, %8.4f) cm, ipart = %d, primary = %d\n",
824// // module,hitXYZ[0],hitXYZ[2],ipart,primary);
825// // }
826// // }
ed4205d8 827
037cc66d 828// }
ed4205d8 829}
830
831//____________________________________________________________________________
832 void AliPHOSAnalyze::AnalyzeEMC(Int_t Nevents)
833{
037cc66d 834// //
835// // Read generated and reconstructed hits in EMC for Nevents events.
836// // Plots the coordinate and energy resolution histograms.
837// // Coordinate resolution is a difference between the reconstructed
838// // coordinate and the exact coordinate on the face of the PHOS
839// // Author: Yuri Kharlov
840// // 27 November 2000
841// //
842
843// // Book histograms
844
845// TH1F *hDx1 = new TH1F("hDx1" ,"EMC x-resolution", 100,-5. , 5.);
846// TH1F *hDz1 = new TH1F("hDz1" ,"EMC z-resolution", 100,-5. , 5.);
847// TH1F *hDE1 = new TH1F("hDE1" ,"EMC E-resolution", 100,-2. , 2.);
848
849// TH2F *hDx2 = new TH2F("hDx2" ,"EMC x-resolution", 100, 0., 10., 100,-5. , 5.);
850// TH2F *hDz2 = new TH2F("hDz2" ,"EMC z-resolution", 100, 0., 10., 100,-5. , 5.);
851// TH2F *hDE2 = new TH2F("hDE2" ,"EMC E-resolution", 100, 0., 10., 100, 0. , 5.);
852
853// cout << "Start EMC Analysis"<< endl ;
854// for (Int_t ievent=0; ievent<Nevents; ievent++) {
ed4205d8 855
037cc66d 856// //========== Event Number>
857// if ( (ievent+1) % (Int_t)TMath::Power( 10, (Int_t)TMath::Log10(ievent+1) ) == 0)
858// cout << "==== AnalyzeEMC ====> Event is " << ievent+1 << endl ;
ed4205d8 859
037cc66d 860// //=========== Connects the various Tree's for evt
861// Int_t ntracks = gAlice->GetEvent(ievent);
ed4205d8 862
037cc66d 863// fPHOS->SetTreeAddress() ;
ed4205d8 864
037cc66d 865// gAlice->TreeD()->GetEvent(0) ;
866// gAlice->TreeR()->GetEvent(0) ;
ed4205d8 867
037cc66d 868// // Create and fill arrays of hits for each EMC module
ed4205d8 869
037cc66d 870// Int_t nOfModules = fGeom->GetNModules();
871// TClonesArray **hitsPerModule = new TClonesArray *[nOfModules];
872// Int_t iModule;
873// for (iModule=0; iModule < nOfModules; iModule++)
874// hitsPerModule[iModule] = new TClonesArray("AliPHOSCPVHit",100);
ed4205d8 875
037cc66d 876// AliPHOSCPVModule emcModule;
877// TClonesArray *emcHits;
878// Int_t nEMChits;
879// AliPHOSCPVHit *emcHit;
ed4205d8 880
037cc66d 881// // First go through all primary tracks and fill the arrays
882// // of hits per each EMC module
883
884// for (Int_t itrack=0; itrack<ntracks; itrack++) {
885// // Get the Hits Tree for the Primary track itrack
886// gAlice->ResetHits();
887// gAlice->TreeH()->GetEvent(itrack);
888// for (Int_t iModule=0; iModule < nOfModules; iModule++) {
889// emcModule = fPHOS->GetEMCModule(iModule);
890// emcHits = emcModule.Hits();
891// nEMChits = emcHits->GetEntriesFast();
892// for (Int_t ihit=0; ihit<nEMChits; ihit++) {
893// emcHit = (AliPHOSCPVHit*)emcHits->UncheckedAt(ihit);
894// TClonesArray &lhits = *(TClonesArray *)hitsPerModule[iModule];
895// new(lhits[hitsPerModule[iModule]->GetEntriesFast()]) AliPHOSCPVHit(*emcHit);
896// }
897// emcModule.Clear();
898// }
899// }
ed4205d8 900
037cc66d 901// // Loop over reconstructed particles
902
903// TClonesArray ** recParticleList = fPHOS->RecParticles() ;
904// AliPHOSRecParticle * recParticle ;
905// Int_t nEMCrecs = (*recParticleList)->GetEntries();
906// if (nEMCrecs == 1) {
907// recParticle = (AliPHOSRecParticle *) (*recParticleList)->At(0) ;
908// Float_t recE = recParticle->Energy();
909// Int_t phosModule;
910// Double_t recX, recZ ;
911// fGeom->ImpactOnEmc(recParticle->Theta(), recParticle->Phi(), phosModule, recX, recZ) ;
912
913// // for this rec.point take the hit list in the same PHOS module
914
915// emcHits = hitsPerModule[phosModule-1];
916// Int_t nEMChits = emcHits->GetEntriesFast();
917// if (nEMChits == 1) {
918// Float_t genX, genZ, genE;
919// for (Int_t ihit=0; ihit<nEMChits; ihit++) {
920// emcHit = (AliPHOSCPVHit*)emcHits->UncheckedAt(ihit);
921// genX = emcHit->X();
922// genZ = emcHit->Y();
923// genE = emcHit->GetMomentum().E();
924// }
925// Float_t dx = recX - genX;
926// Float_t dz = recZ - genZ;
927// Float_t de = recE - genE;
928// hDx1 ->Fill(dx);
929// hDz1 ->Fill(dz);
930// hDE1 ->Fill(de);
931// hDx2 ->Fill(genE,dx);
932// hDz2 ->Fill(genE,dz);
933// hDE2 ->Fill(genE,recE);
934// }
935// }
936// delete [] hitsPerModule;
937// }
938// // Save histograms
939
940// Text_t outputname[80] ;
941// sprintf(outputname,"%s.analyzed",fRootFile->GetName());
942// TFile output(outputname,"RECREATE");
943// output.cd();
944
945// hDx1 ->Write() ;
946// hDz1 ->Write() ;
947// hDE1 ->Write() ;
948// hDx2 ->Write() ;
949// hDz2 ->Write() ;
950// hDE2 ->Write() ;
951
952// // Plot histograms
953
954// TCanvas *emcCanvas = new TCanvas("EMC","EMC analysis",20,20,700,300);
955// gStyle->SetOptStat(111111);
956// gStyle->SetOptFit(1);
957// gStyle->SetOptDate(1);
958// emcCanvas->Divide(3,1);
959
960// emcCanvas->cd(1);
961// gPad->SetFillColor(10);
962// hDx1->SetFillColor(16);
963// hDx1->Draw();
964
965// emcCanvas->cd(2);
966// gPad->SetFillColor(10);
967// hDz1->SetFillColor(16);
968// hDz1->Draw();
969
970// emcCanvas->cd(3);
971// gPad->SetFillColor(10);
972// hDE1->SetFillColor(16);
973// hDE1->Draw();
974
975// emcCanvas->Print("EMC.ps");
ed4205d8 976
977}
978
fc879520 979//____________________________________________________________________________
980 void AliPHOSAnalyze::AnalyzeResolutions(Int_t Nevents )
981{
982 // analyzes Nevents events and calculate Energy and Position resolution as well as
983 // probaility of correct indentifiing of the incident particle
134ce69a 984
fc879520 985 //========== Booking Histograms
986 cout << "AnalyzeResolutions > " << "Booking Histograms" << endl ;
987 BookResolutionHistograms();
134ce69a 988
ed4205d8 989 Int_t counter[9][5] ;
990 Int_t i1,i2,totalInd = 0 ;
fc879520 991 for(i1 = 0; i1<9; i1++)
992 for(i2 = 0; i2<5; i2++)
ed4205d8 993 counter[i1][i2] = 0 ;
fc879520 994
ed4205d8 995 Int_t totalPrimary = 0 ;
996 Int_t totalRecPart = 0 ;
997 Int_t totalRPwithPrim = 0 ;
fc879520 998 Int_t ievent;
999
1000 cout << "Start Analysing"<< endl ;
1001 for ( ievent=0; ievent<Nevents; ievent++)
1002 {
1003
1004 //========== Event Number>
69183710 1005 // if ( ( log10((Float_t)(ievent+1)) - (Int_t)(log10((Float_t)(ievent+1))) ) == 0. )
fc879520 1006 cout << "AnalyzeResolutions > " << "Event is " << ievent << endl ;
1007
1008 //=========== Connects the various Tree's for evt
1009 gAlice->GetEvent(ievent);
134ce69a 1010
69183710 1011 //=========== Gets the Kine TTree
fc879520 1012 gAlice->TreeK()->GetEvent(0) ;
1013
1014 //=========== Gets the list of Primari Particles
fc879520 1015
ed4205d8 1016 TParticle * primary ;
fc879520 1017 Int_t iPrimary ;
2ab0c725 1018 for ( iPrimary = 0 ; iPrimary < gAlice->GetNtrack() ; iPrimary++)
69183710 1019 {
2ab0c725 1020 primary = gAlice->Particle(iPrimary) ;
ed4205d8 1021 Int_t primaryType = primary->GetPdgCode() ;
1022 if( primaryType == 22 ) {
1023 Int_t moduleNumber ;
1024 Double_t primX, primZ ;
1025 fGeom->ImpactOnEmc(primary->Theta(), primary->Phi(), moduleNumber, primX, primZ) ;
1026 if(moduleNumber){
1027 fhPrimary->Fill(primary->Energy()) ;
1028 if(primary->Energy() > 0.3)
1029 totalPrimary++ ;
69183710 1030 }
1031 }
1032 }
fc879520 1033
ed4205d8 1034 fPHOS->SetTreeAddress() ;
69183710 1035
ed4205d8 1036 gAlice->TreeD()->GetEvent(0) ;
fc879520 1037 gAlice->TreeR()->GetEvent(0) ;
69183710 1038
ed4205d8 1039 TClonesArray ** recParticleList = fPHOS->RecParticles() ;
1040
1041 AliPHOSRecParticle * recParticle ;
fc879520 1042 Int_t iRecParticle ;
ed4205d8 1043 for(iRecParticle = 0; iRecParticle < (*recParticleList)->GetEntries() ;iRecParticle++ )
fc879520 1044 {
ed4205d8 1045 recParticle = (AliPHOSRecParticle *) (*recParticleList)->At(iRecParticle) ;
1046 fhAllRP->Fill(CorrectEnergy(recParticle->Energy())) ;
fc879520 1047
ed4205d8 1048 Int_t moduleNumberRec ;
1049 Double_t recX, recZ ;
1050 fGeom->ImpactOnEmc(recParticle->Theta(), recParticle->Phi(), moduleNumberRec, recX, recZ) ;
fc879520 1051
ed4205d8 1052 Double_t minDistance = 100. ;
1053 Int_t closestPrimary = -1 ;
fc879520 1054
1055 Int_t numberofprimaries ;
ed4205d8 1056 Int_t * listofprimaries = recParticle->GetPrimaries(numberofprimaries) ;
fc879520 1057 Int_t index ;
ed4205d8 1058 TParticle * primary ;
1059 Double_t distance = minDistance ;
1060 Double_t dX, dZ;
b3445972 1061 Double_t dXmin = 0.;
1062 Double_t dZmin = 0. ;
69183710 1063 for ( index = 0 ; index < numberofprimaries ; index++){
2ab0c725 1064 primary = gAlice->Particle(listofprimaries[index]) ;
ed4205d8 1065 Int_t moduleNumber ;
1066 Double_t primX, primZ ;
1067 fGeom->ImpactOnEmc(primary->Theta(), primary->Phi(), moduleNumber, primX, primZ) ;
1068 if(moduleNumberRec == moduleNumber) {
1069 dX = recX - primX;
1070 dZ = recZ - primZ;
1071 distance = TMath::Sqrt(dX*dX + dZ*dZ) ;
1072 if(minDistance > distance) {
1073 minDistance = distance ;
1074 dXmin = dX;
1075 dZmin = dZ;
1076 closestPrimary = listofprimaries[index] ;
69183710 1077 }
ed4205d8 1078 }
69183710 1079 }
ed4205d8 1080 totalRecPart++ ;
134ce69a 1081
ed4205d8 1082 if(closestPrimary >=0 ){
1083 totalRPwithPrim++;
69183710 1084
2ab0c725 1085 Int_t primaryType = gAlice->Particle(closestPrimary)->GetPdgCode() ;
1086// TParticlePDG* pDGparticle = gAlice->ParticleAt(closestPrimary)->GetPDG();
2bed9e3e 1087// Double_t charge = PDGparticle->Charge() ;
8f5ada7b 1088// if(charge)
ed4205d8 1089// cout <<"Primary " <<primaryType << " E " << ((TParticle *)primaryList->At(closestPrimary))->Energy() << endl ;
1090 Int_t primaryCode ;
1091 switch(primaryType)
69183710 1092 {
1093 case 22:
ed4205d8 1094 primaryCode = 0; //Photon
2ab0c725 1095 fhAllEnergy ->Fill(gAlice->Particle(closestPrimary)->Energy(), recParticle->Energy()) ;
1096 fhAllPosition ->Fill(gAlice->Particle(closestPrimary)->Energy(), minDistance) ;
ed4205d8 1097 fhAllPositionX->Fill(dXmin);
1098 fhAllPositionZ->Fill(dZmin);
69183710 1099 break;
1100 case 11 :
ed4205d8 1101 primaryCode = 1; //Electron
69183710 1102 break;
1103 case -11 :
ed4205d8 1104 primaryCode = 1; //positron
69183710 1105 break;
1106 case 321 :
ed4205d8 1107 primaryCode = 4; //K+
69183710 1108 break;
1109 case -321 :
ed4205d8 1110 primaryCode = 4; //K-
69183710 1111 break;
1112 case 310 :
ed4205d8 1113 primaryCode = 4; //K0s
69183710 1114 break;
1115 case 130 :
ed4205d8 1116 primaryCode = 4; //K0l
69183710 1117 break;
8f5ada7b 1118 case 211 :
ed4205d8 1119 primaryCode = 2; //K0l
8f5ada7b 1120 break;
1121 case -211 :
ed4205d8 1122 primaryCode = 2; //K0l
8f5ada7b 1123 break;
1124 case 2212 :
ed4205d8 1125 primaryCode = 2; //K0l
8f5ada7b 1126 break;
1127 case -2212 :
ed4205d8 1128 primaryCode = 2; //K0l
8f5ada7b 1129 break;
69183710 1130 default:
ed4205d8 1131 primaryCode = 3; //ELSE
69183710 1132 break;
1133 }
1134
ed4205d8 1135 switch(recParticle->GetType())
69183710 1136 {
1137 case AliPHOSFastRecParticle::kGAMMA:
ed4205d8 1138 if(primaryType == 22){
2ab0c725 1139 fhPhotEnergy->Fill(gAlice->Particle(closestPrimary)->Energy(), recParticle->Energy() ) ;
1140 fhEMEnergy->Fill(gAlice->Particle(closestPrimary)->Energy(), recParticle->Energy() ) ;
1141 fhPPSDEnergy->Fill(gAlice->Particle(closestPrimary)->Energy(), recParticle->Energy() ) ;
69183710 1142
2ab0c725 1143 fhPhotPosition->Fill(gAlice->Particle(closestPrimary)->Energy(),minDistance) ;
1144 fhEMPosition->Fill(gAlice->Particle(closestPrimary)->Energy(),minDistance) ;
1145 fhPPSDPosition->Fill(gAlice->Particle(closestPrimary)->Energy(),minDistance) ;
69183710 1146
ed4205d8 1147 fhPhotReg->Fill(CorrectEnergy(recParticle->Energy()) ) ;
1148 fhPhotEM->Fill(CorrectEnergy(recParticle->Energy()) ) ;
1149 fhPhotPPSD->Fill(CorrectEnergy(recParticle->Energy()) ) ;
69183710 1150
ed4205d8 1151 fhPhotPhot->Fill(CorrectEnergy(recParticle->Energy()) ) ;
69183710 1152 }
ed4205d8 1153 if(primaryType == 2112){ //neutron
1154 fhNReg->Fill(CorrectEnergy(recParticle->Energy()) ) ;
1155 fhNEM->Fill(CorrectEnergy(recParticle->Energy()) ) ;
1156 fhNPPSD->Fill(CorrectEnergy(recParticle->Energy()) ) ;
69183710 1157 }
1158
ed4205d8 1159 if(primaryType == -2112){ //neutron ~
1160 fhNBarReg->Fill(CorrectEnergy(recParticle->Energy()) ) ;
1161 fhNBarEM->Fill(CorrectEnergy(recParticle->Energy()) ) ;
1162 fhNBarPPSD->Fill(CorrectEnergy(recParticle->Energy()) ) ;
69183710 1163
1164 }
ed4205d8 1165 if(primaryCode == 2){
1166 fhChargedReg->Fill(CorrectEnergy(recParticle->Energy()) ) ;
1167 fhChargedEM->Fill(CorrectEnergy(recParticle->Energy()) ) ;
1168 fhChargedPPSD->Fill(CorrectEnergy(recParticle->Energy()) ) ;
69183710 1169 }
1170
ed4205d8 1171 fhAllReg->Fill(CorrectEnergy(recParticle->Energy()) ) ;
1172 fhAllEM->Fill(CorrectEnergy(recParticle->Energy()) ) ;
1173 fhAllPPSD->Fill(CorrectEnergy(recParticle->Energy()) ) ;
1174 fhShape->Fill(CorrectEnergy(recParticle->Energy()) ) ;
1175 fhVeto->Fill(CorrectEnergy(recParticle->Energy()) ) ;
1176 fhPPSD->Fill(CorrectEnergy(recParticle->Energy()) ) ;
1177 counter[0][primaryCode]++;
69183710 1178 break;
1179 case AliPHOSFastRecParticle::kELECTRON:
ed4205d8 1180 if(primaryType == 22){
1181 fhPhotElec->Fill(CorrectEnergy(recParticle->Energy()) ) ;
2ab0c725 1182 fhEMEnergy->Fill(gAlice->Particle(closestPrimary)->Energy(), recParticle->Energy() ) ;
1183 fhEMPosition->Fill(gAlice->Particle(closestPrimary)->Energy(),minDistance) ;
ed4205d8 1184 fhPhotEM->Fill(CorrectEnergy(recParticle->Energy()) ) ;
1185 fhPhotPPSD->Fill(CorrectEnergy(recParticle->Energy()) ) ;
69183710 1186 }
ed4205d8 1187 if(primaryType == 2112){ //neutron
1188 fhNEM->Fill(CorrectEnergy(recParticle->Energy()) ) ;
1189 fhNPPSD->Fill(CorrectEnergy(recParticle->Energy()) ) ;
69183710 1190 }
1191
ed4205d8 1192 if(primaryType == -2112){ //neutron ~
1193 fhNBarEM->Fill(CorrectEnergy(recParticle->Energy()) ) ;
1194 fhNBarPPSD->Fill(CorrectEnergy(recParticle->Energy()) ) ;
69183710 1195
1196 }
ed4205d8 1197 if(primaryCode == 2){
1198 fhChargedEM->Fill(CorrectEnergy(recParticle->Energy()) ) ;
1199 fhChargedPPSD->Fill(CorrectEnergy(recParticle->Energy()) ) ;
69183710 1200 }
1201
ed4205d8 1202 fhAllEM->Fill(CorrectEnergy(recParticle->Energy()) ) ;
1203 fhAllPPSD->Fill(CorrectEnergy(recParticle->Energy()) ) ;
1204 fhShape->Fill(CorrectEnergy(recParticle->Energy()) ) ;
1205 fhPPSD->Fill(CorrectEnergy(recParticle->Energy()) ) ;
1206 counter[1][primaryCode]++;
69183710 1207 break;
1208 case AliPHOSFastRecParticle::kNEUTRALHA:
ed4205d8 1209 if(primaryType == 22)
1210 fhPhotNeuH->Fill(CorrectEnergy(recParticle->Energy()) ) ;
69183710 1211
ed4205d8 1212 fhVeto->Fill(CorrectEnergy(recParticle->Energy()) ) ;
1213 counter[2][primaryCode]++;
69183710 1214 break ;
1215 case AliPHOSFastRecParticle::kNEUTRALEM:
ed4205d8 1216 if(primaryType == 22){
2ab0c725 1217 fhEMEnergy->Fill(gAlice->Particle(closestPrimary)->Energy(),recParticle->Energy() ) ;
1218 fhEMPosition->Fill(gAlice->Particle(closestPrimary)->Energy(),minDistance ) ;
69183710 1219
ed4205d8 1220 fhPhotNuEM->Fill(CorrectEnergy(recParticle->Energy()) ) ;
1221 fhPhotEM->Fill(CorrectEnergy(recParticle->Energy()) ) ;
69183710 1222 }
ed4205d8 1223 if(primaryType == 2112) //neutron
1224 fhNEM->Fill(CorrectEnergy(recParticle->Energy()) ) ;
69183710 1225
ed4205d8 1226 if(primaryType == -2112) //neutron ~
1227 fhNBarEM->Fill(CorrectEnergy(recParticle->Energy()) ) ;
69183710 1228
ed4205d8 1229 if(primaryCode == 2)
1230 fhChargedEM->Fill(CorrectEnergy(recParticle->Energy()) ) ;
69183710 1231
ed4205d8 1232 fhAllEM->Fill(CorrectEnergy(recParticle->Energy()) ) ;
1233 fhShape->Fill(CorrectEnergy(recParticle->Energy()) ) ;
1234 fhVeto->Fill(CorrectEnergy(recParticle->Energy()) ) ;
69183710 1235
ed4205d8 1236 counter[3][primaryCode]++;
69183710 1237 break ;
1238 case AliPHOSFastRecParticle::kCHARGEDHA:
ed4205d8 1239 if(primaryType == 22) //photon
1240 fhPhotChHa->Fill(CorrectEnergy(recParticle->Energy()) ) ;
69183710 1241
ed4205d8 1242 counter[4][primaryCode]++ ;
69183710 1243 break ;
1244 case AliPHOSFastRecParticle::kGAMMAHA:
ed4205d8 1245 if(primaryType == 22){ //photon
1246 fhPhotGaHa->Fill(CorrectEnergy(recParticle->Energy()) ) ;
2ab0c725 1247 fhPPSDEnergy->Fill(gAlice->Particle(closestPrimary)->Energy(), recParticle->Energy() ) ;
1248 fhPPSDPosition->Fill(gAlice->Particle(closestPrimary)->Energy(),minDistance) ;
ed4205d8 1249 fhPhotPPSD->Fill(CorrectEnergy(recParticle->Energy()) ) ;
fc879520 1250 }
ed4205d8 1251 if(primaryType == 2112){ //neutron
1252 fhNPPSD->Fill(CorrectEnergy(recParticle->Energy()) ) ;
fc879520 1253 }
69183710 1254
ed4205d8 1255 if(primaryType == -2112){ //neutron ~
1256 fhNBarPPSD->Fill(CorrectEnergy(recParticle->Energy()) ) ;
fc879520 1257 }
ed4205d8 1258 if(primaryCode == 2){
1259 fhChargedPPSD->Fill(CorrectEnergy(recParticle->Energy()) ) ;
fc879520 1260 }
69183710 1261
ed4205d8 1262 fhAllPPSD->Fill(CorrectEnergy(recParticle->Energy()) ) ;
1263 fhVeto->Fill(CorrectEnergy(recParticle->Energy()) ) ;
1264 fhPPSD->Fill(CorrectEnergy(recParticle->Energy()) ) ;
1265 counter[5][primaryCode]++ ;
fc879520 1266 break ;
69183710 1267 case AliPHOSFastRecParticle::kABSURDEM:
ed4205d8 1268 counter[6][primaryCode]++ ;
1269 fhShape->Fill(CorrectEnergy(recParticle->Energy()) ) ;
69183710 1270 break;
1271 case AliPHOSFastRecParticle::kABSURDHA:
ed4205d8 1272 counter[7][primaryCode]++ ;
69183710 1273 break;
1274 default:
ed4205d8 1275 counter[8][primaryCode]++ ;
69183710 1276 break;
1277 }
1278 }
fc879520 1279 }
1280 } // endfor
46b146ca 1281 SaveHistograms();
fc879520 1282 cout << "Resolutions: Analyzed " << Nevents << " event(s)" << endl ;
ed4205d8 1283 cout << "Resolutions: Total primary " << totalPrimary << endl ;
1284 cout << "Resoluitons: Total reconstracted " << totalRecPart << endl ;
1285 cout << "TotalReconstructed with Primarie " << totalRPwithPrim << endl ;
fc879520 1286 cout << " Primary: Photon Electron Ch. Hadr. Neutr. Hadr Kaons" << endl ;
ed4205d8 1287 cout << " Detected as photon " << counter[0][0] << " " << counter[0][1] << " " << counter[0][2] << " " <<counter[0][3] << " " << counter[0][4] << endl ;
1288 cout << " Detected as electron " << counter[1][0] << " " << counter[1][1] << " " << counter[1][2] << " " <<counter[1][3] << " " << counter[1][4] << endl ;
1289 cout << " Detected as neutral hadron " << counter[2][0] << " " << counter[2][1] << " " << counter[2][2] << " " <<counter[2][3] << " " << counter[2][4] << endl ;
1290 cout << " Detected as neutral EM " << counter[3][0] << " " << counter[3][1] << " " << counter[3][2] << " " <<counter[3][3] << " " << counter[3][4] << endl ;
1291 cout << " Detected as charged hadron " << counter[4][0] << " " << counter[4][1] << " " << counter[4][2] << " " <<counter[4][3] << " " << counter[4][4] << endl ;
1292 cout << " Detected as gamma-hadron " << counter[5][0] << " " << counter[5][1] << " " << counter[5][2] << " " <<counter[5][3] << " " << counter[5][4] << endl ;
1293 cout << " Detected as Absurd EM " << counter[6][0] << " " << counter[6][1] << " " << counter[6][2] << " " <<counter[6][3] << " " << counter[6][4] << endl ;
1294 cout << " Detected as absurd hadron " << counter[7][0] << " " << counter[7][1] << " " << counter[7][2] << " " <<counter[7][3] << " " << counter[7][4] << endl ;
1295 cout << " Detected as undefined " << counter[8][0] << " " << counter[8][1] << " " << counter[8][2] << " " <<counter[8][3] << " " << counter[8][4] << endl ;
fc879520 1296
1297 for(i1 = 0; i1<9; i1++)
1298 for(i2 = 0; i2<5; i2++)
ed4205d8 1299 totalInd+=counter[i1][i2] ;
1300 cout << "Indentified particles " << totalInd << endl ;
fc879520 1301
92862013 1302} // endfunction
1303
1304
1305//____________________________________________________________________________
1306void AliPHOSAnalyze::BookingHistograms()
1307{
b2a60966 1308 // Books the histograms where the results of the analysis are stored (to be changed)
1309
eecb6765 1310 delete fhEmcDigit ;
1311 delete fhVetoDigit ;
1312 delete fhConvertorDigit ;
1313 delete fhEmcCluster ;
1314 delete fhVetoCluster ;
1315 delete fhConvertorCluster ;
1316 delete fhConvertorEmc ;
1317
46b146ca 1318 fhEmcDigit = new TH1F("hEmcDigit", "hEmcDigit", 1000, 0. , 25.);
1319 fhVetoDigit = new TH1F("hVetoDigit", "hVetoDigit", 500, 0. , 3.e-5);
1320 fhConvertorDigit = new TH1F("hConvertorDigit","hConvertorDigit", 500, 0. , 3.e-5);
1321 fhEmcCluster = new TH1F("hEmcCluster", "hEmcCluster", 1000, 0. , 30.);
1322 fhVetoCluster = new TH1F("hVetoCluster", "hVetoCluster", 500, 0. , 3.e-5);
1323 fhConvertorCluster = new TH1F("hConvertorCluster","hConvertorCluster",500, 0. , 3.e-5);
1324 fhConvertorEmc = new TH2F("hConvertorEmc", "hConvertorEmc", 200, 1. , 3., 200, 0., 3.e-5);
92862013 1325
134ce69a 1326}
1327//____________________________________________________________________________
1328void AliPHOSAnalyze::BookResolutionHistograms()
1329{
1330 // Books the histograms where the results of the Resolution analysis are stored
1331
69183710 1332// if(fhAllEnergy)
1333// delete fhAllEnergy ;
1334// if(fhPhotEnergy)
1335// delete fhPhotEnergy ;
1336// if(fhEMEnergy)
1337// delete fhEMEnergy ;
1338// if(fhPPSDEnergy)
1339// delete fhPPSDEnergy ;
1340
1341
1342 fhAllEnergy = new TH2F("hAllEnergy", "Energy of any RP with primary photon",100, 0., 5., 100, 0., 5.);
1343 fhPhotEnergy = new TH2F("hPhotEnergy", "Energy of kGAMMA with primary photon",100, 0., 5., 100, 0., 5.);
1344 fhEMEnergy = new TH2F("hEMEnergy", "Energy of EM with primary photon", 100, 0., 5., 100, 0., 5.);
1345 fhPPSDEnergy = new TH2F("hPPSDEnergy", "Energy of PPSD with primary photon", 100, 0., 5., 100, 0., 5.);
1346
1347// if(fhAllPosition)
1348// delete fhAllPosition ;
1349// if(fhPhotPosition)
1350// delete fhPhotPosition ;
1351// if(fhEMPosition)
1352// delete fhEMPosition ;
1353// if(fhPPSDPosition)
1354// delete fhPPSDPosition ;
1355
1356
1357 fhAllPosition = new TH2F("hAllPosition", "Position of any RP with primary photon",100, 0., 5., 100, 0., 5.);
1358 fhPhotPosition = new TH2F("hPhotPosition", "Position of kGAMMA with primary photon",100, 0., 5., 100, 0., 5.);
1359 fhEMPosition = new TH2F("hEMPosition", "Position of EM with primary photon", 100, 0., 5., 100, 0., 5.);
1360 fhPPSDPosition = new TH2F("hPPSDPosition", "Position of PPSD with primary photon", 100, 0., 5., 100, 0., 5.);
1361
ed4205d8 1362 fhAllPositionX = new TH1F("hAllPositionX", "#Delta X of any RP with primary photon",100, -2., 2.);
1363 fhAllPositionZ = new TH1F("hAllPositionZ", "#Delta X of any RP with primary photon",100, -2., 2.);
1364
69183710 1365// if(fhAllReg)
1366// delete fhAllReg ;
1367// if(fhPhotReg)
1368// delete fhPhotReg ;
1369// if(fhNReg)
1370// delete fhNReg ;
1371// if(fhNBarReg)
1372// delete fhNBarReg ;
1373// if(fhChargedReg)
1374// delete fhChargedReg ;
46b146ca 1375
69183710 1376 fhAllReg = new TH1F("hAllReg", "All primaries registered as photon", 100, 0., 5.);
1377 fhPhotReg = new TH1F("hPhotReg", "Photon registered as photon", 100, 0., 5.);
1378 fhNReg = new TH1F("hNReg", "N registered as photon", 100, 0., 5.);
1379 fhNBarReg = new TH1F("hNBarReg", "NBar registered as photon", 100, 0., 5.);
1380 fhChargedReg= new TH1F("hChargedReg", "Charged hadron registered as photon",100, 0., 5.);
46b146ca 1381
69183710 1382// if(fhAllEM)
1383// delete fhAllEM ;
1384// if(fhPhotEM)
1385// delete fhPhotEM ;
1386// if(fhNEM)
1387// delete fhNEM ;
1388// if(fhNBarEM)
1389// delete fhNBarEM ;
1390// if(fhChargedEM)
1391// delete fhChargedEM ;
46b146ca 1392
69183710 1393 fhAllEM = new TH1F("hAllEM", "All primary registered as EM",100, 0., 5.);
1394 fhPhotEM = new TH1F("hPhotEM", "Photon registered as EM", 100, 0., 5.);
1395 fhNEM = new TH1F("hNEM", "N registered as EM", 100, 0., 5.);
1396 fhNBarEM = new TH1F("hNBarEM", "NBar registered as EM", 100, 0., 5.);
1397 fhChargedEM= new TH1F("hChargedEM","Charged registered as EM",100, 0., 5.);
1398
1399// if(fhAllPPSD)
1400// delete fhAllPPSD ;
1401// if(fhPhotPPSD)
1402// delete fhPhotPPSD ;
1403// if(fhNPPSD)
1404// delete fhNPPSD ;
1405// if(fhNBarPPSD)
1406// delete fhNBarPPSD ;
1407// if(fhChargedPPSD)
1408// delete fhChargedPPSD ;
46b146ca 1409
69183710 1410 fhAllPPSD = new TH1F("hAllPPSD", "All primary registered as PPSD",100, 0., 5.);
1411 fhPhotPPSD = new TH1F("hPhotPPSD", "Photon registered as PPSD", 100, 0., 5.);
1412 fhNPPSD = new TH1F("hNPPSD", "N registered as PPSD", 100, 0., 5.);
1413 fhNBarPPSD = new TH1F("hNBarPPSD", "NBar registered as PPSD", 100, 0., 5.);
1414 fhChargedPPSD= new TH1F("hChargedPPSD","Charged registered as PPSD",100, 0., 5.);
1415
1416// if(fhPrimary)
1417// delete fhPrimary ;
1418 fhPrimary= new TH1F("hPrimary", "hPrimary", 100, 0., 5.);
1419
1420// if(fhAllRP)
1421// delete fhAllRP ;
1422// if(fhVeto)
1423// delete fhVeto ;
1424// if(fhShape)
1425// delete fhShape ;
1426// if(fhPPSD)
1427// delete fhPPSD ;
1428
1429 fhAllRP = new TH1F("hAllRP","All Reconstructed particles", 100, 0., 5.);
1430 fhVeto = new TH1F("hVeto", "All uncharged particles", 100, 0., 5.);
1431 fhShape = new TH1F("hShape","All particles with EM shaower",100, 0., 5.);
1432 fhPPSD = new TH1F("hPPSD", "All PPSD photon particles", 100, 0., 5.);
1433
1434
1435// if(fhPhotPhot)
1436// delete fhPhotPhot ;
1437// if(fhPhotElec)
1438// delete fhPhotElec ;
1439// if(fhPhotNeuH)
1440// delete fhPhotNeuH ;
1441// if(fhPhotNuEM)
1442// delete fhPhotNuEM ;
1443// if(fhPhotChHa)
1444// delete fhPhotChHa ;
1445// if(fhPhotGaHa)
1446// delete fhPhotGaHa ;
1447
1448 fhPhotPhot = new TH1F("hPhotPhot","hPhotPhot", 100, 0., 5.); //Photon registered as photon
1449 fhPhotElec = new TH1F("hPhotElec","hPhotElec", 100, 0., 5.); //Photon registered as Electron
1450 fhPhotNeuH = new TH1F("hPhotNeuH","hPhotNeuH", 100, 0., 5.); //Photon registered as Neutral Hadron
1451 fhPhotNuEM = new TH1F("hPhotNuEM","hPhotNuEM", 100, 0., 5.); //Photon registered as Neutral EM
1452 fhPhotChHa = new TH1F("hPhotChHa","hPhotChHa", 100, 0., 5.); //Photon registered as Charged Hadron
1453 fhPhotGaHa = new TH1F("hPhotGaHa","hPhotGaHa", 100, 0., 5.); //Photon registered as Gamma-Hadron
92862013 1454}
2aad621e 1455
6ad0bfa0 1456//____________________________________________________________________________
1457Bool_t AliPHOSAnalyze::OpenRootFile(Text_t * name)
1458{
b2a60966 1459 // Open the root file named "name"
1460
1461 fRootFile = new TFile(name, "update") ;
6ad0bfa0 1462 return fRootFile->IsOpen() ;
1463}
ed4205d8 1464
92862013 1465//____________________________________________________________________________
46b146ca 1466void AliPHOSAnalyze::SaveHistograms()
134ce69a 1467{
1468 // Saves the histograms in a root file named "name.analyzed"
1469
1470 Text_t outputname[80] ;
1471 sprintf(outputname,"%s.analyzed",fRootFile->GetName());
1472 TFile output(outputname,"RECREATE");
1473 output.cd();
46b146ca 1474
69183710 1475 if (fhAllEnergy)
1476 fhAllEnergy->Write() ;
1477 if (fhPhotEnergy)
1478 fhPhotEnergy->Write() ;
1479 if(fhEMEnergy)
1480 fhEMEnergy->Write() ;
1481 if(fhPPSDEnergy)
1482 fhPPSDEnergy->Write() ;
1483 if(fhAllPosition)
1484 fhAllPosition->Write() ;
ed4205d8 1485 if(fhAllPositionX)
1486 fhAllPositionX->Write() ;
1487 if(fhAllPositionZ)
1488 fhAllPositionZ->Write() ;
69183710 1489 if(fhPhotPosition)
1490 fhPhotPosition->Write() ;
1491 if(fhEMPosition)
1492 fhEMPosition->Write() ;
1493 if(fhPPSDPosition)
1494 fhPPSDPosition->Write() ;
fc879520 1495 if (fhAllReg)
1496 fhAllReg->Write() ;
69183710 1497 if (fhPhotReg)
1498 fhPhotReg->Write() ;
fc879520 1499 if(fhNReg)
1500 fhNReg->Write() ;
1501 if(fhNBarReg)
1502 fhNBarReg->Write() ;
1503 if(fhChargedReg)
1504 fhChargedReg->Write() ;
fc879520 1505 if (fhAllEM)
1506 fhAllEM->Write() ;
69183710 1507 if (fhPhotEM)
1508 fhPhotEM->Write() ;
fc879520 1509 if(fhNEM)
1510 fhNEM->Write() ;
1511 if(fhNBarEM)
1512 fhNBarEM->Write() ;
1513 if(fhChargedEM)
1514 fhChargedEM->Write() ;
69183710 1515 if (fhAllPPSD)
1516 fhAllPPSD->Write() ;
1517 if (fhPhotPPSD)
1518 fhPhotPPSD->Write() ;
1519 if(fhNPPSD)
1520 fhNPPSD->Write() ;
1521 if(fhNBarPPSD)
1522 fhNBarPPSD->Write() ;
1523 if(fhChargedPPSD)
1524 fhChargedPPSD->Write() ;
fc879520 1525 if(fhPrimary)
1526 fhPrimary->Write() ;
69183710 1527 if(fhAllRP)
1528 fhAllRP->Write() ;
1529 if(fhVeto)
1530 fhVeto->Write() ;
1531 if(fhShape)
1532 fhShape->Write() ;
1533 if(fhPPSD)
1534 fhPPSD->Write() ;
fc879520 1535 if(fhPhotPhot)
1536 fhPhotPhot->Write() ;
1537 if(fhPhotElec)
1538 fhPhotElec->Write() ;
1539 if(fhPhotNeuH)
1540 fhPhotNeuH->Write() ;
1541 if(fhPhotNuEM)
1542 fhPhotNuEM->Write() ;
1543 if(fhPhotNuEM)
1544 fhPhotNuEM->Write() ;
1545 if(fhPhotChHa)
1546 fhPhotChHa->Write() ;
1547 if(fhPhotGaHa)
1548 fhPhotGaHa->Write() ;
46b146ca 1549 if(fhEnergyCorrelations)
1550 fhEnergyCorrelations->Write() ;
1551
92862013 1552 output.Write();
1553 output.Close();
1554}
69183710 1555//____________________________________________________________________________
1556Float_t AliPHOSAnalyze::CorrectEnergy(Float_t ERecPart)
1557{
1558 return ERecPart/0.8783 ;
1559}
1560
46b146ca 1561//____________________________________________________________________________
1562void AliPHOSAnalyze::ResetHistograms()
1563{
1564 fhEnergyCorrelations = 0 ; //Energy correlations between Eloss in Convertor and PPSD(2)
1565
1566 fhEmcDigit = 0 ; // Histo of digit energies in the Emc
1567 fhVetoDigit = 0 ; // Histo of digit energies in the Veto
1568 fhConvertorDigit = 0 ; // Histo of digit energies in the Convertor
1569 fhEmcCluster = 0 ; // Histo of Cluster energies in Emc
1570 fhVetoCluster = 0 ; // Histo of Cluster energies in Veto
1571 fhConvertorCluster = 0 ; // Histo of Cluster energies in Convertor
1572 fhConvertorEmc = 0 ; // 2d Convertor versus Emc energies
1573
69183710 1574 fhAllEnergy = 0 ;
1575 fhPhotEnergy = 0 ; // Total spectrum of detected photons
1576 fhEMEnergy = 0 ; // Spectrum of detected electrons with electron primary
1577 fhPPSDEnergy = 0 ;
1578 fhAllPosition = 0 ;
ed4205d8 1579 fhAllPositionX = 0 ;
1580 fhAllPositionZ = 0 ;
69183710 1581 fhPhotPosition = 0 ;
1582 fhEMPosition = 0 ;
1583 fhPPSDPosition = 0 ;
1584
1585 fhPhotReg = 0 ;
46b146ca 1586 fhAllReg = 0 ;
1587 fhNReg = 0 ;
1588 fhNBarReg = 0 ;
1589 fhChargedReg = 0 ;
69183710 1590 fhPhotEM = 0 ;
46b146ca 1591 fhAllEM = 0 ;
1592 fhNEM = 0 ;
1593 fhNBarEM = 0 ;
1594 fhChargedEM = 0 ;
69183710 1595 fhPhotPPSD = 0 ;
1596 fhAllPPSD = 0 ;
1597 fhNPPSD = 0 ;
1598 fhNBarPPSD = 0 ;
1599 fhChargedPPSD = 0 ;
1600
46b146ca 1601 fhPrimary = 0 ;
1602
1603 fhPhotPhot = 0 ;
1604 fhPhotElec = 0 ;
1605 fhPhotNeuH = 0 ;
1606 fhPhotNuEM = 0 ;
1607 fhPhotChHa = 0 ;
1608 fhPhotGaHa = 0 ;
1609
46b146ca 1610}