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