1 /**************************************************************************
2 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
4 * Author: The ALICE Off-line Project. *
5 * Contributors are mentioned in the code where appropriate. *
7 * Permission to use, copy, modify and distribute this software and its *
8 * documentation strictly for non-commercial purposes is hereby granted *
9 * without fee, provided that the above copyright notice appears in all *
10 * copies and that both the copyright notice and this permission notice *
11 * appear in the supporting documentation. The authors make no claims *
12 * about the suitability of this software for any purpose. It is *
13 * provided "as is" without express or implied warranty. *
14 **************************************************************************/
18 //_________________________________________________________________________
19 // Algorythm class to analyze PHOSv0 events:
20 // Construct histograms and displays them.
21 // Use the macro EditorBar.C for best access to the functionnalities
23 //*-- Author: Y. Schutz (SUBATECH)
24 //////////////////////////////////////////////////////////////////////////////
26 // --- ROOT system ---
32 #include "TParticle.h"
33 #include "TClonesArray.h"
38 // --- Standard library ---
43 // --- AliRoot header files ---
46 #include "AliPHOSAnalyze.h"
47 #include "AliPHOSClusterizerv1.h"
48 #include "AliPHOSTrackSegmentMakerv1.h"
49 #include "AliPHOSPIDv1.h"
50 #include "AliPHOSReconstructioner.h"
51 #include "AliPHOSDigit.h"
52 #include "AliPHOSTrackSegment.h"
53 #include "AliPHOSRecParticle.h"
54 #include "AliPHOSIndexToObject.h"
56 ClassImp(AliPHOSAnalyze)
59 //____________________________________________________________________________
60 AliPHOSAnalyze::AliPHOSAnalyze()
62 // default ctor (useless)
67 //____________________________________________________________________________
68 AliPHOSAnalyze::AliPHOSAnalyze(Text_t * name)
70 // ctor: analyze events from root file "name"
72 Bool_t ok = OpenRootFile(name) ;
74 cout << " AliPHOSAnalyze > Error opening " << name << endl ;
77 gAlice = (AliRun*) fRootFile->Get("gAlice");
78 fPHOS = (AliPHOSv0 *)gAlice->GetDetector("PHOS") ;
79 fGeom = AliPHOSGeometry::GetInstance( fPHOS->GetGeometry()->GetName(), fPHOS->GetGeometry()->GetTitle() ) ;
84 //____________________________________________________________________________
85 AliPHOSAnalyze::~AliPHOSAnalyze()
89 if (fRootFile->IsOpen() )
111 //____________________________________________________________________________
112 void AliPHOSAnalyze::AnalyzeOneEvent(Int_t evt)
114 // analyze one single event with id=evt
118 Bool_t ok = Init(evt) ;
121 //=========== Get the number of entries in the Digits array
123 Int_t nId = fPHOS->Digits()->GetEntries();
124 printf("AnalyzeOneEvent > Number of entries in the Digit array is %d \n",nId);
126 //=========== Do the reconstruction
128 cout << "AnalyzeOneEvent > Found " << nId << " digits in PHOS" << endl ;
131 fPHOS->Reconstruction(fRec);
133 // =========== End of reconstruction
135 // =========== Write the root file
139 // =========== Finish
141 cout << "AnalyzeOneEvent > event # " << fEvt << " processed" << endl ;
144 cout << "AnalyzeOneEvent > filed to process event # " << evt << endl ;
146 ts.Stop() ; cout << "CPU time = " << ts.CpuTime() << endl ;
147 cout << "Real time = " << ts.RealTime() << endl ;
150 //____________________________________________________________________________
151 void AliPHOSAnalyze::AnalyzeManyEvents(Int_t Nevents, Int_t module)
153 // analyzes Nevents events in a single PHOS module
155 if ( fRootFile == 0 )
156 cout << "AnalyzeManyEvents > " << "Root File not openned" << endl ;
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 ;
169 AliPHOSDigit * digit ;
170 AliPHOSEmcRecPoint * emc ;
171 AliPHOSPpsdRecPoint * ppsd ;
172 // AliPHOSTrackSegment * tracksegment ;
173 AliPHOSRecParticle * recparticle;
174 for ( ievent=0; ievent<Nevents; ievent++)
176 if (ievent==0) cout << "AnalyzeManyEvents > " << "Starting Analyzing " << endl ;
177 //========== Create the Clusterizer
178 fClu = new AliPHOSClusterizerv1() ;
179 fClu->SetEmcEnergyThreshold(0.025) ;
180 fClu->SetEmcClusteringThreshold(0.50) ;
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
186 fTrs = new AliPHOSTrackSegmentMakerv1() ;
187 fTrs->UnsetUnfoldFlag() ;
188 //========== Creates the particle identifier
189 fPID = new AliPHOSPIDv1() ;
190 fPID->SetShowerProfileCuts(0.3, 1.8, 0.3, 1.8 ) ;
192 //========== Creates the Reconstructioner
193 fRec = new AliPHOSReconstructioner(fClu, fTrs, fPID) ;
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 ;
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() ) )
205 fGeom->AbsToRelNumbering(digit->GetId(), relid) ;
206 if (fClu->IsInEmc(digit)) fhEmcDigit->Fill(fClu->Calibrate(digit->GetAmp())) ;
209 if (relid[1]<17) fhVetoDigit->Fill(fClu->Calibrate(digit->GetAmp()));
210 if (relid[1]>16) fhConvertorDigit->Fill(fClu->Calibrate(digit->GetAmp()));
213 //=========== Do the reconstruction
214 fPHOS->Reconstruction(fRec);
215 //=========== Cluster in module
216 TIter nextEmc(fPHOS->EmcRecPoints() ) ;
217 while((emc = (AliPHOSEmcRecPoint *)nextEmc()))
219 if ( emc->GetPHOSMod() == module )
221 fhEmcCluster->Fill( emc->GetTotalEnergy() );
222 TIter nextPpsd( fPHOS->PpsdRecPoints()) ;
223 while((ppsd = (AliPHOSPpsdRecPoint *)nextPpsd()))
225 if ( ppsd->GetPHOSMod() == module )
227 if (!ppsd->GetUp()) fhConvertorEmc->Fill(emc->GetTotalEnergy(),ppsd->GetTotalEnergy()) ;
232 //=========== Cluster in module PPSD Down
233 TIter nextPpsd(fPHOS->PpsdRecPoints() ) ;
234 while((ppsd = (AliPHOSPpsdRecPoint *)nextPpsd()))
236 if ( ppsd->GetPHOSMod() == module )
238 if (!ppsd->GetUp()) fhConvertorCluster->Fill(ppsd->GetTotalEnergy()) ;
239 if (ppsd->GetUp()) fhVetoCluster ->Fill(ppsd->GetTotalEnergy()) ;
242 //========== TRackSegments in the event
243 TIter nextRecParticle(fPHOS->RecParticles() ) ;
244 while((recparticle = (AliPHOSRecParticle *)nextRecParticle()))
246 if ( recparticle->GetPHOSTrackSegment()->GetPHOSMod() == module )
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 ;
253 for ( index = 0 ; index < numberofprimaries ; index++)
254 cout << " primary # " << index << " = " << listofprimaries[index] << endl ;
255 switch(recparticle->GetType())
258 fhPhotonEnergy->Fill(recparticle->Energy() ) ;
259 //fhPhotonPositionX->Fill(recpart. ) ;
260 //fhPhotonPositionY->Fill(recpart. ) ;
261 cout << "PHOTON" << endl;
264 fhElectronEnergy->Fill(recparticle->Energy() ) ;
265 //fhElectronPositionX->Fill(recpart. ) ;
266 //fhElectronPositionY->Fill(recpart. ) ;
267 cout << "ELECTRON" << endl;
270 fhNeutralHadronEnergy->Fill(recparticle->Energy() ) ;
271 //fhNeutralHadronPositionX->Fill(recpart. ) ;
272 //fhNeutralHadronPositionY->Fill(recpart. ) ;
273 cout << "NEUTRAl HADRON" << endl;
276 fhNeutralEMEnergy->Fill(recparticle->Energy() ) ;
277 //fhNeutralEMPositionX->Fill(recpart. ) ;
278 //fhNeutralEMPositionY->Fill(recpart. ) ;
279 //cout << "NEUTRAL EM" << endl;
282 fhChargedHadronEnergy->Fill(recparticle->Energy() ) ;
283 //fhChargedHadronPositionX->Fill(recpart. ) ;
284 //fhChargedHadronPositionY->Fill(recpart. ) ;
285 cout << "CHARGED HADRON" << endl;
288 fhPhotonHadronEnergy->Fill(recparticle->Energy() ) ;
289 //fhPhotonHadronPositionX->Fill(recpart. ) ;
290 //fhPhotonHadronPositionY->Fill(recpart. ) ;
291 cout << "PHOTON HADRON" << endl;
296 // Deleting fClu, fTrs, fPID et fRec
308 //____________________________________________________________________________
309 void AliPHOSAnalyze::BookingHistograms()
311 // Books the histograms where the results of the analysis are stored (to be changed)
317 if (fhConvertorDigit )
318 delete fhConvertorDigit ;
320 delete fhEmcCluster ;
322 delete fhVetoCluster ;
323 if (fhConvertorCluster )
324 delete fhConvertorCluster ;
326 delete fhConvertorEmc ;
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.);
340 fhPhotonHadronEnergy = new TH1F("hPhotonHadronEnergy","hPhotonHadronEnergy",500,-80. , 80.);
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.);
345 fhChargedHadronPositionX = new TH1F("hChargedHadronPositionX","hChargedHadronPositionX",500,-80. , 80.);
346 fhPhotonHadronPositionX = new TH1F("hPhotonHadronPositionX","hPhotonHadronPositionX",500,-80. , 80.);
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.);
351 fhChargedHadronPositionY = new TH1F("hChargedHadronPositionY","hChargedHadronPositionY",500,-80. , 80.);
352 fhPhotonHadronPositionY = new TH1F("hPhotonHadronPositionY","hPhotonHadronPositionY",500,-80. , 80.);
356 //____________________________________________________________________________
357 Bool_t AliPHOSAnalyze::Init(Int_t evt)
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
366 //========== Open galice root file
368 if ( fRootFile == 0 ) {
369 Text_t * name = new Text_t[80] ;
370 cout << "AnalyzeOneEvent > Enter file root file name : " ;
372 Bool_t ok = OpenRootFile(name) ;
374 cout << " AliPHOSAnalyze > Error opening " << name << endl ;
376 //========== Get AliRun object from file
378 gAlice = (AliRun*) fRootFile->Get("gAlice") ;
380 //=========== Get the PHOS object and associated geometry from the file
382 fPHOS = (AliPHOSv0 *)gAlice->GetDetector("PHOS") ;
383 fGeom = AliPHOSGeometry::GetInstance( fPHOS->GetGeometry()->GetName(), fPHOS->GetGeometry()->GetTitle() );
391 //========== Initializes the Index to Object converter
393 fObjGetter = AliPHOSIndexToObject::GetInstance(fPHOS) ;
395 //========== Create the Clusterizer
397 fClu = new AliPHOSClusterizerv1() ;
398 fClu->SetEmcEnergyThreshold(0.030) ;
399 fClu->SetEmcClusteringThreshold(1.0) ;
400 fClu->SetPpsdEnergyThreshold (0.0000002) ;
401 fClu->SetPpsdClusteringThreshold(0.0000001) ;
402 fClu->SetLocalMaxCut(0.03) ;
403 fClu->SetCalibrationParameters(0., 0.00000001) ;
404 cout << "AnalyzeOneEvent > using clusterizer " << fClu->GetName() << endl ;
405 fClu->PrintParameters() ;
407 //========== Creates the track segment maker
409 fTrs = new AliPHOSTrackSegmentMakerv1() ;
410 cout << "AnalyzeOneEvent > using tack segment maker " << fTrs->GetName() << endl ;
411 fTrs->UnsetUnfoldFlag() ;
413 //========== Creates the particle identifier
415 fPID = new AliPHOSPIDv1() ;
416 cout << "AnalyzeOneEvent > using particle identifier " << fPID->GetName() << endl ;
418 //========== Creates the Reconstructioner
420 fRec = new AliPHOSReconstructioner(fClu, fTrs, fPID) ;
422 //=========== Connect the various Tree's for evt
425 cout << "AnalyzeOneEvent > Enter event number : " ;
427 cout << evt << endl ;
430 gAlice->GetEvent(evt);
432 //=========== Get the Digit TTree
434 gAlice->TreeD()->GetEvent(0) ;
442 //____________________________________________________________________________
443 void AliPHOSAnalyze::DisplayKineEvent(Int_t evt)
445 // Display particles from the Kine Tree in global Alice (theta, phi) coordinates.
446 // One PHOS module at the time.
447 // The particle type can be selected.
453 cout << "DisplayKineEvent > which module (1-5, -1: all) ? " ;
454 cin >> module ; cout << module << endl ;
457 cout << " 22 : PHOTON " << endl
458 << " (-)11 : (POSITRON)ELECTRON " << endl
459 << " (-)2112 : (ANTI)NEUTRON " << endl
460 << " -999 : Everything else " << endl ;
461 cout << "DisplayKineEvent > enter PDG particle code to display " ;
462 cin >> testparticle ; cout << testparticle << endl ;
464 Text_t histoname[80] ;
465 sprintf(histoname,"Event %d: Incident particles in module %d", evt, module) ;
467 Double_t tm, tM, pm, pM ; // min and Max theta and phi covered by module
468 fGeom->EmcModuleCoverage(module, tm, tM, pm, pM, kDegre) ;
470 Double_t theta, phi ;
471 fGeom->EmcXtalCoverage(theta, phi, kDegre) ;
473 Int_t tdim = (Int_t)( (tM - tm) / theta ) ;
474 Int_t pdim = (Int_t)( (pM - pm) / phi ) ;
481 TH2F * histoparticle = new TH2F("histoparticle", histoname,
482 pdim, pm, pM, tdim, tm, tM) ;
483 histoparticle->SetStats(kFALSE) ;
485 // Get pointers to Alice Particle TClonesArray
487 TParticle * particle;
488 TClonesArray * particlearray = gAlice->Particles();
490 Text_t canvasname[80];
491 sprintf(canvasname,"Particles incident in PHOS/EMC module # %d",module) ;
492 TCanvas * kinecanvas = new TCanvas("kinecanvas", canvasname, 650, 500) ;
496 TTree * kine = gAlice->TreeK() ;
497 Stat_t nParticles = kine->GetEntries() ;
498 cout << "DisplayKineEvent > events in kine " << nParticles << endl ;
500 // loop over particles
502 Double_t kRADDEG = 180. / TMath::Pi() ;
504 Int_t nparticlein = 0 ;
505 for (index1 = 0 ; index1 < nParticles ; index1++){
506 Int_t nparticle = particlearray->GetEntriesFast() ;
508 for( index2 = 0 ; index2 < nparticle ; index2++) {
509 particle = (TParticle*)particlearray->UncheckedAt(index2) ;
510 Int_t particletype = particle->GetPdgCode() ;
511 if (testparticle == -999 || testparticle == particletype) {
512 Double_t phi = particle->Phi() ;
513 Double_t theta = particle->Theta() ;
516 fGeom->ImpactOnEmc(theta, phi, mod, z, x) ;
517 if ( mod == module ) {
519 if (particle->Energy() > fClu->GetEmcClusteringThreshold() )
520 histoparticle->Fill(phi*kRADDEG, theta*kRADDEG, particle->Energy() ) ;
526 histoparticle->Draw("color") ;
527 TPaveText * pavetext = new TPaveText(294, 100, 300, 101);
529 sprintf(text, "Particles: %d ", nparticlein) ;
530 pavetext->AddText(text) ;
532 kinecanvas->Update();
535 //____________________________________________________________________________
536 void AliPHOSAnalyze::DisplayRecParticles()
538 // Display reconstructed particles in global Alice(theta, phi) coordinates.
539 // One PHOS module at the time.
540 // Click on symbols indicate the reconstructed particle type.
543 cout << "DisplayRecParticles > Analyze an event first ... (y/n) " ;
545 cin >> answer ; cout << answer ;
552 cout << "DisplayRecParticles > which module (1-5, -1: all) ? " ;
553 cin >> module ; cout << module << endl ;
554 Text_t histoname[80] ;
555 sprintf(histoname,"Event %d: Reconstructed particles in module %d", fEvt, module) ;
556 Double_t tm, tM, pm, pM ; // min and Max theta and phi covered by module
557 fGeom->EmcModuleCoverage(module, tm, tM, pm, pM, kDegre) ;
558 Double_t theta, phi ;
559 fGeom->EmcXtalCoverage(theta, phi, kDegre) ;
560 Int_t tdim = (Int_t)( (tM - tm) / theta ) ;
561 Int_t pdim = (Int_t)( (pM - pm) / phi ) ;
565 TH2F * histoRparticle = new TH2F("histoRparticle", histoname,
566 pdim, pm, pM, tdim, tm, tM) ;
567 histoRparticle->SetStats(kFALSE) ;
568 Text_t canvasname[80] ;
569 sprintf(canvasname, "Reconstructed particles in PHOSmodule # %d", module) ;
570 TCanvas * rparticlecanvas = new TCanvas("RparticleCanvas", canvasname, 650, 500) ;
571 RecParticlesList * rpl = fPHOS->RecParticles() ;
572 Int_t nRecParticles = rpl->GetEntries() ;
573 Int_t nRecParticlesInModule = 0 ;
574 TIter nextRecPart(rpl) ;
575 AliPHOSRecParticle * rp ;
576 cout << "DisplayRecParticles > " << nRecParticles << " reconstructed particles " << endl ;
577 Double_t kRADDEG = 180. / TMath::Pi() ;
578 while ( (rp = (AliPHOSRecParticle *)nextRecPart() ) ) {
579 AliPHOSTrackSegment * ts = rp->GetPHOSTrackSegment() ;
580 if ( ts->GetPHOSMod() == module ) {
581 Int_t numberofprimaries = 0 ;
582 Int_t * listofprimaries = rp->GetPrimaries(numberofprimaries) ;
583 cout << "Number of primaries = " << numberofprimaries << endl ;
585 for ( index = 0 ; index < numberofprimaries ; index++)
586 cout << " primary # " << index << " = " << listofprimaries[index] << endl ;
588 nRecParticlesInModule++ ;
589 Double_t theta = rp->Theta() * kRADDEG ;
590 Double_t phi = rp->Phi() * kRADDEG ;
591 Double_t energy = rp->Energy() ;
592 histoRparticle->Fill(phi, theta, energy) ;
595 histoRparticle->Draw("color") ;
597 nextRecPart.Reset() ;
598 while ( (rp = (AliPHOSRecParticle *)nextRecPart() ) ) {
599 AliPHOSTrackSegment * ts = rp->GetPHOSTrackSegment() ;
600 if ( ts->GetPHOSMod() == module )
605 sprintf(text, "reconstructed particles: %d", nRecParticlesInModule) ;
606 TPaveText * pavetext = new TPaveText(292, 100, 300, 101);
607 pavetext->AddText(text) ;
609 rparticlecanvas->Update() ;
613 //____________________________________________________________________________
614 void AliPHOSAnalyze::DisplayRecPoints()
616 // Display reconstructed points in local PHOS-module (x, z) coordinates.
617 // One PHOS module at the time.
618 // Click on symbols displays the EMC cluster, or PPSD information.
621 cout << "DisplayRecPoints > Analyze an event first ... (y/n) " ;
623 cin >> answer ; cout << answer ;
630 cout << "DisplayRecPoints > which module (1-5, -1: all) ? " ;
631 cin >> module ; cout << module << endl ;
633 Text_t canvasname[80];
634 sprintf(canvasname,"Digits in PHOS/EMC module # %d",module) ;
635 TCanvas * modulecanvas = new TCanvas("module", canvasname, 650, 500) ;
636 modulecanvas->Draw() ;
638 //=========== Creating 2d-histogram of the PHOS module
639 // a little bit junkie but is used to test Geom functinalities
641 Double_t tm, tM, pm, pM ; // min and Max theta and phi covered by module
643 fGeom->EmcModuleCoverage(module, tm, tM, pm, pM);
644 // convert angles into coordinates local to the EMC module of interest
646 Int_t emcModuleNumber ;
647 Double_t emcModulexm, emcModulezm ; // minimum local coordinate in a given EMCA module
648 Double_t emcModulexM, emcModulezM ; // maximum local coordinate in a given EMCA module
649 fGeom->ImpactOnEmc(tm, pm, emcModuleNumber, emcModulezm, emcModulexm) ;
650 fGeom->ImpactOnEmc(tM, pM, emcModuleNumber, emcModulezM, emcModulexM) ;
651 Int_t xdim = (Int_t)( ( emcModulexM - emcModulexm ) / fGeom->GetCrystalSize(0) ) ;
652 Int_t zdim = (Int_t)( ( emcModulezM - emcModulezm ) / fGeom->GetCrystalSize(2) ) ;
653 Float_t xmin = emcModulexm - fGeom->GetCrystalSize(0) ;
654 Float_t xMax = emcModulexM + fGeom->GetCrystalSize(0) ;
655 Float_t zmin = emcModulezm - fGeom->GetCrystalSize(2) ;
656 Float_t zMax = emcModulezM + fGeom->GetCrystalSize(2) ;
657 Text_t histoname[80];
658 sprintf(histoname,"Event %d: Digits and RecPoints in module %d", fEvt, module) ;
659 TH2F * hModule = new TH2F("HistoReconstructed", histoname,
660 xdim, xmin, xMax, zdim, zmin, zMax) ;
661 hModule->SetMaximum(2.0);
662 hModule->SetMinimum(0.0);
663 hModule->SetStats(kFALSE);
665 TIter next(fPHOS->Digits()) ;
666 Float_t energy, y, z;
668 Int_t relid[4]; Int_t nDigits = 0 ;
669 AliPHOSDigit * digit ;
671 // Making 2D histogram of the EMC module
672 while((digit = (AliPHOSDigit *)next()))
674 fGeom->AbsToRelNumbering(digit->GetId(), relid) ;
675 if (relid[0] == module && relid[1] == 0)
677 energy = fClu->Calibrate(digit->GetAmp()) ;
678 cout << "Energy is " << energy << " and threshold is " << fClu->GetEmcEnergyThreshold() << endl;
679 if (energy > fClu->GetEmcEnergyThreshold() ){
682 fGeom->RelPosInModule(relid,y,z) ;
683 hModule->Fill(y, z, energy) ;
687 cout <<"DrawRecPoints > Found in module "
688 << module << " " << nDigits << " digits with total energy " << etot << endl ;
689 hModule->Draw("col2") ;
691 //=========== Cluster in module
693 // TClonesArray * emcRP = fPHOS->EmcClusters() ;
694 TObjArray * emcRP = fPHOS->EmcRecPoints() ;
697 Int_t totalnClusters = 0 ;
698 Int_t nClusters = 0 ;
699 TIter nextemc(emcRP) ;
700 AliPHOSEmcRecPoint * emc ;
701 while((emc = (AliPHOSEmcRecPoint *)nextemc()))
703 // Int_t numberofprimaries ;
704 // Int_t * primariesarray = new Int_t[10] ;
705 // emc->GetPrimaries(numberofprimaries, primariesarray) ;
707 if ( emc->GetPHOSMod() == module )
710 energy = emc->GetTotalEnergy() ;
715 cout << "DrawRecPoints > Found " << totalnClusters << " EMC Clusters in PHOS" << endl ;
716 cout << "DrawRecPoints > Found in module " << module << " " << nClusters << " EMC Clusters " << endl ;
717 cout << "DrawRecPoints > total energy " << etot << endl ;
719 TPaveText * pavetext = new TPaveText(22, 80, 83, 90);
721 sprintf(text, "digits: %d; clusters: %d", nDigits, nClusters) ;
722 pavetext->AddText(text) ;
724 modulecanvas->Update();
726 //=========== Cluster in module PPSD Down
728 // TClonesArray * ppsdRP = fPHOS->PpsdClusters() ;
729 TObjArray * ppsdRP = fPHOS->PpsdRecPoints() ;
732 TIter nextPpsd(ppsdRP) ;
733 AliPHOSPpsdRecPoint * ppsd ;
734 while((ppsd = (AliPHOSPpsdRecPoint *)nextPpsd()))
737 if ( ppsd->GetPHOSMod() == module )
740 energy = ppsd->GetEnergy() ;
742 if (!ppsd->GetUp()) ppsd->Draw("P") ;
745 cout << "DrawRecPoints > Found " << totalnClusters << " Ppsd Down Clusters in PHOS" << endl ;
746 cout << "DrawRecPoints > Found in module " << module << " " << nClusters << " Ppsd Down Clusters " << endl ;
747 cout << "DrawRecPoints > total energy " << etot << endl ;
749 //=========== Cluster in module PPSD Up
751 ppsdRP = fPHOS->PpsdRecPoints() ;
754 TIter nextPpsdUp(ppsdRP) ;
755 while((ppsd = (AliPHOSPpsdRecPoint *)nextPpsdUp()))
758 if ( ppsd->GetPHOSMod() == module )
761 energy = ppsd->GetEnergy() ;
763 if (ppsd->GetUp()) ppsd->Draw("P") ;
766 cout << "DrawRecPoints > Found " << totalnClusters << " Ppsd Up Clusters in PHOS" << endl ;
767 cout << "DrawRecPoints > Found in module " << module << " " << nClusters << " Ppsd Up Clusters " << endl ;
768 cout << "DrawRecPoints > total energy " << etot << endl ;
773 //____________________________________________________________________________
774 void AliPHOSAnalyze::DisplayTrackSegments()
776 // Display track segments in local PHOS-module (x, z) coordinates.
777 // One PHOS module at the time.
778 // One symbol per PHOS subsystem: EMC, upper PPSD, lower PPSD.
781 cout << "DisplayTrackSegments > Analyze an event first ... (y/n) " ;
783 cin >> answer ; cout << answer ;
790 cout << "DisplayTrackSegments > which module (1-5, -1: all) ? " ;
791 cin >> module ; cout << module << endl ;
792 //=========== Creating 2d-histogram of the PHOS module
793 // a little bit junkie but is used to test Geom functinalities
795 Double_t tm, tM, pm, pM ; // min and Max theta and phi covered by module
797 fGeom->EmcModuleCoverage(module, tm, tM, pm, pM);
798 // convert angles into coordinates local to the EMC module of interest
800 Int_t emcModuleNumber ;
801 Double_t emcModulexm, emcModulezm ; // minimum local coordinate in a given EMCA module
802 Double_t emcModulexM, emcModulezM ; // maximum local coordinate in a given EMCA module
803 fGeom->ImpactOnEmc(tm, pm, emcModuleNumber, emcModulezm, emcModulexm) ;
804 fGeom->ImpactOnEmc(tM, pM, emcModuleNumber, emcModulezM, emcModulexM) ;
805 Int_t xdim = (Int_t)( ( emcModulexM - emcModulexm ) / fGeom->GetCrystalSize(0) ) ;
806 Int_t zdim = (Int_t)( ( emcModulezM - emcModulezm ) / fGeom->GetCrystalSize(2) ) ;
807 Float_t xmin = emcModulexm - fGeom->GetCrystalSize(0) ;
808 Float_t xMax = emcModulexM + fGeom->GetCrystalSize(0) ;
809 Float_t zmin = emcModulezm - fGeom->GetCrystalSize(2) ;
810 Float_t zMax = emcModulezM + fGeom->GetCrystalSize(2) ;
811 Text_t histoname[80];
812 sprintf(histoname,"Event %d: Track Segments in module %d", fEvt, module) ;
813 TH2F * histotrack = new TH2F("histotrack", histoname,
814 xdim, xmin, xMax, zdim, zmin, zMax) ;
815 histotrack->SetStats(kFALSE);
816 Text_t canvasname[80];
817 sprintf(canvasname,"Track segments in PHOS/EMC-PPSD module # %d", module) ;
818 TCanvas * trackcanvas = new TCanvas("TrackSegmentCanvas", canvasname, 650, 500) ;
821 TrackSegmentsList * trsegl = fPHOS->TrackSegments() ;
822 AliPHOSTrackSegment * trseg ;
824 Int_t nTrackSegments = trsegl->GetEntries() ;
827 Int_t nTrackSegmentsInModule = 0 ;
828 for(index = 0; index < nTrackSegments ; index++){
829 trseg = (AliPHOSTrackSegment * )trsegl->At(index) ;
830 etot+= trseg->GetEnergy() ;
831 if ( trseg->GetPHOSMod() == module ) {
832 nTrackSegmentsInModule++ ;
837 sprintf(text, "track segments: %d", nTrackSegmentsInModule) ;
838 TPaveText * pavetext = new TPaveText(22, 80, 83, 90);
839 pavetext->AddText(text) ;
841 trackcanvas->Update() ;
842 cout << "DisplayTrackSegments > Found " << trsegl->GetEntries() << " Track segments with total energy "<< etot << endl ;
846 //____________________________________________________________________________
847 Bool_t AliPHOSAnalyze::OpenRootFile(Text_t * name)
849 // Open the root file named "name"
851 fRootFile = new TFile(name, "update") ;
852 return fRootFile->IsOpen() ;
854 //____________________________________________________________________________
855 void AliPHOSAnalyze::SavingHistograms()
857 // Saves the histograms in a root file named "name.analyzed"
859 Text_t outputname[80] ;
860 sprintf(outputname,"%s.analyzed",fRootFile->GetName());
861 TFile output(outputname,"RECREATE");
864 fhEmcDigit->Write() ;
866 fhVetoDigit->Write() ;
867 if (fhConvertorDigit )
868 fhConvertorDigit->Write() ;
870 fhEmcCluster->Write() ;
872 fhVetoCluster->Write() ;
873 if (fhConvertorCluster )
874 fhConvertorCluster->Write() ;
876 fhConvertorEmc->Write() ;
878 fhPhotonEnergy->Write() ;
879 if (fhPhotonPositionX)
880 fhPhotonPositionX->Write() ;
881 if (fhPhotonPositionY)
882 fhPhotonPositionX->Write() ;
883 if (fhElectronEnergy)
884 fhElectronEnergy->Write() ;
885 if (fhElectronPositionX)
886 fhElectronPositionX->Write() ;
887 if (fhElectronPositionY)
888 fhElectronPositionX->Write() ;
889 if (fhNeutralHadronEnergy)
890 fhNeutralHadronEnergy->Write() ;
891 if (fhNeutralHadronPositionX)
892 fhNeutralHadronPositionX->Write() ;
893 if (fhNeutralHadronPositionY)
894 fhNeutralHadronPositionX->Write() ;
895 if (fhNeutralEMEnergy)
896 fhNeutralEMEnergy->Write() ;
897 if (fhNeutralEMPositionX)
898 fhNeutralEMPositionX->Write() ;
899 if (fhNeutralEMPositionY)
900 fhNeutralEMPositionX->Write() ;
901 if (fhChargedHadronEnergy)
902 fhChargedHadronEnergy->Write() ;
903 if (fhChargedHadronPositionX)
904 fhChargedHadronPositionX->Write() ;
905 if (fhChargedHadronPositionY)
906 fhChargedHadronPositionX->Write() ;
907 if (fhPhotonHadronEnergy)
908 fhPhotonHadronEnergy->Write() ;
909 if (fhPhotonHadronPositionX)
910 fhPhotonHadronPositionX->Write() ;
911 if (fhPhotonHadronPositionY)
912 fhPhotonHadronPositionX->Write() ;