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