]> git.uio.no Git - u/mrichter/AliRoot.git/blame - PHOS/AliPHOSAnalyze.cxx
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
Name-convetion rule violation corrected
[u/mrichter/AliRoot.git] / PHOS / AliPHOSAnalyze.cxx
CommitLineData
6ad0bfa0 1/**************************************************************************
2 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
3 * *
4 * Author: The ALICE Off-line Project. *
5 * Contributors are mentioned in the code where appropriate. *
6 * *
7 * Permission to use, copy, modify and distribute this software and its *
8 * documentation strictly for non-commercial purposes is hereby granted *
9 * without fee, provided that the above copyright notice appears in all *
10 * copies and that both the copyright notice and this permission notice *
11 * appear in the supporting documentation. The authors make no claims *
12 * about the suitability of this software for any purpose. It is *
13 * provided "as is" without express or implied warranty. *
14 **************************************************************************/
15
b2a60966 16/* $Id$ */
17
6ad0bfa0 18//_________________________________________________________________________
5f20d3fb 19// Algorythm class to analyze PHOSv1 events:
b2a60966 20// Construct histograms and displays them.
21// Use the macro EditorBar.C for best access to the functionnalities
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 {
eecb6765 77 //========== Get AliRun object from file
78 gAlice = (AliRun*) fRootFile->Get("gAlice") ;
fc879520 79
eecb6765 80 //=========== Get the PHOS object and associated geometry from the file
788dcca4 81 fPHOS = (AliPHOSv1 *)gAlice->GetDetector("PHOS") ;
eecb6765 82 fGeom = AliPHOSGeometry::GetInstance( fPHOS->GetGeometry()->GetName(), fPHOS->GetGeometry()->GetTitle() );
fc879520 83
eecb6765 84 //========== Initializes the Index to Object converter
85 fObjGetter = AliPHOSIndexToObject::GetInstance(fPHOS) ;
86 //========== Current event number
87 fEvt = -999 ;
fc879520 88
6ad0bfa0 89 }
69183710 90 fClu = 0 ;
91 fPID = 0 ;
92 fTrs = 0 ;
93 fRec = 0 ;
46b146ca 94 ResetHistograms() ;
6ad0bfa0 95}
96
88714635 97//____________________________________________________________________________
98AliPHOSAnalyze::AliPHOSAnalyze(const AliPHOSAnalyze & ana)
99{
100 // copy ctor
101 ( (AliPHOSAnalyze &)ana ).Copy(*this) ;
102}
103
104//____________________________________________________________________________
191c1c41 105void AliPHOSAnalyze::Copy(TObject & obj)
88714635 106{
107 // copy an analysis into an other one
108 TObject::Copy(obj) ;
109 // I do nothing more because the copy is silly but the Code checkers requires one
110}
111
6ad0bfa0 112//____________________________________________________________________________
113AliPHOSAnalyze::~AliPHOSAnalyze()
114{
115 // dtor
116
83974468 117 if (fRootFile->IsOpen() )
118 fRootFile->Close() ;
69183710 119 if(fRootFile)
120 delete fRootFile ;
121
122 if(fPHOS)
123 delete fPHOS ;
124
125 if(fClu)
126 delete fClu ;
127
128 if(fPID)
129 delete fPID ;
130
131 if(fRec)
132 delete fRec ;
133
134 if(fTrs)
135 delete fTrs ;
6ad0bfa0 136
137}
138
139//____________________________________________________________________________
46b146ca 140void AliPHOSAnalyze::ActivePPSD(Int_t Nevents=1){
6ad0bfa0 141
46b146ca 142 fhEnergyCorrelations = new TH2F("hEnergyCorrelations","hEnergyCorrelations",40, 0., 0.15, 30, 0., 3.e-5);
143 //========== Create the Clusterizer
144 fClu = new AliPHOSClusterizerv1() ;
69183710 145 fClu->SetEmcEnergyThreshold(0.01) ;
46b146ca 146 fClu->SetEmcClusteringThreshold(0.20) ;
147 fClu->SetPpsdEnergyThreshold (0.0000002) ;
148 fClu->SetPpsdClusteringThreshold(0.0000001) ;
69183710 149 fClu->SetLocalMaxCut(0.02) ;
46b146ca 150 fClu->SetCalibrationParameters(0., 0.00000001) ;
b2a60966 151
46b146ca 152 Int_t ievent;
83974468 153
46b146ca 154 for ( ievent=0; ievent<Nevents; ievent++)
155 {
156
157 //========== Event Number>
158 if ( ( log10((Float_t)(ievent+1)) - (Int_t)(log10((Float_t)(ievent+1))) ) == 0. )
159 cout << "AnalyzeResolutions > " << "Event is " << ievent << endl ;
160
161 //=========== Connects the various Tree's for evt
162 gAlice->GetEvent(ievent);
b2a60966 163
46b146ca 164 //=========== Gets the Kine TTree
165 gAlice->TreeK()->GetEvent(0) ;
69183710 166
46b146ca 167 //=========== Get the Digit Tree
168 gAlice->TreeD()->GetEvent(0) ;
169
170 //========== Creating branches ===================================
171 AliPHOSRecPoint::RecPointsList ** EmcRecPoints = fPHOS->EmcRecPoints() ;
172 gAlice->TreeR()->SetBranchAddress( "PHOSEmcRP", EmcRecPoints ) ;
173
174 AliPHOSRecPoint::RecPointsList ** PpsdRecPoints = fPHOS->PpsdRecPoints() ;
175 gAlice->TreeR()->SetBranchAddress( "PHOSPpsdRP", PpsdRecPoints ) ;
176
177 AliPHOSTrackSegment::TrackSegmentsList ** TrackSegmentsList = fPHOS->TrackSegments() ;
178 if( (*TrackSegmentsList) )
179 (*TrackSegmentsList)->Clear() ;
180 gAlice->TreeR()->SetBranchAddress( "PHOSTS", TrackSegmentsList ) ;
181
182 AliPHOSRecParticle::RecParticlesList ** RecParticleList = fPHOS->RecParticles() ;
183 if( (*RecParticleList) )
184 (*RecParticleList)->Clear() ;
185 gAlice->TreeR()->SetBranchAddress( "PHOSRP", RecParticleList ) ;
186
187
188 //=========== Gets the Reconstraction TTree
189 gAlice->TreeR()->GetEvent(0) ;
190
191 AliPHOSPpsdRecPoint * RecPoint ;
192 Int_t relid[4] ;
193 TIter nextRP(*fPHOS->PpsdRecPoints() ) ;
194 while( ( RecPoint = (AliPHOSPpsdRecPoint *)nextRP() ) )
195 {
196 if(!(RecPoint->GetUp()) ) {
197 AliPHOSDigit *digit ;
198 Int_t iDigit ;
199 for(iDigit = 0; iDigit < fPHOS->Digits()->GetEntries(); iDigit++)
200 {
201 digit = (AliPHOSDigit *) fPHOS->Digits()->At(iDigit) ;
202 fGeom->AbsToRelNumbering(digit->GetId(), relid) ;
203 if((relid[2]==1)&&(relid[3]==1)&&(relid[0]==RecPoint->GetPHOSMod())){
204 Float_t ConvertorEnergy = fClu->Calibrate(digit->GetAmp()) ;
205 fhEnergyCorrelations->Fill(ConvertorEnergy,RecPoint->GetTotalEnergy() );
206 break ;
207 }
208 }
209 break ;
210 }
211 }
212 }
213 SaveHistograms() ;
214 fhEnergyCorrelations->Draw("BOX") ;
6ad0bfa0 215}
216
69183710 217
92862013 218//____________________________________________________________________________
c3b9b3f9 219void AliPHOSAnalyze::AnalyzeManyEvents(Int_t Nevents, Int_t module)
92862013 220{
b2a60966 221 // analyzes Nevents events in a single PHOS module
46b146ca 222 // Events should be reconstructed by Reconstruct()
92862013 223
224 if ( fRootFile == 0 )
225 cout << "AnalyzeManyEvents > " << "Root File not openned" << endl ;
226 else
227 {
92862013 228 //========== Booking Histograms
229 cout << "AnalyzeManyEvents > " << "Booking Histograms" << endl ;
230 BookingHistograms();
46b146ca 231
92862013 232 Int_t ievent;
6a3f1304 233 Int_t relid[4] ;
92862013 234 AliPHOSDigit * digit ;
46b146ca 235 AliPHOSEmcRecPoint * emc ;
236 AliPHOSPpsdRecPoint * ppsd ;
38f8ae6d 237 // AliPHOSTrackSegment * tracksegment ;
238 AliPHOSRecParticle * recparticle;
46b146ca 239
92862013 240 for ( ievent=0; ievent<Nevents; ievent++)
241 {
de9ec31b 242 //========== Event Number>
46b146ca 243 if ( ( log10((Float_t)(ievent+1)) - (Int_t)(log10((Float_t)(ievent+1))) ) == 0. )
244 cout << "AnalyzeManyEvents > " << "Event is " << ievent << endl ;
245
92862013 246 //=========== Connects the various Tree's for evt
247 gAlice->GetEvent(ievent);
46b146ca 248
92862013 249 //=========== Gets the Digit TTree
250 gAlice->TreeD()->GetEvent(0) ;
46b146ca 251
92862013 252 //=========== Gets the number of entries in the Digits array
253 TIter nextdigit(fPHOS->Digits()) ;
254 while( ( digit = (AliPHOSDigit *)nextdigit() ) )
255 {
256 fGeom->AbsToRelNumbering(digit->GetId(), relid) ;
46b146ca 257 if (fClu->IsInEmc(digit)) fhEmcDigit->Fill(fClu->Calibrate(digit->GetAmp())) ;
258 else
c3b9b3f9 259 {
260 if (relid[1]<17) fhVetoDigit->Fill(fClu->Calibrate(digit->GetAmp()));
261 if (relid[1]>16) fhConvertorDigit->Fill(fClu->Calibrate(digit->GetAmp()));
262 }
92862013 263 }
11e586af 264
46b146ca 265
266 //=========== Cluster in module
267 TIter nextEmc(*fPHOS->EmcRecPoints() ) ;
268 while((emc = (AliPHOSEmcRecPoint *)nextEmc()))
269 {
270 if ( emc->GetPHOSMod() == module )
271 {
272 fhEmcCluster->Fill( emc->GetTotalEnergy() );
273 TIter nextPpsd( *fPHOS->PpsdRecPoints()) ;
274 while((ppsd = (AliPHOSPpsdRecPoint *)nextPpsd()))
275 {
276 if ( ppsd->GetPHOSMod() == module )
277 {
278 if (!ppsd->GetUp()) fhConvertorEmc->Fill(emc->GetTotalEnergy(),ppsd->GetTotalEnergy()) ;
279 }
280 }
281 }
282 }
283
284 //=========== Cluster in module PPSD Down
285 TIter nextPpsd(*fPHOS->PpsdRecPoints() ) ;
286 while((ppsd = (AliPHOSPpsdRecPoint *)nextPpsd()))
287 {
288 if ( ppsd->GetPHOSMod() == module )
289 {
290 if (!ppsd->GetUp()) fhConvertorCluster->Fill(ppsd->GetTotalEnergy()) ;
291 if (ppsd->GetUp()) fhVetoCluster ->Fill(ppsd->GetTotalEnergy()) ;
292 }
293 }
294
92862013 295 //========== TRackSegments in the event
fc879520 296 TIter nextRecParticle(*fPHOS->RecParticles() ) ;
38f8ae6d 297 while((recparticle = (AliPHOSRecParticle *)nextRecParticle()))
92862013 298 {
38f8ae6d 299 if ( recparticle->GetPHOSTrackSegment()->GetPHOSMod() == module )
92862013 300 {
b2a60966 301 cout << "Particle type is " << recparticle->GetType() << endl ;
302 Int_t numberofprimaries = 0 ;
303 Int_t * listofprimaries = recparticle->GetPrimaries(numberofprimaries) ;
304 cout << "Number of primaries = " << numberofprimaries << endl ;
305 Int_t index ;
306 for ( index = 0 ; index < numberofprimaries ; index++)
307 cout << " primary # " << index << " = " << listofprimaries[index] << endl ;
134ce69a 308 }
309 }
134ce69a 310 } // endfor
46b146ca 311 SaveHistograms();
134ce69a 312 } // endif
313} // endfunction
314
315//____________________________________________________________________________
69183710 316 void AliPHOSAnalyze::Reconstruct(Int_t Nevents,Int_t FirstEvent )
c3b9b3f9 317{
c3b9b3f9 318 Int_t ievent ;
69183710 319 for ( ievent=FirstEvent; ievent<Nevents; ievent++)
fc879520 320 {
69183710 321 if (ievent==FirstEvent)
fc879520 322 {
46b146ca 323 cout << "Analyze > Starting Reconstructing " << endl ;
fc879520 324 //========== Create the Clusterizer
325 fClu = new AliPHOSClusterizerv1() ;
69183710 326 fClu->SetEmcEnergyThreshold(0.03) ;
fc879520 327 fClu->SetEmcClusteringThreshold(0.20) ;
69183710 328 fClu->SetPpsdEnergyThreshold (0.0000001) ;
329 fClu->SetPpsdClusteringThreshold(0.0000001) ;
8f5ada7b 330 fClu->SetLocalMaxCut(0.03) ;
fc879520 331 fClu->SetCalibrationParameters(0., 0.00000001) ;
332
333 //========== Creates the track segment maker
334 fTrs = new AliPHOSTrackSegmentMakerv1() ;
335 // fTrs->UnsetUnfoldFlag() ;
336
337 //========== Creates the particle identifier
338 fPID = new AliPHOSPIDv1() ;
339 fPID->SetShowerProfileCuts(0.3, 1.8, 0.3, 1.8 ) ;
69183710 340 fPID->SetDispersionCutOff(2.0) ;
341 fPID->SetRelativeDistanceCut(3.) ;
fc879520 342
343 //========== Creates the Reconstructioner
344 fRec = new AliPHOSReconstructioner(fClu, fTrs, fPID) ;
69183710 345 // fRec -> SetDebugReconstruction(kTRUE);
134ce69a 346
fc879520 347 }
348
349 //========== Event Number>
69183710 350 // if ( ( log10((Float_t)(ievent+1)) - (Int_t)(log10((Float_t)(ievent+1))) ) == 0. )
351 cout << "Reconstruct > Event is " << ievent << endl ;
fc879520 352
353 //=========== Connects the various Tree's for evt
354 gAlice->GetEvent(ievent);
134ce69a 355
fc879520 356 //=========== Gets the Digit TTree
357 gAlice->TreeD()->GetEvent(0) ;
134ce69a 358
fc879520 359 //=========== Do the reconstruction
360 fPHOS->Reconstruction(fRec);
361 }
134ce69a 362
69183710 363 fClu->Delete();
364 fClu=0 ;
365 fTrs->Delete();
366 fTrs = 0 ;
367 fPID->Delete();
368 fPID = 0 ;
369 fRec->Delete();
370 fRec = 0 ;
c3b9b3f9 371
c3b9b3f9 372}
c3b9b3f9 373//____________________________________________________________________________
69183710 374 void AliPHOSAnalyze::InvariantMass(Int_t Nevents )
c3b9b3f9 375{
69183710 376 // Calculates Real and Mixed invariant mass distributions
377 Int_t NMixedEvents = 4 ; //# of events used for calculation of 'mixed' distribution
378 Int_t MixedLoops = (Int_t )TMath::Ceil(Nevents/NMixedEvents) ;
379
380 //========== Booking Histograms
381 TH2D * hRealEM = new TH2D("hRealEM", "Real for EM particles", 250,0.,1.,40,0.,4.) ;
382 TH2D * hRealPhot = new TH2D("hRealPhot", "Real for kPhoton particles", 250,0.,1.,40,0.,4.) ;
383 TH2D * hMixedEM = new TH2D("hMixedEM", "Mixed for EM particles", 250,0.,1.,40,0.,4.) ;
384 TH2D * hMixedPhot= new TH2D("hMixedPhot","Mixed for kPhoton particles",250,0.,1.,40,0.,4.) ;
385
386 Int_t ievent;
387 Int_t EventInMixedLoop ;
388
389 Int_t NRecParticles[NMixedEvents] ;
390
391 AliPHOSRecParticle::RecParticlesList * AllRecParticleList = new TClonesArray("AliPHOSRecParticle", NMixedEvents*1000) ;
392
393 for(EventInMixedLoop = 0; EventInMixedLoop < MixedLoops; EventInMixedLoop++ ){
394 Int_t iRecPhot = 0 ;
c3b9b3f9 395
69183710 396 for ( ievent=0; ievent < NMixedEvents; ievent++){
397
398 Int_t AbsEventNumber = EventInMixedLoop*NMixedEvents + ievent ;
399
400 //=========== Connects the various Tree's for evt
401 gAlice->GetEvent(AbsEventNumber);
402
403 //=========== Get the Digit Tree
404 gAlice->TreeD()->GetEvent(0) ;
405
406 //========== Creating branches ===================================
407
408 AliPHOSRecParticle::RecParticlesList ** RecParticleList = fPHOS->RecParticles() ;
409 if( (*RecParticleList) )
410 (*RecParticleList)->Clear() ;
411 gAlice->TreeR()->SetBranchAddress( "PHOSRP", RecParticleList ) ;
412
413 //=========== Gets the Reconstraction TTree
414 gAlice->TreeR()->GetEvent(0) ;
415
416 AliPHOSRecParticle * RecParticle ;
417 Int_t iRecParticle ;
418 for(iRecParticle = 0; iRecParticle < (*RecParticleList)->GetEntries() ;iRecParticle++ )
419 {
420 RecParticle = (AliPHOSRecParticle *) (*RecParticleList)->At(iRecParticle) ;
421 if((RecParticle->GetType() == AliPHOSFastRecParticle::kGAMMA)||
422 (RecParticle->GetType() == AliPHOSFastRecParticle::kNEUTRALEM)){
423 new( (*AllRecParticleList)[iRecPhot] ) AliPHOSRecParticle(*RecParticle) ;
424 iRecPhot++;
425 }
426 }
427
428 NRecParticles[ievent] = iRecPhot-1 ;
c3b9b3f9 429 }
69183710 430
c3b9b3f9 431
69183710 432 //Now calculate invariant mass:
433 Int_t irp1,irp2 ;
434 Int_t NCurEvent = 0 ;
c3b9b3f9 435
69183710 436 for(irp1 = 0; irp1 < AllRecParticleList->GetEntries()-1; irp1++){
437 AliPHOSRecParticle * rp1 = (AliPHOSRecParticle *)AllRecParticleList->At(irp1) ;
c3b9b3f9 438
69183710 439 for(irp2 = irp1+1; irp2 < AllRecParticleList->GetEntries(); irp2++){
440 AliPHOSRecParticle * rp2 = (AliPHOSRecParticle *)AllRecParticleList->At(irp2) ;
441
442 Double_t InvMass ;
443 InvMass = (rp1->Energy()+rp2->Energy())*(rp1->Energy()+rp2->Energy())-
444 (rp1->Px()+rp2->Px())*(rp1->Px()+rp2->Px())-
445 (rp1->Py()+rp2->Py())*(rp1->Py()+rp2->Py())-
446 (rp1->Pz()+rp2->Pz())*(rp1->Pz()+rp2->Pz()) ;
447
448 if(InvMass> 0)
449 InvMass = TMath::Sqrt(InvMass);
450
451 Double_t Pt ;
452 Pt = TMath::Sqrt((rp1->Px()+rp2->Px() )*( rp1->Px()+rp2->Px() ) +(rp1->Py()+rp2->Py())*(rp1->Py()+rp2->Py()));
453
454 if(irp1 > NRecParticles[NCurEvent])
455 NCurEvent++;
456
457 if(irp2 <= NRecParticles[NCurEvent]){ //'Real' event
458 hRealEM->Fill(InvMass,Pt);
459 if((rp1->GetType() == AliPHOSFastRecParticle::kGAMMA)&&(rp2->GetType() == AliPHOSFastRecParticle::kGAMMA))
460 hRealPhot->Fill(InvMass,Pt);
461 }
462 else{
463 hMixedEM->Fill(InvMass,Pt);
464 if((rp1->GetType() == AliPHOSFastRecParticle::kGAMMA)&&(rp2->GetType() == AliPHOSFastRecParticle::kGAMMA))
465 hMixedPhot->Fill(InvMass,Pt);
466 } //real-mixed
467
468 } //loop over second rp
469 }//loop over first rp
c3b9b3f9 470
c3b9b3f9 471
69183710 472 AllRecParticleList->Delete() ;
473 } //Loop over events
474
475 delete AllRecParticleList ;
476
477 //writing output
478 TFile output("invmass.root","RECREATE");
479 output.cd();
480
481 hRealEM->Write() ;
482 hRealPhot->Write() ;
483 hMixedEM->Write() ;
484 hMixedPhot->Write() ;
485
486 output.Write();
487 output.Close();
c3b9b3f9 488
489}
490
fc879520 491//____________________________________________________________________________
492 void AliPHOSAnalyze::AnalyzeResolutions(Int_t Nevents )
493{
494 // analyzes Nevents events and calculate Energy and Position resolution as well as
495 // probaility of correct indentifiing of the incident particle
134ce69a 496
fc879520 497 //========== Booking Histograms
498 cout << "AnalyzeResolutions > " << "Booking Histograms" << endl ;
499 BookResolutionHistograms();
134ce69a 500
fc879520 501 Int_t Counter[9][5] ;
502 Int_t i1,i2,TotalInd = 0 ;
503 for(i1 = 0; i1<9; i1++)
504 for(i2 = 0; i2<5; i2++)
505 Counter[i1][i2] = 0 ;
506
507 Int_t TotalPrimary = 0 ;
508 Int_t TotalRecPart = 0 ;
509 Int_t TotalRPwithPrim = 0 ;
510 Int_t ievent;
511
512 cout << "Start Analysing"<< endl ;
513 for ( ievent=0; ievent<Nevents; ievent++)
514 {
515
516 //========== Event Number>
69183710 517 // if ( ( log10((Float_t)(ievent+1)) - (Int_t)(log10((Float_t)(ievent+1))) ) == 0. )
fc879520 518 cout << "AnalyzeResolutions > " << "Event is " << ievent << endl ;
519
520 //=========== Connects the various Tree's for evt
521 gAlice->GetEvent(ievent);
134ce69a 522
69183710 523 //=========== Gets the Kine TTree
fc879520 524 gAlice->TreeK()->GetEvent(0) ;
525
526 //=========== Gets the list of Primari Particles
527 TClonesArray * PrimaryList = gAlice->Particles();
528
529 TParticle * Primary ;
530 Int_t iPrimary ;
531 for ( iPrimary = 0 ; iPrimary < PrimaryList->GetEntries() ; iPrimary++)
69183710 532 {
533 Primary = (TParticle*)PrimaryList->UncheckedAt(iPrimary) ;
534 Int_t PrimaryType = Primary->GetPdgCode() ;
535 if( PrimaryType == 22 ) {
536 Int_t ModuleNumber ;
537 Double_t PrimX, PrimZ ;
538 fGeom->ImpactOnEmc(Primary->Theta(), Primary->Phi(), ModuleNumber, PrimX, PrimZ) ;
539 if(ModuleNumber){
540 fhPrimary->Fill(Primary->Energy()) ;
541 if(Primary->Energy() > 0.3)
542 TotalPrimary++ ;
543 }
544 }
545 }
fc879520 546
547 //=========== Get the Digit Tree
548 gAlice->TreeD()->GetEvent(0) ;
69183710 549
fc879520 550 //========== Creating branches ===================================
551 AliPHOSRecPoint::RecPointsList ** EmcRecPoints = fPHOS->EmcRecPoints() ;
552 gAlice->TreeR()->SetBranchAddress( "PHOSEmcRP", EmcRecPoints ) ;
69183710 553
fc879520 554 AliPHOSRecPoint::RecPointsList ** PpsdRecPoints = fPHOS->PpsdRecPoints() ;
555 gAlice->TreeR()->SetBranchAddress( "PHOSPpsdRP", PpsdRecPoints ) ;
69183710 556
fc879520 557 AliPHOSTrackSegment::TrackSegmentsList ** TrackSegmentsList = fPHOS->TrackSegments() ;
558 if( (*TrackSegmentsList) )
559 (*TrackSegmentsList)->Clear() ;
560 gAlice->TreeR()->SetBranchAddress( "PHOSTS", TrackSegmentsList ) ;
561
562 AliPHOSRecParticle::RecParticlesList ** RecParticleList = fPHOS->RecParticles() ;
563 if( (*RecParticleList) )
564 (*RecParticleList)->Clear() ;
565 gAlice->TreeR()->SetBranchAddress( "PHOSRP", RecParticleList ) ;
69183710 566
fc879520 567 //=========== Gets the Reconstraction TTree
568 gAlice->TreeR()->GetEvent(0) ;
69183710 569
fc879520 570 AliPHOSRecParticle * RecParticle ;
571 Int_t iRecParticle ;
572 for(iRecParticle = 0; iRecParticle < (*RecParticleList)->GetEntries() ;iRecParticle++ )
573 {
574 RecParticle = (AliPHOSRecParticle *) (*RecParticleList)->At(iRecParticle) ;
69183710 575 fhAllRP->Fill(CorrectEnergy(RecParticle->Energy())) ;
fc879520 576
577 Int_t ModuleNumberRec ;
578 Double_t RecX, RecZ ;
579 fGeom->ImpactOnEmc(RecParticle->Theta(), RecParticle->Phi(), ModuleNumberRec, RecX, RecZ) ;
580
8f5ada7b 581 Double_t MinDistance = 2. ;
fc879520 582 Int_t ClosestPrimary = -1 ;
583
584 Int_t numberofprimaries ;
585 Int_t * listofprimaries = RecParticle->GetPrimaries(numberofprimaries) ;
586 Int_t index ;
587 TParticle * Primary ;
588 Double_t Distance = MinDistance ;
69183710 589 for ( index = 0 ; index < numberofprimaries ; index++){
590 Primary = (TParticle*)PrimaryList->UncheckedAt(listofprimaries[index]) ;
591 Int_t ModuleNumber ;
592 Double_t PrimX, PrimZ ;
593 fGeom->ImpactOnEmc(Primary->Theta(), Primary->Phi(), ModuleNumber, PrimX, PrimZ) ;
594 if(ModuleNumberRec == ModuleNumber)
595 Distance = TMath::Sqrt((RecX-PrimX)*(RecX-PrimX)+(RecZ-PrimZ)*(RecZ-PrimZ) ) ;
596 if(MinDistance > Distance)
597 {
598 MinDistance = Distance ;
599 ClosestPrimary = listofprimaries[index] ;
600 }
601 }
fc879520 602 TotalRecPart++ ;
134ce69a 603
69183710 604 if(ClosestPrimary >=0 ){
605 TotalRPwithPrim++;
606
607 Int_t PrimaryType = ((TParticle *)PrimaryList->At(ClosestPrimary))->GetPdgCode() ;
608 TParticlePDG* PDGparticle = ((TParticle *)PrimaryList->At(ClosestPrimary))->GetPDG();
609 Double_t charge = PDGparticle->Charge() ;
8f5ada7b 610// if(charge)
611// cout <<"Primary " <<PrimaryType << " E " << ((TParticle *)PrimaryList->At(ClosestPrimary))->Energy() << endl ;
69183710 612 Int_t PrimaryCode ;
613 switch(PrimaryType)
614 {
615 case 22:
616 PrimaryCode = 0; //Photon
617 fhAllEnergy->Fill(((TParticle *) PrimaryList->At(ClosestPrimary))->Energy(), RecParticle->Energy()) ;
618 fhAllPosition->Fill(((TParticle *) PrimaryList->At(ClosestPrimary))->Energy(),MinDistance) ;
619 break;
620 case 11 :
621 PrimaryCode = 1; //Electron
622 break;
623 case -11 :
624 PrimaryCode = 1; //positron
625 break;
626 case 321 :
627 PrimaryCode = 4; //K+
628 break;
629 case -321 :
630 PrimaryCode = 4; //K-
631 break;
632 case 310 :
633 PrimaryCode = 4; //K0s
634 break;
635 case 130 :
636 PrimaryCode = 4; //K0l
637 break;
8f5ada7b 638 case 211 :
639 PrimaryCode = 2; //K0l
640 break;
641 case -211 :
642 PrimaryCode = 2; //K0l
643 break;
644 case 2212 :
645 PrimaryCode = 2; //K0l
646 break;
647 case -2212 :
648 PrimaryCode = 2; //K0l
649 break;
69183710 650 default:
8f5ada7b 651 PrimaryCode = 3; //ELSE
69183710 652 break;
653 }
654
655 switch(RecParticle->GetType())
656 {
657 case AliPHOSFastRecParticle::kGAMMA:
658 if(PrimaryType == 22){
659 fhPhotEnergy->Fill(((TParticle *) PrimaryList->At(ClosestPrimary))->Energy(), RecParticle->Energy() ) ;
660 fhEMEnergy->Fill(((TParticle *) PrimaryList->At(ClosestPrimary))->Energy(), RecParticle->Energy() ) ;
661 fhPPSDEnergy->Fill(((TParticle *) PrimaryList->At(ClosestPrimary))->Energy(), RecParticle->Energy() ) ;
662
663 fhPhotPosition->Fill(((TParticle *) PrimaryList->At(ClosestPrimary))->Energy(),MinDistance) ;
664 fhEMPosition->Fill(((TParticle *) PrimaryList->At(ClosestPrimary))->Energy(),MinDistance) ;
665 fhPPSDPosition->Fill(((TParticle *) PrimaryList->At(ClosestPrimary))->Energy(),MinDistance) ;
666
667 fhPhotReg->Fill(CorrectEnergy(RecParticle->Energy()) ) ;
668 fhPhotEM->Fill(CorrectEnergy(RecParticle->Energy()) ) ;
669 fhPhotPPSD->Fill(CorrectEnergy(RecParticle->Energy()) ) ;
670
671 fhPhotPhot->Fill(CorrectEnergy(RecParticle->Energy()) ) ;
672 }
673 if(PrimaryType == 2112){ //neutron
674 fhNReg->Fill(CorrectEnergy(RecParticle->Energy()) ) ;
675 fhNEM->Fill(CorrectEnergy(RecParticle->Energy()) ) ;
676 fhNPPSD->Fill(CorrectEnergy(RecParticle->Energy()) ) ;
677 }
678
679 if(PrimaryType == -2112){ //neutron ~
680 fhNBarReg->Fill(CorrectEnergy(RecParticle->Energy()) ) ;
681 fhNBarEM->Fill(CorrectEnergy(RecParticle->Energy()) ) ;
682 fhNBarPPSD->Fill(CorrectEnergy(RecParticle->Energy()) ) ;
683
684 }
685 if(PrimaryCode == 2){
686 fhChargedReg->Fill(CorrectEnergy(RecParticle->Energy()) ) ;
687 fhChargedEM->Fill(CorrectEnergy(RecParticle->Energy()) ) ;
688 fhChargedPPSD->Fill(CorrectEnergy(RecParticle->Energy()) ) ;
689 }
690
691 fhAllReg->Fill(CorrectEnergy(RecParticle->Energy()) ) ;
692 fhAllEM->Fill(CorrectEnergy(RecParticle->Energy()) ) ;
693 fhAllPPSD->Fill(CorrectEnergy(RecParticle->Energy()) ) ;
694 fhShape->Fill(CorrectEnergy(RecParticle->Energy()) ) ;
695 fhVeto->Fill(CorrectEnergy(RecParticle->Energy()) ) ;
696 fhPPSD->Fill(CorrectEnergy(RecParticle->Energy()) ) ;
697 Counter[0][PrimaryCode]++;
698 break;
699 case AliPHOSFastRecParticle::kELECTRON:
700 if(PrimaryType == 22){
701 fhPhotElec->Fill(CorrectEnergy(RecParticle->Energy()) ) ;
702 fhEMEnergy->Fill(((TParticle *) PrimaryList->At(ClosestPrimary))->Energy(), RecParticle->Energy() ) ;
703 fhEMPosition->Fill(((TParticle *) PrimaryList->At(ClosestPrimary))->Energy(),MinDistance) ;
704 fhPhotEM->Fill(CorrectEnergy(RecParticle->Energy()) ) ;
705 fhPhotPPSD->Fill(CorrectEnergy(RecParticle->Energy()) ) ;
706 }
707 if(PrimaryType == 2112){ //neutron
708 fhNEM->Fill(CorrectEnergy(RecParticle->Energy()) ) ;
709 fhNPPSD->Fill(CorrectEnergy(RecParticle->Energy()) ) ;
710 }
711
712 if(PrimaryType == -2112){ //neutron ~
713 fhNBarEM->Fill(CorrectEnergy(RecParticle->Energy()) ) ;
714 fhNBarPPSD->Fill(CorrectEnergy(RecParticle->Energy()) ) ;
715
716 }
717 if(PrimaryCode == 2){
718 fhChargedEM->Fill(CorrectEnergy(RecParticle->Energy()) ) ;
719 fhChargedPPSD->Fill(CorrectEnergy(RecParticle->Energy()) ) ;
720 }
721
722 fhAllEM->Fill(CorrectEnergy(RecParticle->Energy()) ) ;
723 fhAllPPSD->Fill(CorrectEnergy(RecParticle->Energy()) ) ;
724 fhShape->Fill(CorrectEnergy(RecParticle->Energy()) ) ;
725 fhPPSD->Fill(CorrectEnergy(RecParticle->Energy()) ) ;
726 Counter[1][PrimaryCode]++;
727 break;
728 case AliPHOSFastRecParticle::kNEUTRALHA:
729 if(PrimaryType == 22)
730 fhPhotNeuH->Fill(CorrectEnergy(RecParticle->Energy()) ) ;
731
732 fhVeto->Fill(CorrectEnergy(RecParticle->Energy()) ) ;
733 Counter[2][PrimaryCode]++;
734 break ;
735 case AliPHOSFastRecParticle::kNEUTRALEM:
736 if(PrimaryType == 22){
737 fhEMEnergy->Fill(((TParticle *)PrimaryList->At(ClosestPrimary))->Energy(),RecParticle->Energy() ) ;
738 fhEMPosition->Fill(((TParticle *)PrimaryList->At(ClosestPrimary))->Energy(),MinDistance ) ;
739
740 fhPhotNuEM->Fill(CorrectEnergy(RecParticle->Energy()) ) ;
741 fhPhotEM->Fill(CorrectEnergy(RecParticle->Energy()) ) ;
742 }
743 if(PrimaryType == 2112) //neutron
744 fhNEM->Fill(CorrectEnergy(RecParticle->Energy()) ) ;
745
746 if(PrimaryType == -2112) //neutron ~
747 fhNBarEM->Fill(CorrectEnergy(RecParticle->Energy()) ) ;
748
749 if(PrimaryCode == 2)
750 fhChargedEM->Fill(CorrectEnergy(RecParticle->Energy()) ) ;
751
752 fhAllEM->Fill(CorrectEnergy(RecParticle->Energy()) ) ;
753 fhShape->Fill(CorrectEnergy(RecParticle->Energy()) ) ;
754 fhVeto->Fill(CorrectEnergy(RecParticle->Energy()) ) ;
755
756 Counter[3][PrimaryCode]++;
757 break ;
758 case AliPHOSFastRecParticle::kCHARGEDHA:
759 if(PrimaryType == 22) //photon
760 fhPhotChHa->Fill(CorrectEnergy(RecParticle->Energy()) ) ;
761
762 Counter[4][PrimaryCode]++ ;
763 break ;
764 case AliPHOSFastRecParticle::kGAMMAHA:
765 if(PrimaryType == 22){ //photon
766 fhPhotGaHa->Fill(CorrectEnergy(RecParticle->Energy()) ) ;
767 fhPPSDEnergy->Fill(((TParticle *) PrimaryList->At(ClosestPrimary))->Energy(), RecParticle->Energy() ) ;
768 fhPPSDPosition->Fill(((TParticle *) PrimaryList->At(ClosestPrimary))->Energy(),MinDistance) ;
769 fhPhotPPSD->Fill(CorrectEnergy(RecParticle->Energy()) ) ;
fc879520 770 }
771 if(PrimaryType == 2112){ //neutron
69183710 772 fhNPPSD->Fill(CorrectEnergy(RecParticle->Energy()) ) ;
fc879520 773 }
69183710 774
fc879520 775 if(PrimaryType == -2112){ //neutron ~
69183710 776 fhNBarPPSD->Fill(CorrectEnergy(RecParticle->Energy()) ) ;
fc879520 777 }
778 if(PrimaryCode == 2){
69183710 779 fhChargedPPSD->Fill(CorrectEnergy(RecParticle->Energy()) ) ;
fc879520 780 }
69183710 781
782 fhAllPPSD->Fill(CorrectEnergy(RecParticle->Energy()) ) ;
783 fhVeto->Fill(CorrectEnergy(RecParticle->Energy()) ) ;
784 fhPPSD->Fill(CorrectEnergy(RecParticle->Energy()) ) ;
fc879520 785
fc879520 786 Counter[5][PrimaryCode]++ ;
787 break ;
69183710 788 case AliPHOSFastRecParticle::kABSURDEM:
789 Counter[6][PrimaryCode]++ ;
790 fhShape->Fill(CorrectEnergy(RecParticle->Energy()) ) ;
791 break;
792 case AliPHOSFastRecParticle::kABSURDHA:
793 Counter[7][PrimaryCode]++ ;
794 break;
795 default:
796 Counter[8][PrimaryCode]++ ;
797 break;
798 }
799 }
fc879520 800 }
801 } // endfor
46b146ca 802 SaveHistograms();
fc879520 803 cout << "Resolutions: Analyzed " << Nevents << " event(s)" << endl ;
804 cout << "Resolutions: Total primary " << TotalPrimary << endl ;
805 cout << "Resoluitons: Total reconstracted " << TotalRecPart << endl ;
806 cout << "TotalReconstructed with Primarie " << TotalRPwithPrim << endl ;
807 cout << " Primary: Photon Electron Ch. Hadr. Neutr. Hadr Kaons" << endl ;
808 cout << " Detected as photon " << Counter[0][0] << " " << Counter[0][1] << " " << Counter[0][2] << " " <<Counter[0][3] << " " << Counter[0][4] << endl ;
809 cout << " Detected as electron " << Counter[1][0] << " " << Counter[1][1] << " " << Counter[1][2] << " " <<Counter[1][3] << " " << Counter[1][4] << endl ;
810 cout << " Detected as neutral hadron " << Counter[2][0] << " " << Counter[2][1] << " " << Counter[2][2] << " " <<Counter[2][3] << " " << Counter[2][4] << endl ;
811 cout << " Detected as neutral EM " << Counter[3][0] << " " << Counter[3][1] << " " << Counter[3][2] << " " <<Counter[3][3] << " " << Counter[3][4] << endl ;
812 cout << " Detected as charged hadron " << Counter[4][0] << " " << Counter[4][1] << " " << Counter[4][2] << " " <<Counter[4][3] << " " << Counter[4][4] << endl ;
813 cout << " Detected as gamma-hadron " << Counter[5][0] << " " << Counter[5][1] << " " << Counter[5][2] << " " <<Counter[5][3] << " " << Counter[5][4] << endl ;
814 cout << " Detected as Absurd EM " << Counter[6][0] << " " << Counter[6][1] << " " << Counter[6][2] << " " <<Counter[6][3] << " " << Counter[6][4] << endl ;
815 cout << " Detected as absurd hadron " << Counter[7][0] << " " << Counter[7][1] << " " << Counter[7][2] << " " <<Counter[7][3] << " " << Counter[7][4] << endl ;
816 cout << " Detected as undefined " << Counter[8][0] << " " << Counter[8][1] << " " << Counter[8][2] << " " <<Counter[8][3] << " " << Counter[8][4] << endl ;
817
818 for(i1 = 0; i1<9; i1++)
819 for(i2 = 0; i2<5; i2++)
820 TotalInd+=Counter[i1][i2] ;
821 cout << "Indentified particles " << TotalInd << endl ;
822
92862013 823} // endfunction
824
825
826//____________________________________________________________________________
827void AliPHOSAnalyze::BookingHistograms()
828{
b2a60966 829 // Books the histograms where the results of the analysis are stored (to be changed)
830
eecb6765 831 delete fhEmcDigit ;
832 delete fhVetoDigit ;
833 delete fhConvertorDigit ;
834 delete fhEmcCluster ;
835 delete fhVetoCluster ;
836 delete fhConvertorCluster ;
837 delete fhConvertorEmc ;
838
46b146ca 839 fhEmcDigit = new TH1F("hEmcDigit", "hEmcDigit", 1000, 0. , 25.);
840 fhVetoDigit = new TH1F("hVetoDigit", "hVetoDigit", 500, 0. , 3.e-5);
841 fhConvertorDigit = new TH1F("hConvertorDigit","hConvertorDigit", 500, 0. , 3.e-5);
842 fhEmcCluster = new TH1F("hEmcCluster", "hEmcCluster", 1000, 0. , 30.);
843 fhVetoCluster = new TH1F("hVetoCluster", "hVetoCluster", 500, 0. , 3.e-5);
844 fhConvertorCluster = new TH1F("hConvertorCluster","hConvertorCluster",500, 0. , 3.e-5);
845 fhConvertorEmc = new TH2F("hConvertorEmc", "hConvertorEmc", 200, 1. , 3., 200, 0., 3.e-5);
92862013 846
134ce69a 847}
848//____________________________________________________________________________
849void AliPHOSAnalyze::BookResolutionHistograms()
850{
851 // Books the histograms where the results of the Resolution analysis are stored
852
69183710 853// if(fhAllEnergy)
854// delete fhAllEnergy ;
855// if(fhPhotEnergy)
856// delete fhPhotEnergy ;
857// if(fhEMEnergy)
858// delete fhEMEnergy ;
859// if(fhPPSDEnergy)
860// delete fhPPSDEnergy ;
861
862
863 fhAllEnergy = new TH2F("hAllEnergy", "Energy of any RP with primary photon",100, 0., 5., 100, 0., 5.);
864 fhPhotEnergy = new TH2F("hPhotEnergy", "Energy of kGAMMA with primary photon",100, 0., 5., 100, 0., 5.);
865 fhEMEnergy = new TH2F("hEMEnergy", "Energy of EM with primary photon", 100, 0., 5., 100, 0., 5.);
866 fhPPSDEnergy = new TH2F("hPPSDEnergy", "Energy of PPSD with primary photon", 100, 0., 5., 100, 0., 5.);
867
868// if(fhAllPosition)
869// delete fhAllPosition ;
870// if(fhPhotPosition)
871// delete fhPhotPosition ;
872// if(fhEMPosition)
873// delete fhEMPosition ;
874// if(fhPPSDPosition)
875// delete fhPPSDPosition ;
876
877
878 fhAllPosition = new TH2F("hAllPosition", "Position of any RP with primary photon",100, 0., 5., 100, 0., 5.);
879 fhPhotPosition = new TH2F("hPhotPosition", "Position of kGAMMA with primary photon",100, 0., 5., 100, 0., 5.);
880 fhEMPosition = new TH2F("hEMPosition", "Position of EM with primary photon", 100, 0., 5., 100, 0., 5.);
881 fhPPSDPosition = new TH2F("hPPSDPosition", "Position of PPSD with primary photon", 100, 0., 5., 100, 0., 5.);
882
883// if(fhAllReg)
884// delete fhAllReg ;
885// if(fhPhotReg)
886// delete fhPhotReg ;
887// if(fhNReg)
888// delete fhNReg ;
889// if(fhNBarReg)
890// delete fhNBarReg ;
891// if(fhChargedReg)
892// delete fhChargedReg ;
46b146ca 893
69183710 894 fhAllReg = new TH1F("hAllReg", "All primaries registered as photon", 100, 0., 5.);
895 fhPhotReg = new TH1F("hPhotReg", "Photon registered as photon", 100, 0., 5.);
896 fhNReg = new TH1F("hNReg", "N registered as photon", 100, 0., 5.);
897 fhNBarReg = new TH1F("hNBarReg", "NBar registered as photon", 100, 0., 5.);
898 fhChargedReg= new TH1F("hChargedReg", "Charged hadron registered as photon",100, 0., 5.);
46b146ca 899
69183710 900// if(fhAllEM)
901// delete fhAllEM ;
902// if(fhPhotEM)
903// delete fhPhotEM ;
904// if(fhNEM)
905// delete fhNEM ;
906// if(fhNBarEM)
907// delete fhNBarEM ;
908// if(fhChargedEM)
909// delete fhChargedEM ;
46b146ca 910
69183710 911 fhAllEM = new TH1F("hAllEM", "All primary registered as EM",100, 0., 5.);
912 fhPhotEM = new TH1F("hPhotEM", "Photon registered as EM", 100, 0., 5.);
913 fhNEM = new TH1F("hNEM", "N registered as EM", 100, 0., 5.);
914 fhNBarEM = new TH1F("hNBarEM", "NBar registered as EM", 100, 0., 5.);
915 fhChargedEM= new TH1F("hChargedEM","Charged registered as EM",100, 0., 5.);
916
917// if(fhAllPPSD)
918// delete fhAllPPSD ;
919// if(fhPhotPPSD)
920// delete fhPhotPPSD ;
921// if(fhNPPSD)
922// delete fhNPPSD ;
923// if(fhNBarPPSD)
924// delete fhNBarPPSD ;
925// if(fhChargedPPSD)
926// delete fhChargedPPSD ;
46b146ca 927
69183710 928 fhAllPPSD = new TH1F("hAllPPSD", "All primary registered as PPSD",100, 0., 5.);
929 fhPhotPPSD = new TH1F("hPhotPPSD", "Photon registered as PPSD", 100, 0., 5.);
930 fhNPPSD = new TH1F("hNPPSD", "N registered as PPSD", 100, 0., 5.);
931 fhNBarPPSD = new TH1F("hNBarPPSD", "NBar registered as PPSD", 100, 0., 5.);
932 fhChargedPPSD= new TH1F("hChargedPPSD","Charged registered as PPSD",100, 0., 5.);
933
934// if(fhPrimary)
935// delete fhPrimary ;
936 fhPrimary= new TH1F("hPrimary", "hPrimary", 100, 0., 5.);
937
938// if(fhAllRP)
939// delete fhAllRP ;
940// if(fhVeto)
941// delete fhVeto ;
942// if(fhShape)
943// delete fhShape ;
944// if(fhPPSD)
945// delete fhPPSD ;
946
947 fhAllRP = new TH1F("hAllRP","All Reconstructed particles", 100, 0., 5.);
948 fhVeto = new TH1F("hVeto", "All uncharged particles", 100, 0., 5.);
949 fhShape = new TH1F("hShape","All particles with EM shaower",100, 0., 5.);
950 fhPPSD = new TH1F("hPPSD", "All PPSD photon particles", 100, 0., 5.);
951
952
953// if(fhPhotPhot)
954// delete fhPhotPhot ;
955// if(fhPhotElec)
956// delete fhPhotElec ;
957// if(fhPhotNeuH)
958// delete fhPhotNeuH ;
959// if(fhPhotNuEM)
960// delete fhPhotNuEM ;
961// if(fhPhotChHa)
962// delete fhPhotChHa ;
963// if(fhPhotGaHa)
964// delete fhPhotGaHa ;
965
966 fhPhotPhot = new TH1F("hPhotPhot","hPhotPhot", 100, 0., 5.); //Photon registered as photon
967 fhPhotElec = new TH1F("hPhotElec","hPhotElec", 100, 0., 5.); //Photon registered as Electron
968 fhPhotNeuH = new TH1F("hPhotNeuH","hPhotNeuH", 100, 0., 5.); //Photon registered as Neutral Hadron
969 fhPhotNuEM = new TH1F("hPhotNuEM","hPhotNuEM", 100, 0., 5.); //Photon registered as Neutral EM
970 fhPhotChHa = new TH1F("hPhotChHa","hPhotChHa", 100, 0., 5.); //Photon registered as Charged Hadron
971 fhPhotGaHa = new TH1F("hPhotGaHa","hPhotGaHa", 100, 0., 5.); //Photon registered as Gamma-Hadron
fc879520 972
3fc07a80 973
92862013 974}
6ad0bfa0 975//____________________________________________________________________________
976Bool_t AliPHOSAnalyze::Init(Int_t evt)
977{
b2a60966 978 // Do a few initializations: open the root file
979 // get the AliRun object
980 // defines the clusterizer, tracksegment maker and particle identifier
981 // sets the associated parameters
6ad0bfa0 982
92862013 983 Bool_t ok = kTRUE ;
6ad0bfa0 984
985 //========== Open galice root file
986
987 if ( fRootFile == 0 ) {
988 Text_t * name = new Text_t[80] ;
989 cout << "AnalyzeOneEvent > Enter file root file name : " ;
990 cin >> name ;
92862013 991 Bool_t ok = OpenRootFile(name) ;
992 if ( !ok )
6ad0bfa0 993 cout << " AliPHOSAnalyze > Error opening " << name << endl ;
994 else {
995 //========== Get AliRun object from file
996
997 gAlice = (AliRun*) fRootFile->Get("gAlice") ;
998
999 //=========== Get the PHOS object and associated geometry from the file
1000
788dcca4 1001 fPHOS = (AliPHOSv1 *)gAlice->GetDetector("PHOS") ;
aafe457d 1002 fGeom = fPHOS->GetGeometry();
1003 // fGeom = AliPHOSGeometry::GetInstance( fPHOS->GetGeometry()->GetName(), fPHOS->GetGeometry()->GetTitle() );
83974468 1004
92862013 1005 } // else !ok
6ad0bfa0 1006 } // if fRootFile
1007
92862013 1008 if ( ok ) {
6ad0bfa0 1009
1010 //========== Create the Clusterizer
1011
774e78c2 1012 fClu = new AliPHOSClusterizerv1() ;
1013 fClu->SetEmcEnergyThreshold(0.030) ;
aafe457d 1014 fClu->SetEmcClusteringThreshold(0.20) ;
92862013 1015 fClu->SetPpsdEnergyThreshold (0.0000002) ;
1016 fClu->SetPpsdClusteringThreshold(0.0000001) ;
6ad0bfa0 1017 fClu->SetLocalMaxCut(0.03) ;
92862013 1018 fClu->SetCalibrationParameters(0., 0.00000001) ;
6ad0bfa0 1019 cout << "AnalyzeOneEvent > using clusterizer " << fClu->GetName() << endl ;
1020 fClu->PrintParameters() ;
1021
1022 //========== Creates the track segment maker
1023
1024 fTrs = new AliPHOSTrackSegmentMakerv1() ;
1025 cout << "AnalyzeOneEvent > using tack segment maker " << fTrs->GetName() << endl ;
134ce69a 1026 // fTrs->UnsetUnfoldFlag() ;
6ad0bfa0 1027
26d4b141 1028 //========== Creates the particle identifier
6ad0bfa0 1029
26d4b141 1030 fPID = new AliPHOSPIDv1() ;
1031 cout << "AnalyzeOneEvent > using particle identifier " << fPID->GetName() << endl ;
2aad621e 1032 //fPID->SetShowerProfileCuts(Float_t l1m, Float_t l1M, Float_t l2m, Float_t l2M) ;
1033 fPID->SetShowerProfileCuts(0.7, 2.0 , 0.6 , 1.5) ;
1034
6ad0bfa0 1035 //========== Creates the Reconstructioner
1036
2aad621e 1037 fRec = new AliPHOSReconstructioner(fClu, fTrs, fPID) ;
69183710 1038 fRec -> SetDebugReconstruction(kFALSE);
6ad0bfa0 1039
1040 //=========== Connect the various Tree's for evt
1041
1042 if ( evt == -999 ) {
1043 cout << "AnalyzeOneEvent > Enter event number : " ;
1044 cin >> evt ;
1045 cout << evt << endl ;
1046 }
1047 fEvt = evt ;
1048 gAlice->GetEvent(evt);
1049
1050 //=========== Get the Digit TTree
1051
1052 gAlice->TreeD()->GetEvent(0) ;
1053
92862013 1054 } // ok
6ad0bfa0 1055
92862013 1056 return ok ;
6ad0bfa0 1057}
1058
92862013 1059
6ad0bfa0 1060//____________________________________________________________________________
92862013 1061void AliPHOSAnalyze::DisplayKineEvent(Int_t evt)
6ad0bfa0 1062{
b2a60966 1063 // Display particles from the Kine Tree in global Alice (theta, phi) coordinates.
1064 // One PHOS module at the time.
1065 // The particle type can be selected.
1066
6ad0bfa0 1067 if (evt == -999)
1068 evt = fEvt ;
1069
92862013 1070 Int_t module ;
6ad0bfa0 1071 cout << "DisplayKineEvent > which module (1-5, -1: all) ? " ;
92862013 1072 cin >> module ; cout << module << endl ;
6ad0bfa0 1073
92862013 1074 Int_t testparticle ;
6ad0bfa0 1075 cout << " 22 : PHOTON " << endl
1076 << " (-)11 : (POSITRON)ELECTRON " << endl
1077 << " (-)2112 : (ANTI)NEUTRON " << endl
1078 << " -999 : Everything else " << endl ;
1079 cout << "DisplayKineEvent > enter PDG particle code to display " ;
92862013 1080 cin >> testparticle ; cout << testparticle << endl ;
6ad0bfa0 1081
92862013 1082 Text_t histoname[80] ;
1083 sprintf(histoname,"Event %d: Incident particles in module %d", evt, module) ;
6ad0bfa0 1084
92862013 1085 Double_t tm, tM, pm, pM ; // min and Max theta and phi covered by module
cf0c2bc1 1086 fGeom->EmcModuleCoverage(module, tm, tM, pm, pM, AliPHOSGeometry::Degre() ) ;
6ad0bfa0 1087
1088 Double_t theta, phi ;
cf0c2bc1 1089 fGeom->EmcXtalCoverage(theta, phi, AliPHOSGeometry::Degre() ) ;
6ad0bfa0 1090
1091 Int_t tdim = (Int_t)( (tM - tm) / theta ) ;
1092 Int_t pdim = (Int_t)( (pM - pm) / phi ) ;
1093
1094 tm -= theta ;
1095 tM += theta ;
1096 pm -= phi ;
1097 pM += phi ;
1098
92862013 1099 TH2F * histoparticle = new TH2F("histoparticle", histoname,
6ad0bfa0 1100 pdim, pm, pM, tdim, tm, tM) ;
92862013 1101 histoparticle->SetStats(kFALSE) ;
6ad0bfa0 1102
1103 // Get pointers to Alice Particle TClonesArray
1104
92862013 1105 TParticle * particle;
1106 TClonesArray * particlearray = gAlice->Particles();
6ad0bfa0 1107
92862013 1108 Text_t canvasname[80];
1109 sprintf(canvasname,"Particles incident in PHOS/EMC module # %d",module) ;
1110 TCanvas * kinecanvas = new TCanvas("kinecanvas", canvasname, 650, 500) ;
6ad0bfa0 1111
1112 // get the KINE Tree
1113
92862013 1114 TTree * kine = gAlice->TreeK() ;
1115 Stat_t nParticles = kine->GetEntries() ;
1116 cout << "DisplayKineEvent > events in kine " << nParticles << endl ;
6ad0bfa0 1117
1118 // loop over particles
1119
92862013 1120 Double_t kRADDEG = 180. / TMath::Pi() ;
6ad0bfa0 1121 Int_t index1 ;
1122 Int_t nparticlein = 0 ;
92862013 1123 for (index1 = 0 ; index1 < nParticles ; index1++){
1124 Int_t nparticle = particlearray->GetEntriesFast() ;
6ad0bfa0 1125 Int_t index2 ;
1126 for( index2 = 0 ; index2 < nparticle ; index2++) {
92862013 1127 particle = (TParticle*)particlearray->UncheckedAt(index2) ;
1128 Int_t particletype = particle->GetPdgCode() ;
1129 if (testparticle == -999 || testparticle == particletype) {
1130 Double_t phi = particle->Phi() ;
1131 Double_t theta = particle->Theta() ;
6ad0bfa0 1132 Int_t mod ;
1133 Double_t x, z ;
92862013 1134 fGeom->ImpactOnEmc(theta, phi, mod, z, x) ;
1135 if ( mod == module ) {
6ad0bfa0 1136 nparticlein++ ;
b2a60966 1137 if (particle->Energy() > fClu->GetEmcClusteringThreshold() )
1138 histoparticle->Fill(phi*kRADDEG, theta*kRADDEG, particle->Energy() ) ;
6ad0bfa0 1139 }
1140 }
1141 }
1142 }
92862013 1143 kinecanvas->Draw() ;
1144 histoparticle->Draw("color") ;
1145 TPaveText * pavetext = new TPaveText(294, 100, 300, 101);
6ad0bfa0 1146 Text_t text[40] ;
1147 sprintf(text, "Particles: %d ", nparticlein) ;
92862013 1148 pavetext->AddText(text) ;
1149 pavetext->Draw() ;
1150 kinecanvas->Update();
6ad0bfa0 1151
1152}
1153//____________________________________________________________________________
1154void AliPHOSAnalyze::DisplayRecParticles()
1155{
b2a60966 1156 // Display reconstructed particles in global Alice(theta, phi) coordinates.
1157 // One PHOS module at the time.
1158 // Click on symbols indicate the reconstructed particle type.
1159
6ad0bfa0 1160 if (fEvt == -999) {
15605d3c 1161 cout << "DisplayRecParticles > Analyze an event first ... (y/n) " ;
6ad0bfa0 1162 Text_t answer[1] ;
1163 cin >> answer ; cout << answer ;
46b146ca 1164// if ( answer == "y" )
1165// AnalyzeOneEvent() ;
6ad0bfa0 1166 }
1167 if (fEvt != -999) {
1168
92862013 1169 Int_t module ;
15605d3c 1170 cout << "DisplayRecParticles > which module (1-5, -1: all) ? " ;
92862013 1171 cin >> module ; cout << module << endl ;
1172 Text_t histoname[80] ;
1173 sprintf(histoname,"Event %d: Reconstructed particles in module %d", fEvt, module) ;
1174 Double_t tm, tM, pm, pM ; // min and Max theta and phi covered by module
cf0c2bc1 1175 fGeom->EmcModuleCoverage(module, tm, tM, pm, pM, AliPHOSGeometry::Degre() ) ;
6ad0bfa0 1176 Double_t theta, phi ;
cf0c2bc1 1177 fGeom->EmcXtalCoverage(theta, phi, AliPHOSGeometry::Degre() ) ;
6ad0bfa0 1178 Int_t tdim = (Int_t)( (tM - tm) / theta ) ;
1179 Int_t pdim = (Int_t)( (pM - pm) / phi ) ;
1180 tm -= theta ;
1181 tM += theta ;
1182 pm -= phi ;
92862013 1183 TH2F * histoRparticle = new TH2F("histoRparticle", histoname,
6ad0bfa0 1184 pdim, pm, pM, tdim, tm, tM) ;
92862013 1185 histoRparticle->SetStats(kFALSE) ;
1186 Text_t canvasname[80] ;
1187 sprintf(canvasname, "Reconstructed particles in PHOSmodule # %d", module) ;
1188 TCanvas * rparticlecanvas = new TCanvas("RparticleCanvas", canvasname, 650, 500) ;
fc879520 1189 AliPHOSRecParticle::RecParticlesList * rpl = *fPHOS->RecParticles() ;
92862013 1190 Int_t nRecParticles = rpl->GetEntries() ;
1191 Int_t nRecParticlesInModule = 0 ;
6ad0bfa0 1192 TIter nextRecPart(rpl) ;
1193 AliPHOSRecParticle * rp ;
92862013 1194 cout << "DisplayRecParticles > " << nRecParticles << " reconstructed particles " << endl ;
1195 Double_t kRADDEG = 180. / TMath::Pi() ;
6ad0bfa0 1196 while ( (rp = (AliPHOSRecParticle *)nextRecPart() ) ) {
1197 AliPHOSTrackSegment * ts = rp->GetPHOSTrackSegment() ;
b2a60966 1198 if ( ts->GetPHOSMod() == module ) {
1199 Int_t numberofprimaries = 0 ;
134ce69a 1200 Int_t * listofprimaries = 0;
1201 rp->GetPrimaries(numberofprimaries) ;
b2a60966 1202 cout << "Number of primaries = " << numberofprimaries << endl ;
1203 Int_t index ;
1204 for ( index = 0 ; index < numberofprimaries ; index++)
1205 cout << " primary # " << index << " = " << listofprimaries[index] << endl ;
1206
92862013 1207 nRecParticlesInModule++ ;
1208 Double_t theta = rp->Theta() * kRADDEG ;
1209 Double_t phi = rp->Phi() * kRADDEG ;
6ad0bfa0 1210 Double_t energy = rp->Energy() ;
92862013 1211 histoRparticle->Fill(phi, theta, energy) ;
6ad0bfa0 1212 }
1213 }
92862013 1214 histoRparticle->Draw("color") ;
15605d3c 1215
1216 nextRecPart.Reset() ;
1217 while ( (rp = (AliPHOSRecParticle *)nextRecPart() ) ) {
1218 AliPHOSTrackSegment * ts = rp->GetPHOSTrackSegment() ;
1219 if ( ts->GetPHOSMod() == module )
1220 rp->Draw("P") ;
1221 }
1222
6ad0bfa0 1223 Text_t text[80] ;
92862013 1224 sprintf(text, "reconstructed particles: %d", nRecParticlesInModule) ;
1225 TPaveText * pavetext = new TPaveText(292, 100, 300, 101);
1226 pavetext->AddText(text) ;
1227 pavetext->Draw() ;
1228 rparticlecanvas->Update() ;
6ad0bfa0 1229 }
1230}
1231
1232//____________________________________________________________________________
1233void AliPHOSAnalyze::DisplayRecPoints()
1234{
b2a60966 1235 // Display reconstructed points in local PHOS-module (x, z) coordinates.
1236 // One PHOS module at the time.
1237 // Click on symbols displays the EMC cluster, or PPSD information.
1238
6ad0bfa0 1239 if (fEvt == -999) {
1240 cout << "DisplayRecPoints > Analyze an event first ... (y/n) " ;
1241 Text_t answer[1] ;
1242 cin >> answer ; cout << answer ;
46b146ca 1243// if ( answer == "y" )
1244// AnalyzeOneEvent() ;
6ad0bfa0 1245 }
1246 if (fEvt != -999) {
1247
92862013 1248 Int_t module ;
6ad0bfa0 1249 cout << "DisplayRecPoints > which module (1-5, -1: all) ? " ;
92862013 1250 cin >> module ; cout << module << endl ;
6ad0bfa0 1251
92862013 1252 Text_t canvasname[80];
1253 sprintf(canvasname,"Digits in PHOS/EMC module # %d",module) ;
1254 TCanvas * modulecanvas = new TCanvas("module", canvasname, 650, 500) ;
1255 modulecanvas->Draw() ;
6ad0bfa0 1256
92862013 1257 //=========== Creating 2d-histogram of the PHOS module
6ad0bfa0 1258 // a little bit junkie but is used to test Geom functinalities
1259
92862013 1260 Double_t tm, tM, pm, pM ; // min and Max theta and phi covered by module
6ad0bfa0 1261
92862013 1262 fGeom->EmcModuleCoverage(module, tm, tM, pm, pM);
6ad0bfa0 1263 // convert angles into coordinates local to the EMC module of interest
1264
92862013 1265 Int_t emcModuleNumber ;
1266 Double_t emcModulexm, emcModulezm ; // minimum local coordinate in a given EMCA module
1267 Double_t emcModulexM, emcModulezM ; // maximum local coordinate in a given EMCA module
1268 fGeom->ImpactOnEmc(tm, pm, emcModuleNumber, emcModulezm, emcModulexm) ;
1269 fGeom->ImpactOnEmc(tM, pM, emcModuleNumber, emcModulezM, emcModulexM) ;
1270 Int_t xdim = (Int_t)( ( emcModulexM - emcModulexm ) / fGeom->GetCrystalSize(0) ) ;
1271 Int_t zdim = (Int_t)( ( emcModulezM - emcModulezm ) / fGeom->GetCrystalSize(2) ) ;
1272 Float_t xmin = emcModulexm - fGeom->GetCrystalSize(0) ;
1273 Float_t xMax = emcModulexM + fGeom->GetCrystalSize(0) ;
1274 Float_t zmin = emcModulezm - fGeom->GetCrystalSize(2) ;
1275 Float_t zMax = emcModulezM + fGeom->GetCrystalSize(2) ;
1276 Text_t histoname[80];
1277 sprintf(histoname,"Event %d: Digits and RecPoints in module %d", fEvt, module) ;
1278 TH2F * hModule = new TH2F("HistoReconstructed", histoname,
6ad0bfa0 1279 xdim, xmin, xMax, zdim, zmin, zMax) ;
1280 hModule->SetMaximum(2.0);
1281 hModule->SetMinimum(0.0);
1282 hModule->SetStats(kFALSE);
1283
1284 TIter next(fPHOS->Digits()) ;
92862013 1285 Float_t energy, y, z;
1286 Float_t etot=0.;
1287 Int_t relid[4]; Int_t nDigits = 0 ;
6ad0bfa0 1288 AliPHOSDigit * digit ;
78ccc411 1289
1290 // Making 2D histogram of the EMC module
6ad0bfa0 1291 while((digit = (AliPHOSDigit *)next()))
1292 {
92862013 1293 fGeom->AbsToRelNumbering(digit->GetId(), relid) ;
78ccc411 1294 if (relid[0] == module && relid[1] == 0)
5959e315 1295 {
92862013 1296 energy = fClu->Calibrate(digit->GetAmp()) ;
774e78c2 1297 cout << "Energy is " << energy << " and threshold is " << fClu->GetEmcEnergyThreshold() << endl;
78ccc411 1298 if (energy > fClu->GetEmcEnergyThreshold() ){
5959e315 1299 nDigits++ ;
1300 etot += energy ;
1301 fGeom->RelPosInModule(relid,y,z) ;
1302 hModule->Fill(y, z, energy) ;
3fc07a80 1303 }
6ad0bfa0 1304 }
1305 }
92862013 1306 cout <<"DrawRecPoints > Found in module "
1307 << module << " " << nDigits << " digits with total energy " << etot << endl ;
6ad0bfa0 1308 hModule->Draw("col2") ;
1309
92862013 1310 //=========== Cluster in module
6ad0bfa0 1311
83974468 1312 // TClonesArray * emcRP = fPHOS->EmcClusters() ;
fc879520 1313 TObjArray * emcRP = *(fPHOS->EmcRecPoints()) ;
83974468 1314
92862013 1315 etot = 0.;
1316 Int_t totalnClusters = 0 ;
1317 Int_t nClusters = 0 ;
1318 TIter nextemc(emcRP) ;
6ad0bfa0 1319 AliPHOSEmcRecPoint * emc ;
1320 while((emc = (AliPHOSEmcRecPoint *)nextemc()))
1321 {
3fc07a80 1322 // Int_t numberofprimaries ;
1323 // Int_t * primariesarray = new Int_t[10] ;
1324 // emc->GetPrimaries(numberofprimaries, primariesarray) ;
92862013 1325 totalnClusters++ ;
1326 if ( emc->GetPHOSMod() == module )
6ad0bfa0 1327 {
92862013 1328 nClusters++ ;
1329 energy = emc->GetTotalEnergy() ;
1330 etot+= energy ;
26d4b141 1331 emc->Draw("M") ;
6ad0bfa0 1332 }
1333 }
92862013 1334 cout << "DrawRecPoints > Found " << totalnClusters << " EMC Clusters in PHOS" << endl ;
1335 cout << "DrawRecPoints > Found in module " << module << " " << nClusters << " EMC Clusters " << endl ;
1336 cout << "DrawRecPoints > total energy " << etot << endl ;
6ad0bfa0 1337
92862013 1338 TPaveText * pavetext = new TPaveText(22, 80, 83, 90);
6ad0bfa0 1339 Text_t text[40] ;
92862013 1340 sprintf(text, "digits: %d; clusters: %d", nDigits, nClusters) ;
1341 pavetext->AddText(text) ;
1342 pavetext->Draw() ;
1343 modulecanvas->Update();
6ad0bfa0 1344
92862013 1345 //=========== Cluster in module PPSD Down
6ad0bfa0 1346
83974468 1347 // TClonesArray * ppsdRP = fPHOS->PpsdClusters() ;
fc879520 1348 TObjArray * ppsdRP = *(fPHOS->PpsdRecPoints() );
83974468 1349
92862013 1350 etot = 0.;
1351 TIter nextPpsd(ppsdRP) ;
1352 AliPHOSPpsdRecPoint * ppsd ;
1353 while((ppsd = (AliPHOSPpsdRecPoint *)nextPpsd()))
6ad0bfa0 1354 {
92862013 1355 totalnClusters++ ;
1356 if ( ppsd->GetPHOSMod() == module )
6ad0bfa0 1357 {
92862013 1358 nClusters++ ;
1359 energy = ppsd->GetEnergy() ;
1360 etot+=energy ;
1361 if (!ppsd->GetUp()) ppsd->Draw("P") ;
6ad0bfa0 1362 }
1363 }
92862013 1364 cout << "DrawRecPoints > Found " << totalnClusters << " Ppsd Down Clusters in PHOS" << endl ;
1365 cout << "DrawRecPoints > Found in module " << module << " " << nClusters << " Ppsd Down Clusters " << endl ;
1366 cout << "DrawRecPoints > total energy " << etot << endl ;
6ad0bfa0 1367
92862013 1368 //=========== Cluster in module PPSD Up
6ad0bfa0 1369
fc879520 1370 ppsdRP = *(fPHOS->PpsdRecPoints()) ;
83974468 1371
92862013 1372 etot = 0.;
1373 TIter nextPpsdUp(ppsdRP) ;
1374 while((ppsd = (AliPHOSPpsdRecPoint *)nextPpsdUp()))
6ad0bfa0 1375 {
92862013 1376 totalnClusters++ ;
1377 if ( ppsd->GetPHOSMod() == module )
6ad0bfa0 1378 {
92862013 1379 nClusters++ ;
1380 energy = ppsd->GetEnergy() ;
1381 etot+=energy ;
1382 if (ppsd->GetUp()) ppsd->Draw("P") ;
6ad0bfa0 1383 }
1384 }
92862013 1385 cout << "DrawRecPoints > Found " << totalnClusters << " Ppsd Up Clusters in PHOS" << endl ;
1386 cout << "DrawRecPoints > Found in module " << module << " " << nClusters << " Ppsd Up Clusters " << endl ;
1387 cout << "DrawRecPoints > total energy " << etot << endl ;
6ad0bfa0 1388
1389 } // if !-999
1390}
1391
1392//____________________________________________________________________________
1393void AliPHOSAnalyze::DisplayTrackSegments()
1394{
b2a60966 1395 // Display track segments in local PHOS-module (x, z) coordinates.
1396 // One PHOS module at the time.
1397 // One symbol per PHOS subsystem: EMC, upper PPSD, lower PPSD.
1398
6ad0bfa0 1399 if (fEvt == -999) {
1400 cout << "DisplayTrackSegments > Analyze an event first ... (y/n) " ;
1401 Text_t answer[1] ;
1402 cin >> answer ; cout << answer ;
46b146ca 1403// if ( answer == "y" )
1404// AnalyzeOneEvent() ;
6ad0bfa0 1405 }
1406 if (fEvt != -999) {
1407
92862013 1408 Int_t module ;
6ad0bfa0 1409 cout << "DisplayTrackSegments > which module (1-5, -1: all) ? " ;
92862013 1410 cin >> module ; cout << module << endl ;
1411 //=========== Creating 2d-histogram of the PHOS module
6ad0bfa0 1412 // a little bit junkie but is used to test Geom functinalities
1413
92862013 1414 Double_t tm, tM, pm, pM ; // min and Max theta and phi covered by module
6ad0bfa0 1415
92862013 1416 fGeom->EmcModuleCoverage(module, tm, tM, pm, pM);
6ad0bfa0 1417 // convert angles into coordinates local to the EMC module of interest
1418
92862013 1419 Int_t emcModuleNumber ;
1420 Double_t emcModulexm, emcModulezm ; // minimum local coordinate in a given EMCA module
1421 Double_t emcModulexM, emcModulezM ; // maximum local coordinate in a given EMCA module
1422 fGeom->ImpactOnEmc(tm, pm, emcModuleNumber, emcModulezm, emcModulexm) ;
1423 fGeom->ImpactOnEmc(tM, pM, emcModuleNumber, emcModulezM, emcModulexM) ;
1424 Int_t xdim = (Int_t)( ( emcModulexM - emcModulexm ) / fGeom->GetCrystalSize(0) ) ;
1425 Int_t zdim = (Int_t)( ( emcModulezM - emcModulezm ) / fGeom->GetCrystalSize(2) ) ;
1426 Float_t xmin = emcModulexm - fGeom->GetCrystalSize(0) ;
1427 Float_t xMax = emcModulexM + fGeom->GetCrystalSize(0) ;
1428 Float_t zmin = emcModulezm - fGeom->GetCrystalSize(2) ;
1429 Float_t zMax = emcModulezM + fGeom->GetCrystalSize(2) ;
1430 Text_t histoname[80];
1431 sprintf(histoname,"Event %d: Track Segments in module %d", fEvt, module) ;
1432 TH2F * histotrack = new TH2F("histotrack", histoname,
6ad0bfa0 1433 xdim, xmin, xMax, zdim, zmin, zMax) ;
92862013 1434 histotrack->SetStats(kFALSE);
1435 Text_t canvasname[80];
1436 sprintf(canvasname,"Track segments in PHOS/EMC-PPSD module # %d", module) ;
1437 TCanvas * trackcanvas = new TCanvas("TrackSegmentCanvas", canvasname, 650, 500) ;
1438 histotrack->Draw() ;
6ad0bfa0 1439
fc879520 1440 AliPHOSTrackSegment::TrackSegmentsList * trsegl = *(fPHOS->TrackSegments()) ;
6ad0bfa0 1441 AliPHOSTrackSegment * trseg ;
1442
92862013 1443 Int_t nTrackSegments = trsegl->GetEntries() ;
6ad0bfa0 1444 Int_t index ;
92862013 1445 Float_t etot = 0 ;
1446 Int_t nTrackSegmentsInModule = 0 ;
1447 for(index = 0; index < nTrackSegments ; index++){
6ad0bfa0 1448 trseg = (AliPHOSTrackSegment * )trsegl->At(index) ;
92862013 1449 etot+= trseg->GetEnergy() ;
1450 if ( trseg->GetPHOSMod() == module ) {
1451 nTrackSegmentsInModule++ ;
6ad0bfa0 1452 trseg->Draw("P");
1453 }
1454 }
1455 Text_t text[80] ;
92862013 1456 sprintf(text, "track segments: %d", nTrackSegmentsInModule) ;
1457 TPaveText * pavetext = new TPaveText(22, 80, 83, 90);
1458 pavetext->AddText(text) ;
1459 pavetext->Draw() ;
1460 trackcanvas->Update() ;
1461 cout << "DisplayTrackSegments > Found " << trsegl->GetEntries() << " Track segments with total energy "<< etot << endl ;
6ad0bfa0 1462
1463 }
1464}
1465//____________________________________________________________________________
1466Bool_t AliPHOSAnalyze::OpenRootFile(Text_t * name)
1467{
b2a60966 1468 // Open the root file named "name"
1469
1470 fRootFile = new TFile(name, "update") ;
6ad0bfa0 1471 return fRootFile->IsOpen() ;
1472}
92862013 1473//____________________________________________________________________________
46b146ca 1474void AliPHOSAnalyze::SaveHistograms()
134ce69a 1475{
1476 // Saves the histograms in a root file named "name.analyzed"
1477
1478 Text_t outputname[80] ;
1479 sprintf(outputname,"%s.analyzed",fRootFile->GetName());
1480 TFile output(outputname,"RECREATE");
1481 output.cd();
46b146ca 1482
69183710 1483 if (fhAllEnergy)
1484 fhAllEnergy->Write() ;
1485 if (fhPhotEnergy)
1486 fhPhotEnergy->Write() ;
1487 if(fhEMEnergy)
1488 fhEMEnergy->Write() ;
1489 if(fhPPSDEnergy)
1490 fhPPSDEnergy->Write() ;
1491 if(fhAllPosition)
1492 fhAllPosition->Write() ;
1493 if(fhPhotPosition)
1494 fhPhotPosition->Write() ;
1495 if(fhEMPosition)
1496 fhEMPosition->Write() ;
1497 if(fhPPSDPosition)
1498 fhPPSDPosition->Write() ;
fc879520 1499 if (fhAllReg)
1500 fhAllReg->Write() ;
69183710 1501 if (fhPhotReg)
1502 fhPhotReg->Write() ;
fc879520 1503 if(fhNReg)
1504 fhNReg->Write() ;
1505 if(fhNBarReg)
1506 fhNBarReg->Write() ;
1507 if(fhChargedReg)
1508 fhChargedReg->Write() ;
fc879520 1509 if (fhAllEM)
1510 fhAllEM->Write() ;
69183710 1511 if (fhPhotEM)
1512 fhPhotEM->Write() ;
fc879520 1513 if(fhNEM)
1514 fhNEM->Write() ;
1515 if(fhNBarEM)
1516 fhNBarEM->Write() ;
1517 if(fhChargedEM)
1518 fhChargedEM->Write() ;
69183710 1519 if (fhAllPPSD)
1520 fhAllPPSD->Write() ;
1521 if (fhPhotPPSD)
1522 fhPhotPPSD->Write() ;
1523 if(fhNPPSD)
1524 fhNPPSD->Write() ;
1525 if(fhNBarPPSD)
1526 fhNBarPPSD->Write() ;
1527 if(fhChargedPPSD)
1528 fhChargedPPSD->Write() ;
fc879520 1529 if(fhPrimary)
1530 fhPrimary->Write() ;
69183710 1531 if(fhAllRP)
1532 fhAllRP->Write() ;
1533 if(fhVeto)
1534 fhVeto->Write() ;
1535 if(fhShape)
1536 fhShape->Write() ;
1537 if(fhPPSD)
1538 fhPPSD->Write() ;
fc879520 1539 if(fhPhotPhot)
1540 fhPhotPhot->Write() ;
1541 if(fhPhotElec)
1542 fhPhotElec->Write() ;
1543 if(fhPhotNeuH)
1544 fhPhotNeuH->Write() ;
1545 if(fhPhotNuEM)
1546 fhPhotNuEM->Write() ;
1547 if(fhPhotNuEM)
1548 fhPhotNuEM->Write() ;
1549 if(fhPhotChHa)
1550 fhPhotChHa->Write() ;
1551 if(fhPhotGaHa)
1552 fhPhotGaHa->Write() ;
46b146ca 1553 if(fhEnergyCorrelations)
1554 fhEnergyCorrelations->Write() ;
1555
92862013 1556 output.Write();
1557 output.Close();
1558}
69183710 1559//____________________________________________________________________________
1560Float_t AliPHOSAnalyze::CorrectEnergy(Float_t ERecPart)
1561{
1562 return ERecPart/0.8783 ;
1563}
1564
46b146ca 1565//____________________________________________________________________________
1566void AliPHOSAnalyze::ResetHistograms()
1567{
1568 fhEnergyCorrelations = 0 ; //Energy correlations between Eloss in Convertor and PPSD(2)
1569
1570 fhEmcDigit = 0 ; // Histo of digit energies in the Emc
1571 fhVetoDigit = 0 ; // Histo of digit energies in the Veto
1572 fhConvertorDigit = 0 ; // Histo of digit energies in the Convertor
1573 fhEmcCluster = 0 ; // Histo of Cluster energies in Emc
1574 fhVetoCluster = 0 ; // Histo of Cluster energies in Veto
1575 fhConvertorCluster = 0 ; // Histo of Cluster energies in Convertor
1576 fhConvertorEmc = 0 ; // 2d Convertor versus Emc energies
1577
69183710 1578 fhAllEnergy = 0 ;
1579 fhPhotEnergy = 0 ; // Total spectrum of detected photons
1580 fhEMEnergy = 0 ; // Spectrum of detected electrons with electron primary
1581 fhPPSDEnergy = 0 ;
1582 fhAllPosition = 0 ;
1583 fhPhotPosition = 0 ;
1584 fhEMPosition = 0 ;
1585 fhPPSDPosition = 0 ;
1586
1587 fhPhotReg = 0 ;
46b146ca 1588 fhAllReg = 0 ;
1589 fhNReg = 0 ;
1590 fhNBarReg = 0 ;
1591 fhChargedReg = 0 ;
69183710 1592 fhPhotEM = 0 ;
46b146ca 1593 fhAllEM = 0 ;
1594 fhNEM = 0 ;
1595 fhNBarEM = 0 ;
1596 fhChargedEM = 0 ;
69183710 1597 fhPhotPPSD = 0 ;
1598 fhAllPPSD = 0 ;
1599 fhNPPSD = 0 ;
1600 fhNBarPPSD = 0 ;
1601 fhChargedPPSD = 0 ;
1602
46b146ca 1603 fhPrimary = 0 ;
1604
1605 fhPhotPhot = 0 ;
1606 fhPhotElec = 0 ;
1607 fhPhotNeuH = 0 ;
1608 fhPhotNuEM = 0 ;
1609 fhPhotChHa = 0 ;
1610 fhPhotGaHa = 0 ;
1611
134ce69a 1612
46b146ca 1613}
134ce69a 1614
1615