Full PID delegated to AliPHOSPID
[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
16//_________________________________________________________________________
17// Algorythm class to analyze PHOS events
18//*-- Y. Schutz : SUBATECH
19//////////////////////////////////////////////////////////////////////////////
20
21// --- ROOT system ---
22
92862013 23#include "TFile.h"
24#include "TH1.h"
25#include "TPad.h"
6ad0bfa0 26#include "TH2.h"
27#include "TParticle.h"
28#include "TClonesArray.h"
29#include "TTree.h"
30#include "TMath.h"
31#include "TCanvas.h"
32
33// --- Standard library ---
34
92862013 35#include <iostream>
36#include <cstdio>
37
6ad0bfa0 38// --- AliRoot header files ---
39
40#include "AliRun.h"
41#include "AliPHOSAnalyze.h"
42#include "AliPHOSClusterizerv1.h"
43#include "AliPHOSTrackSegmentMakerv1.h"
26d4b141 44#include "AliPHOSPIDv1.h"
6ad0bfa0 45#include "AliPHOSReconstructioner.h"
46#include "AliPHOSDigit.h"
47#include "AliPHOSTrackSegment.h"
48#include "AliPHOSRecParticle.h"
49
50ClassImp(AliPHOSAnalyze)
51
52
53//____________________________________________________________________________
54 AliPHOSAnalyze::AliPHOSAnalyze()
55{
56 // ctor
57
58 fRootFile = 0 ;
59}
60
61//____________________________________________________________________________
62AliPHOSAnalyze::AliPHOSAnalyze(Text_t * name)
63{
64 // ctor
65
92862013 66 Bool_t ok = OpenRootFile(name) ;
67 if ( !ok ) {
6ad0bfa0 68 cout << " AliPHOSAnalyze > Error opening " << name << endl ;
69 }
70 else {
71 gAlice = (AliRun*) fRootFile->Get("gAlice");
72 fPHOS = (AliPHOSv0 *)gAlice->GetDetector("PHOS") ;
73 fGeom = AliPHOSGeometry::GetInstance( fPHOS->GetGeometry()->GetName(), fPHOS->GetGeometry()->GetTitle() ) ;
74 fEvt = -999 ;
75 }
76}
77
78//____________________________________________________________________________
79AliPHOSAnalyze::~AliPHOSAnalyze()
80{
81 // dtor
82
83 fRootFile->Close() ;
84 delete fRootFile ;
85 fRootFile = 0 ;
86
87 delete fPHOS ;
88 fPHOS = 0 ;
89
90 delete fClu ;
91 fClu = 0 ;
92
26d4b141 93 delete fPID ;
94 fPID = 0 ;
6ad0bfa0 95
96 delete fRec ;
97 fRec = 0 ;
98
99 delete fTrs ;
100 fTrs = 0 ;
101
102}
103
104//____________________________________________________________________________
105void AliPHOSAnalyze::AnalyzeOneEvent(Int_t evt)
106{
92862013 107 Bool_t ok = Init(evt) ;
6ad0bfa0 108
92862013 109 if ( ok ) {
6ad0bfa0 110 //=========== Get the number of entries in the Digits array
111
112 Int_t nId = fPHOS->Digits()->GetEntries();
113 printf("AnalyzeOneEvent > Number of entries in the Digit array is %d \n",nId);
114
115 //=========== Do the reconstruction
116
117 cout << "AnalyzeOneEvent > Found " << nId << " digits in PHOS" << endl ;
118
119 fPHOS->Reconstruction(fRec);
120
121 // =========== End of reconstruction
122
123 cout << "AnalyzeOneEvent > event # " << fEvt << " processed" << endl ;
92862013 124 } // ok
6ad0bfa0 125 else
126 cout << "AnalyzeOneEvent > filed to process event # " << evt << endl ;
127
128}
129
92862013 130//____________________________________________________________________________
131 void AliPHOSAnalyze::AnalyzeManyEvents(Int_t Nevents, Int_t module) // analyzes many events
132{
133
134 if ( fRootFile == 0 )
135 cout << "AnalyzeManyEvents > " << "Root File not openned" << endl ;
136 else
137 {
138 //========== Get AliRun object from file
139 gAlice = (AliRun*) fRootFile->Get("gAlice") ;
140 //=========== Get the PHOS object and associated geometry from the file
141 fPHOS = (AliPHOSv0 *)gAlice->GetDetector("PHOS") ;
142 fGeom = AliPHOSGeometry::GetInstance( fPHOS->GetGeometry()->GetName(), fPHOS->GetGeometry()->GetTitle() );
143 //========== Booking Histograms
144 cout << "AnalyzeManyEvents > " << "Booking Histograms" << endl ;
145 BookingHistograms();
146 Int_t ievent;
147 Int_t relid[4] ;
148 AliPHOSDigit * digit ;
149 AliPHOSEmcRecPoint * emc ;
150 AliPHOSPpsdRecPoint * ppsd ;
38f8ae6d 151 // AliPHOSTrackSegment * tracksegment ;
152 AliPHOSRecParticle * recparticle;
92862013 153 for ( ievent=0; ievent<Nevents; ievent++)
154 {
155 if (ievent==0) cout << "AnalyzeManyEvents > " << "Starting Analyzing " << endl ;
156 //========== Create the Clusterizer
157 fClu = new AliPHOSClusterizerv1() ;
158 fClu->SetEmcEnergyThreshold(0.025) ;
159 fClu->SetEmcClusteringThreshold(0.75) ;
160 fClu->SetPpsdEnergyThreshold (0.0000002) ;
161 fClu->SetPpsdClusteringThreshold(0.0000001) ;
162 fClu->SetLocalMaxCut(0.03) ;
163 fClu->SetCalibrationParameters(0., 0.00000001) ;
164 //========== Creates the track segment maker
165 fTrs = new AliPHOSTrackSegmentMakerv1() ;
26d4b141 166 //========== Creates the particle identifier
167 fPID = new AliPHOSPIDv1() ;
09fc14a0 168 fPID->SetShowerProfileCuts(0.5, 1.5, 0.5, 1.5 ) ;
169 fPID->Print() ;
92862013 170 //========== Creates the Reconstructioner
26d4b141 171 fRec = new AliPHOSReconstructioner(fClu, fTrs, fPID) ;
92862013 172 //========== Event Number
173 if ( ( log10(ievent+1) - (Int_t)(log10(ievent+1)) ) == 0. ) cout << "AnalyzeManyEvents > " << "Event is " << ievent << endl ;
174 //=========== Connects the various Tree's for evt
175 gAlice->GetEvent(ievent);
176 //=========== Gets the Digit TTree
177 gAlice->TreeD()->GetEvent(0) ;
178 //=========== Gets the number of entries in the Digits array
179 TIter nextdigit(fPHOS->Digits()) ;
180 while( ( digit = (AliPHOSDigit *)nextdigit() ) )
181 {
182 fGeom->AbsToRelNumbering(digit->GetId(), relid) ;
183 if (fClu->IsInEmc(digit)) fhEmcDigit->Fill(fClu->Calibrate(digit->GetAmp())) ;
184 else
185 {
186 if (relid[1]<17) fhVetoDigit->Fill(fClu->Calibrate(digit->GetAmp()));
187 if (relid[1]>16) fhConvertorDigit->Fill(fClu->Calibrate(digit->GetAmp()));
188 }
189 }
190 //=========== Do the reconstruction
191 fPHOS->Reconstruction(fRec);
192 //=========== Cluster in module
193 TIter nextEmc(fPHOS->EmcClusters() ) ;
194 while((emc = (AliPHOSEmcRecPoint *)nextEmc()))
195 {
196 if ( emc->GetPHOSMod() == module )
197 {
198 fhEmcCluster->Fill( emc->GetTotalEnergy() );
199 TIter nextPpsd( fPHOS->PpsdClusters()) ;
200 while((ppsd = (AliPHOSPpsdRecPoint *)nextPpsd()))
201 {
202 if ( ppsd->GetPHOSMod() == module )
203 {
204 if (!ppsd->GetUp()) fhConvertorEmc->Fill(emc->GetTotalEnergy(),ppsd->GetTotalEnergy()) ;
205 }
206 }
207 }
208 }
209 //=========== Cluster in module PPSD Down
210 TIter nextPpsd(fPHOS->PpsdClusters() ) ;
211 while((ppsd = (AliPHOSPpsdRecPoint *)nextPpsd()))
212 {
213 if ( ppsd->GetPHOSMod() == module )
214 {
215 if (!ppsd->GetUp()) fhConvertorCluster->Fill(ppsd->GetTotalEnergy()) ;
216 if (ppsd->GetUp()) fhVetoCluster ->Fill(ppsd->GetTotalEnergy()) ;
217 }
218 }
219 //========== TRackSegments in the event
38f8ae6d 220 TIter nextRecParticle(fPHOS->RecParticles() ) ;
221 while((recparticle = (AliPHOSRecParticle *)nextRecParticle()))
92862013 222 {
38f8ae6d 223 if ( recparticle->GetPHOSTrackSegment()->GetPHOSMod() == module )
92862013 224 {
38f8ae6d 225 cout << "Particle type is " << recparticle->GetType() << endl ;
226 switch(recparticle->GetType())
92862013 227 {
09fc14a0 228 case kGAMMA:
38f8ae6d 229 fhPhotonEnergy->Fill(recparticle->Energy() ) ;
09fc14a0 230 //fhPhotonPositionX->Fill(recpart. ) ;
231 //fhPhotonPositionY->Fill(recpart. ) ;
38f8ae6d 232 cout << "PHOTON" << endl;
92862013 233 break;
09fc14a0 234 case kELECTRON:
38f8ae6d 235 fhElectronEnergy->Fill(recparticle->Energy() ) ;
09fc14a0 236 //fhElectronPositionX->Fill(recpart. ) ;
237 //fhElectronPositionY->Fill(recpart. ) ;
38f8ae6d 238 cout << "ELECTRON" << endl;
92862013 239 break;
b9bbdad1 240 case kNEUTRALHADRON:
38f8ae6d 241 fhNeutralHadronEnergy->Fill(recparticle->Energy() ) ;
b9bbdad1 242 //fhNeutralHadronPositionX->Fill(recpart. ) ;
243 //fhNeutralHadronPositionY->Fill(recpart. ) ;
38f8ae6d 244 cout << "NEUTRAl HADRON" << endl;
b9bbdad1 245 break ;
246 case kNEUTRALEM:
38f8ae6d 247 fhNeutralEMEnergy->Fill(recparticle->Energy() ) ;
b9bbdad1 248 //fhNeutralEMPositionX->Fill(recpart. ) ;
249 //fhNeutralEMPositionY->Fill(recpart. ) ;
250 //cout << "NEUTRAL EM" << endl;
92862013 251 break ;
0dd37dda 252 case kCHARGEDHADRON:
38f8ae6d 253 fhChargedHadronEnergy->Fill(recparticle->Energy() ) ;
09fc14a0 254 //fhChargedHadronPositionX->Fill(recpart. ) ;
255 //fhChargedHadronPositionY->Fill(recpart. ) ;
38f8ae6d 256 cout << "CHARGED HADRON" << endl;
92862013 257 break ;
258
259 }
260 }
261 }
26d4b141 262 // Deleting fClu, fTrs, fPID et fRec
92862013 263 fClu->Delete();
264 fTrs->Delete();
26d4b141 265 fPID->Delete();
92862013 266 fRec->Delete();
267
268 } // endfor
269 SavingHistograms();
270 } // endif
271} // endfunction
272
273
274//____________________________________________________________________________
275void AliPHOSAnalyze::BookingHistograms()
276{
b9bbdad1 277 if (fhEmcDigit )
278 delete fhEmcDigit ;
279 if (fhVetoDigit )
280 delete fhVetoDigit ;
281 if (fhConvertorDigit )
282 delete fhConvertorDigit ;
283 if (fhEmcCluster )
284 delete fhEmcCluster ;
285 if (fhVetoCluster )
286 delete fhVetoCluster ;
287 if (fhConvertorCluster )
288 delete fhConvertorCluster ;
289 if (fhConvertorEmc )
290 delete fhConvertorEmc ;
291
292 fhEmcDigit = new TH1F("hEmcDigit", "hEmcDigit", 1000, 0. , 25.);
293 fhVetoDigit = new TH1F("hVetoDigit", "hVetoDigit", 500, 0. , 3.e-5);
294 fhConvertorDigit = new TH1F("hConvertorDigit","hConvertorDigit", 500, 0. , 3.e-5);
295 fhEmcCluster = new TH1F("hEmcCluster", "hEmcCluster", 1000, 0. , 30.);
296 fhVetoCluster = new TH1F("hVetoCluster", "hVetoCluster", 500, 0. , 3.e-5);
297 fhConvertorCluster = new TH1F("hConvertorCluster","hConvertorCluster",500, 0. , 3.e-5);
298 fhConvertorEmc = new TH2F("hConvertorEmc", "hConvertorEmc", 200, 1. , 3., 200, 0., 3.e-5);
299 fhPhotonEnergy = new TH1F("hPhotonEnergy", "hPhotonEnergy", 1000, 0. , 30.);
300 fhElectronEnergy = new TH1F("hElectronEnergy","hElectronEnergy", 1000, 0. , 30.);
301 fhNeutralHadronEnergy = new TH1F("hNeutralHadronEnergy", "hNeutralHadronEnergy", 1000, 0. , 30.);
302 fhNeutralEMEnergy = new TH1F("hNeutralEMEnergy", "hNeutralEMEnergy", 1000, 0. , 30.);
303 fhChargedHadronEnergy = new TH1F("hChargedHadronEnergy", "hChargedHadronEnergy", 1000, 0. , 30.);
304 fhPhotonPositionX = new TH1F("hPhotonPositionX","hPhotonPositionX", 500,-80. , 80.);
305 fhElectronPositionX = new TH1F("hElectronPositionX","hElectronPositionX",500,-80. , 80.);
306 fhNeutralHadronPositionX = new TH1F("hNeutralHadronPositionX","hNeutralHadronPositionX",500,-80. , 80.);
307 fhNeutralEMPositionX = new TH1F("hNeutralEMPositionX","hNeutralEMPositionX",500,-80. , 80.);
09fc14a0 308 fhChargedHadronPositionX = new TH1F("hChargedHadronPositionX","hChargedHadronPositionX",500,-80. , 80.);
b9bbdad1 309 fhPhotonPositionY = new TH1F("hPhotonPositionY","hPhotonPositionY", 500,-80. , 80.);
310 fhElectronPositionY = new TH1F("hElectronPositionY","hElectronPositionY",500,-80. , 80.);
311 fhNeutralHadronPositionY = new TH1F("hNeutralHadronPositionY","hNeutralHadronPositionY",500,-80. , 80.);
312 fhNeutralEMPositionY = new TH1F("hNeutralEMPositionY","hNeutralEMPositionY",500,-80. , 80.);
09fc14a0 313 fhChargedHadronPositionY = new TH1F("hChargedHadronPositionY","hChargedHadronPositionY",500,-80. , 80.);
92862013 314
315}
6ad0bfa0 316//____________________________________________________________________________
317Bool_t AliPHOSAnalyze::Init(Int_t evt)
318{
319
92862013 320 Bool_t ok = kTRUE ;
6ad0bfa0 321
322 //========== Open galice root file
323
324 if ( fRootFile == 0 ) {
325 Text_t * name = new Text_t[80] ;
326 cout << "AnalyzeOneEvent > Enter file root file name : " ;
327 cin >> name ;
92862013 328 Bool_t ok = OpenRootFile(name) ;
329 if ( !ok )
6ad0bfa0 330 cout << " AliPHOSAnalyze > Error opening " << name << endl ;
331 else {
332 //========== Get AliRun object from file
333
334 gAlice = (AliRun*) fRootFile->Get("gAlice") ;
335
336 //=========== Get the PHOS object and associated geometry from the file
337
338 fPHOS = (AliPHOSv0 *)gAlice->GetDetector("PHOS") ;
339 fGeom = AliPHOSGeometry::GetInstance( fPHOS->GetGeometry()->GetName(), fPHOS->GetGeometry()->GetTitle() );
92862013 340 } // else !ok
6ad0bfa0 341 } // if fRootFile
342
92862013 343 if ( ok ) {
6ad0bfa0 344
345 //========== Create the Clusterizer
346
347 fClu = new AliPHOSClusterizerv1() ;
92862013 348 fClu->SetEmcEnergyThreshold(0.025) ;
349 fClu->SetEmcClusteringThreshold(0.75) ;
350 fClu->SetPpsdEnergyThreshold (0.0000002) ;
351 fClu->SetPpsdClusteringThreshold(0.0000001) ;
6ad0bfa0 352 fClu->SetLocalMaxCut(0.03) ;
92862013 353 fClu->SetCalibrationParameters(0., 0.00000001) ;
6ad0bfa0 354 cout << "AnalyzeOneEvent > using clusterizer " << fClu->GetName() << endl ;
355 fClu->PrintParameters() ;
356
357 //========== Creates the track segment maker
358
359 fTrs = new AliPHOSTrackSegmentMakerv1() ;
360 cout << "AnalyzeOneEvent > using tack segment maker " << fTrs->GetName() << endl ;
361
26d4b141 362 //========== Creates the particle identifier
6ad0bfa0 363
26d4b141 364 fPID = new AliPHOSPIDv1() ;
365 cout << "AnalyzeOneEvent > using particle identifier " << fPID->GetName() << endl ;
6ad0bfa0 366
367 //========== Creates the Reconstructioner
368
26d4b141 369 fRec = new AliPHOSReconstructioner(fClu, fTrs, fPID) ;
6ad0bfa0 370
371 //=========== Connect the various Tree's for evt
372
373 if ( evt == -999 ) {
374 cout << "AnalyzeOneEvent > Enter event number : " ;
375 cin >> evt ;
376 cout << evt << endl ;
377 }
378 fEvt = evt ;
379 gAlice->GetEvent(evt);
380
381 //=========== Get the Digit TTree
382
383 gAlice->TreeD()->GetEvent(0) ;
384
92862013 385 } // ok
6ad0bfa0 386
92862013 387 return ok ;
6ad0bfa0 388}
389
92862013 390
6ad0bfa0 391//____________________________________________________________________________
92862013 392void AliPHOSAnalyze::DisplayKineEvent(Int_t evt)
6ad0bfa0 393{
394 if (evt == -999)
395 evt = fEvt ;
396
92862013 397 Int_t module ;
6ad0bfa0 398 cout << "DisplayKineEvent > which module (1-5, -1: all) ? " ;
92862013 399 cin >> module ; cout << module << endl ;
6ad0bfa0 400
92862013 401 Int_t testparticle ;
6ad0bfa0 402 cout << " 22 : PHOTON " << endl
403 << " (-)11 : (POSITRON)ELECTRON " << endl
404 << " (-)2112 : (ANTI)NEUTRON " << endl
405 << " -999 : Everything else " << endl ;
406 cout << "DisplayKineEvent > enter PDG particle code to display " ;
92862013 407 cin >> testparticle ; cout << testparticle << endl ;
6ad0bfa0 408
92862013 409 Text_t histoname[80] ;
410 sprintf(histoname,"Event %d: Incident particles in module %d", evt, module) ;
6ad0bfa0 411
92862013 412 Double_t tm, tM, pm, pM ; // min and Max theta and phi covered by module
413 fGeom->EmcModuleCoverage(module, tm, tM, pm, pM, kDegre) ;
6ad0bfa0 414
415 Double_t theta, phi ;
416 fGeom->EmcXtalCoverage(theta, phi, kDegre) ;
417
418 Int_t tdim = (Int_t)( (tM - tm) / theta ) ;
419 Int_t pdim = (Int_t)( (pM - pm) / phi ) ;
420
421 tm -= theta ;
422 tM += theta ;
423 pm -= phi ;
424 pM += phi ;
425
92862013 426 TH2F * histoparticle = new TH2F("histoparticle", histoname,
6ad0bfa0 427 pdim, pm, pM, tdim, tm, tM) ;
92862013 428 histoparticle->SetStats(kFALSE) ;
6ad0bfa0 429
430 // Get pointers to Alice Particle TClonesArray
431
92862013 432 TParticle * particle;
433 TClonesArray * particlearray = gAlice->Particles();
6ad0bfa0 434
92862013 435 Text_t canvasname[80];
436 sprintf(canvasname,"Particles incident in PHOS/EMC module # %d",module) ;
437 TCanvas * kinecanvas = new TCanvas("kinecanvas", canvasname, 650, 500) ;
6ad0bfa0 438
439 // get the KINE Tree
440
92862013 441 TTree * kine = gAlice->TreeK() ;
442 Stat_t nParticles = kine->GetEntries() ;
443 cout << "DisplayKineEvent > events in kine " << nParticles << endl ;
6ad0bfa0 444
445 // loop over particles
446
92862013 447 Double_t kRADDEG = 180. / TMath::Pi() ;
6ad0bfa0 448 Int_t index1 ;
449 Int_t nparticlein = 0 ;
92862013 450 for (index1 = 0 ; index1 < nParticles ; index1++){
451 Int_t nparticle = particlearray->GetEntriesFast() ;
6ad0bfa0 452 Int_t index2 ;
453 for( index2 = 0 ; index2 < nparticle ; index2++) {
92862013 454 particle = (TParticle*)particlearray->UncheckedAt(index2) ;
455 Int_t particletype = particle->GetPdgCode() ;
456 if (testparticle == -999 || testparticle == particletype) {
457 Double_t phi = particle->Phi() ;
458 Double_t theta = particle->Theta() ;
6ad0bfa0 459 Int_t mod ;
460 Double_t x, z ;
92862013 461 fGeom->ImpactOnEmc(theta, phi, mod, z, x) ;
462 if ( mod == module ) {
6ad0bfa0 463 nparticlein++ ;
92862013 464 histoparticle->Fill(phi*kRADDEG, theta*kRADDEG, particle->Energy() ) ;
6ad0bfa0 465 }
466 }
467 }
468 }
92862013 469 kinecanvas->Draw() ;
470 histoparticle->Draw("color") ;
471 TPaveText * pavetext = new TPaveText(294, 100, 300, 101);
6ad0bfa0 472 Text_t text[40] ;
473 sprintf(text, "Particles: %d ", nparticlein) ;
92862013 474 pavetext->AddText(text) ;
475 pavetext->Draw() ;
476 kinecanvas->Update();
6ad0bfa0 477
478}
479//____________________________________________________________________________
480void AliPHOSAnalyze::DisplayRecParticles()
481{
482 if (fEvt == -999) {
483 cout << "DisplayRecPoints > Analyze an event first ... (y/n) " ;
484 Text_t answer[1] ;
485 cin >> answer ; cout << answer ;
486 if ( answer == "y" )
487 AnalyzeOneEvent() ;
488 }
489 if (fEvt != -999) {
490
92862013 491 Int_t module ;
6ad0bfa0 492 cout << "DisplayRecPoints > which module (1-5, -1: all) ? " ;
92862013 493 cin >> module ; cout << module << endl ;
494 Text_t histoname[80] ;
495 sprintf(histoname,"Event %d: Reconstructed particles in module %d", fEvt, module) ;
496 Double_t tm, tM, pm, pM ; // min and Max theta and phi covered by module
497 fGeom->EmcModuleCoverage(module, tm, tM, pm, pM, kDegre) ;
6ad0bfa0 498 Double_t theta, phi ;
499 fGeom->EmcXtalCoverage(theta, phi, kDegre) ;
500 Int_t tdim = (Int_t)( (tM - tm) / theta ) ;
501 Int_t pdim = (Int_t)( (pM - pm) / phi ) ;
502 tm -= theta ;
503 tM += theta ;
504 pm -= phi ;
92862013 505 TH2F * histoRparticle = new TH2F("histoRparticle", histoname,
6ad0bfa0 506 pdim, pm, pM, tdim, tm, tM) ;
92862013 507 histoRparticle->SetStats(kFALSE) ;
508 Text_t canvasname[80] ;
509 sprintf(canvasname, "Reconstructed particles in PHOSmodule # %d", module) ;
510 TCanvas * rparticlecanvas = new TCanvas("RparticleCanvas", canvasname, 650, 500) ;
6ad0bfa0 511 RecParticlesList * rpl = fPHOS->RecParticles() ;
92862013 512 Int_t nRecParticles = rpl->GetEntries() ;
513 Int_t nRecParticlesInModule = 0 ;
6ad0bfa0 514 TIter nextRecPart(rpl) ;
515 AliPHOSRecParticle * rp ;
92862013 516 cout << "DisplayRecParticles > " << nRecParticles << " reconstructed particles " << endl ;
517 Double_t kRADDEG = 180. / TMath::Pi() ;
6ad0bfa0 518 while ( (rp = (AliPHOSRecParticle *)nextRecPart() ) ) {
519 AliPHOSTrackSegment * ts = rp->GetPHOSTrackSegment() ;
92862013 520 if ( ts->GetPHOSMod() == module ) {
521 nRecParticlesInModule++ ;
522 Double_t theta = rp->Theta() * kRADDEG ;
523 Double_t phi = rp->Phi() * kRADDEG ;
6ad0bfa0 524 Double_t energy = rp->Energy() ;
92862013 525 histoRparticle->Fill(phi, theta, energy) ;
6ad0bfa0 526 }
527 }
92862013 528 histoRparticle->Draw("color") ;
6ad0bfa0 529 Text_t text[80] ;
92862013 530 sprintf(text, "reconstructed particles: %d", nRecParticlesInModule) ;
531 TPaveText * pavetext = new TPaveText(292, 100, 300, 101);
532 pavetext->AddText(text) ;
533 pavetext->Draw() ;
534 rparticlecanvas->Update() ;
6ad0bfa0 535 }
536}
537
538//____________________________________________________________________________
539void AliPHOSAnalyze::DisplayRecPoints()
540{
541 if (fEvt == -999) {
542 cout << "DisplayRecPoints > Analyze an event first ... (y/n) " ;
543 Text_t answer[1] ;
544 cin >> answer ; cout << answer ;
545 if ( answer == "y" )
546 AnalyzeOneEvent() ;
547 }
548 if (fEvt != -999) {
549
92862013 550 Int_t module ;
6ad0bfa0 551 cout << "DisplayRecPoints > which module (1-5, -1: all) ? " ;
92862013 552 cin >> module ; cout << module << endl ;
6ad0bfa0 553
92862013 554 Text_t canvasname[80];
555 sprintf(canvasname,"Digits in PHOS/EMC module # %d",module) ;
556 TCanvas * modulecanvas = new TCanvas("module", canvasname, 650, 500) ;
557 modulecanvas->Draw() ;
6ad0bfa0 558
92862013 559 //=========== Creating 2d-histogram of the PHOS module
6ad0bfa0 560 // a little bit junkie but is used to test Geom functinalities
561
92862013 562 Double_t tm, tM, pm, pM ; // min and Max theta and phi covered by module
6ad0bfa0 563
92862013 564 fGeom->EmcModuleCoverage(module, tm, tM, pm, pM);
6ad0bfa0 565 // convert angles into coordinates local to the EMC module of interest
566
92862013 567 Int_t emcModuleNumber ;
568 Double_t emcModulexm, emcModulezm ; // minimum local coordinate in a given EMCA module
569 Double_t emcModulexM, emcModulezM ; // maximum local coordinate in a given EMCA module
570 fGeom->ImpactOnEmc(tm, pm, emcModuleNumber, emcModulezm, emcModulexm) ;
571 fGeom->ImpactOnEmc(tM, pM, emcModuleNumber, emcModulezM, emcModulexM) ;
572 Int_t xdim = (Int_t)( ( emcModulexM - emcModulexm ) / fGeom->GetCrystalSize(0) ) ;
573 Int_t zdim = (Int_t)( ( emcModulezM - emcModulezm ) / fGeom->GetCrystalSize(2) ) ;
574 Float_t xmin = emcModulexm - fGeom->GetCrystalSize(0) ;
575 Float_t xMax = emcModulexM + fGeom->GetCrystalSize(0) ;
576 Float_t zmin = emcModulezm - fGeom->GetCrystalSize(2) ;
577 Float_t zMax = emcModulezM + fGeom->GetCrystalSize(2) ;
578 Text_t histoname[80];
579 sprintf(histoname,"Event %d: Digits and RecPoints in module %d", fEvt, module) ;
580 TH2F * hModule = new TH2F("HistoReconstructed", histoname,
6ad0bfa0 581 xdim, xmin, xMax, zdim, zmin, zMax) ;
582 hModule->SetMaximum(2.0);
583 hModule->SetMinimum(0.0);
584 hModule->SetStats(kFALSE);
585
586 TIter next(fPHOS->Digits()) ;
92862013 587 Float_t energy, y, z;
588 Float_t etot=0.;
589 Int_t relid[4]; Int_t nDigits = 0 ;
6ad0bfa0 590 AliPHOSDigit * digit ;
6ad0bfa0 591 while((digit = (AliPHOSDigit *)next()))
592 {
92862013 593 fGeom->AbsToRelNumbering(digit->GetId(), relid) ;
594 if (relid[0] == module)
6ad0bfa0 595 {
92862013 596 nDigits++ ;
597 energy = fClu->Calibrate(digit->GetAmp()) ;
598 etot += energy ;
599 fGeom->RelPosInModule(relid,y,z) ;
600 if (energy > 0.01 )
601 hModule->Fill(y, z, energy) ;
6ad0bfa0 602 }
603 }
92862013 604 cout <<"DrawRecPoints > Found in module "
605 << module << " " << nDigits << " digits with total energy " << etot << endl ;
6ad0bfa0 606 hModule->Draw("col2") ;
607
92862013 608 //=========== Cluster in module
6ad0bfa0 609
92862013 610 TClonesArray * emcRP = fPHOS->EmcClusters() ;
611 etot = 0.;
612 Int_t totalnClusters = 0 ;
613 Int_t nClusters = 0 ;
614 TIter nextemc(emcRP) ;
6ad0bfa0 615 AliPHOSEmcRecPoint * emc ;
616 while((emc = (AliPHOSEmcRecPoint *)nextemc()))
617 {
cf239357 618 Int_t numberofprimaries ;
619 Int_t * primariesarray = new Int_t[10] ;
620 emc->GetPrimaries(numberofprimaries, primariesarray) ;
92862013 621 totalnClusters++ ;
622 if ( emc->GetPHOSMod() == module )
6ad0bfa0 623 {
92862013 624 nClusters++ ;
625 energy = emc->GetTotalEnergy() ;
626 etot+= energy ;
26d4b141 627 emc->Draw("M") ;
6ad0bfa0 628 }
629 }
92862013 630 cout << "DrawRecPoints > Found " << totalnClusters << " EMC Clusters in PHOS" << endl ;
631 cout << "DrawRecPoints > Found in module " << module << " " << nClusters << " EMC Clusters " << endl ;
632 cout << "DrawRecPoints > total energy " << etot << endl ;
6ad0bfa0 633
92862013 634 TPaveText * pavetext = new TPaveText(22, 80, 83, 90);
6ad0bfa0 635 Text_t text[40] ;
92862013 636 sprintf(text, "digits: %d; clusters: %d", nDigits, nClusters) ;
637 pavetext->AddText(text) ;
638 pavetext->Draw() ;
639 modulecanvas->Update();
6ad0bfa0 640
92862013 641 //=========== Cluster in module PPSD Down
6ad0bfa0 642
92862013 643 TClonesArray * ppsdRP = fPHOS->PpsdClusters() ;
644 etot = 0.;
645 TIter nextPpsd(ppsdRP) ;
646 AliPHOSPpsdRecPoint * ppsd ;
647 while((ppsd = (AliPHOSPpsdRecPoint *)nextPpsd()))
6ad0bfa0 648 {
92862013 649 totalnClusters++ ;
650 if ( ppsd->GetPHOSMod() == module )
6ad0bfa0 651 {
92862013 652 nClusters++ ;
653 energy = ppsd->GetEnergy() ;
654 etot+=energy ;
655 if (!ppsd->GetUp()) ppsd->Draw("P") ;
6ad0bfa0 656 }
657 }
92862013 658 cout << "DrawRecPoints > Found " << totalnClusters << " Ppsd Down Clusters in PHOS" << endl ;
659 cout << "DrawRecPoints > Found in module " << module << " " << nClusters << " Ppsd Down Clusters " << endl ;
660 cout << "DrawRecPoints > total energy " << etot << endl ;
6ad0bfa0 661
92862013 662 //=========== Cluster in module PPSD Up
6ad0bfa0 663
92862013 664 ppsdRP = fPHOS->PpsdClusters() ;
665 etot = 0.;
666 TIter nextPpsdUp(ppsdRP) ;
667 while((ppsd = (AliPHOSPpsdRecPoint *)nextPpsdUp()))
6ad0bfa0 668 {
92862013 669 totalnClusters++ ;
670 if ( ppsd->GetPHOSMod() == module )
6ad0bfa0 671 {
92862013 672 nClusters++ ;
673 energy = ppsd->GetEnergy() ;
674 etot+=energy ;
675 if (ppsd->GetUp()) ppsd->Draw("P") ;
6ad0bfa0 676 }
677 }
92862013 678 cout << "DrawRecPoints > Found " << totalnClusters << " Ppsd Up Clusters in PHOS" << endl ;
679 cout << "DrawRecPoints > Found in module " << module << " " << nClusters << " Ppsd Up Clusters " << endl ;
680 cout << "DrawRecPoints > total energy " << etot << endl ;
6ad0bfa0 681
682 } // if !-999
683}
684
685//____________________________________________________________________________
686void AliPHOSAnalyze::DisplayTrackSegments()
687{
688 if (fEvt == -999) {
689 cout << "DisplayTrackSegments > Analyze an event first ... (y/n) " ;
690 Text_t answer[1] ;
691 cin >> answer ; cout << answer ;
692 if ( answer == "y" )
693 AnalyzeOneEvent() ;
694 }
695 if (fEvt != -999) {
696
92862013 697 Int_t module ;
6ad0bfa0 698 cout << "DisplayTrackSegments > which module (1-5, -1: all) ? " ;
92862013 699 cin >> module ; cout << module << endl ;
700 //=========== Creating 2d-histogram of the PHOS module
6ad0bfa0 701 // a little bit junkie but is used to test Geom functinalities
702
92862013 703 Double_t tm, tM, pm, pM ; // min and Max theta and phi covered by module
6ad0bfa0 704
92862013 705 fGeom->EmcModuleCoverage(module, tm, tM, pm, pM);
6ad0bfa0 706 // convert angles into coordinates local to the EMC module of interest
707
92862013 708 Int_t emcModuleNumber ;
709 Double_t emcModulexm, emcModulezm ; // minimum local coordinate in a given EMCA module
710 Double_t emcModulexM, emcModulezM ; // maximum local coordinate in a given EMCA module
711 fGeom->ImpactOnEmc(tm, pm, emcModuleNumber, emcModulezm, emcModulexm) ;
712 fGeom->ImpactOnEmc(tM, pM, emcModuleNumber, emcModulezM, emcModulexM) ;
713 Int_t xdim = (Int_t)( ( emcModulexM - emcModulexm ) / fGeom->GetCrystalSize(0) ) ;
714 Int_t zdim = (Int_t)( ( emcModulezM - emcModulezm ) / fGeom->GetCrystalSize(2) ) ;
715 Float_t xmin = emcModulexm - fGeom->GetCrystalSize(0) ;
716 Float_t xMax = emcModulexM + fGeom->GetCrystalSize(0) ;
717 Float_t zmin = emcModulezm - fGeom->GetCrystalSize(2) ;
718 Float_t zMax = emcModulezM + fGeom->GetCrystalSize(2) ;
719 Text_t histoname[80];
720 sprintf(histoname,"Event %d: Track Segments in module %d", fEvt, module) ;
721 TH2F * histotrack = new TH2F("histotrack", histoname,
6ad0bfa0 722 xdim, xmin, xMax, zdim, zmin, zMax) ;
92862013 723 histotrack->SetStats(kFALSE);
724 Text_t canvasname[80];
725 sprintf(canvasname,"Track segments in PHOS/EMC-PPSD module # %d", module) ;
726 TCanvas * trackcanvas = new TCanvas("TrackSegmentCanvas", canvasname, 650, 500) ;
727 histotrack->Draw() ;
6ad0bfa0 728
729 TrackSegmentsList * trsegl = fPHOS->TrackSegments() ;
730 AliPHOSTrackSegment * trseg ;
731
92862013 732 Int_t nTrackSegments = trsegl->GetEntries() ;
6ad0bfa0 733 Int_t index ;
92862013 734 Float_t etot = 0 ;
735 Int_t nTrackSegmentsInModule = 0 ;
736 for(index = 0; index < nTrackSegments ; index++){
6ad0bfa0 737 trseg = (AliPHOSTrackSegment * )trsegl->At(index) ;
92862013 738 etot+= trseg->GetEnergy() ;
739 if ( trseg->GetPHOSMod() == module ) {
740 nTrackSegmentsInModule++ ;
6ad0bfa0 741 trseg->Draw("P");
742 }
743 }
744 Text_t text[80] ;
92862013 745 sprintf(text, "track segments: %d", nTrackSegmentsInModule) ;
746 TPaveText * pavetext = new TPaveText(22, 80, 83, 90);
747 pavetext->AddText(text) ;
748 pavetext->Draw() ;
749 trackcanvas->Update() ;
750 cout << "DisplayTrackSegments > Found " << trsegl->GetEntries() << " Track segments with total energy "<< etot << endl ;
6ad0bfa0 751
752 }
753}
754//____________________________________________________________________________
755Bool_t AliPHOSAnalyze::OpenRootFile(Text_t * name)
756{
757 fRootFile = new TFile(name) ;
758 return fRootFile->IsOpen() ;
759}
92862013 760//____________________________________________________________________________
761void AliPHOSAnalyze::SavingHistograms()
762{
3862b681 763 Text_t outputname[80] ;// = fRootFile->GetName();
764 sprintf(outputname,"%s.analyzed",fRootFile->GetName());
765 TFile output(outputname,"RECREATE");
92862013 766 output.cd();
b9bbdad1 767 if (fhEmcDigit )
768 fhEmcDigit->Write() ;
769 if (fhVetoDigit )
770 fhVetoDigit->Write() ;
771 if (fhConvertorDigit )
772 fhConvertorDigit->Write() ;
773 if (fhEmcCluster )
774 fhEmcCluster->Write() ;
775 if (fhVetoCluster )
776 fhVetoCluster->Write() ;
777 if (fhConvertorCluster )
778 fhConvertorCluster->Write() ;
779 if (fhConvertorEmc )
780 fhConvertorEmc->Write() ;
781 if (fhPhotonEnergy)
782 fhPhotonEnergy->Write() ;
783 if (fhPhotonPositionX)
784 fhPhotonPositionX->Write() ;
785 if (fhPhotonPositionY)
786 fhPhotonPositionX->Write() ;
787 if (fhElectronEnergy)
788 fhElectronEnergy->Write() ;
789 if (fhElectronPositionX)
790 fhElectronPositionX->Write() ;
791 if (fhElectronPositionY)
792 fhElectronPositionX->Write() ;
793 if (fhNeutralHadronEnergy)
794 fhNeutralHadronEnergy->Write() ;
795 if (fhNeutralHadronPositionX)
796 fhNeutralHadronPositionX->Write() ;
797 if (fhNeutralHadronPositionY)
798 fhNeutralHadronPositionX->Write() ;
799 if (fhNeutralEMEnergy)
800 fhNeutralEMEnergy->Write() ;
801 if (fhNeutralEMPositionX)
802 fhNeutralEMPositionX->Write() ;
803 if (fhNeutralEMPositionY)
804 fhNeutralEMPositionX->Write() ;
805 if (fhChargedHadronEnergy)
806 fhChargedHadronEnergy->Write() ;
807 if (fhChargedHadronPositionX)
808 fhChargedHadronPositionX->Write() ;
809 if (fhChargedHadronPositionY)
810 fhChargedHadronPositionX->Write() ;
92862013 811
812 output.Write();
813 output.Close();
814}