]> 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
Minor protection in GetPHOSMod()
[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"
c3b9b3f9 37#include "TStyle.h"
6ad0bfa0 38
39// --- Standard library ---
40
de9ec31b 41#include <iostream.h>
42#include <stdio.h>
92862013 43
6ad0bfa0 44// --- AliRoot header files ---
45
46#include "AliRun.h"
47#include "AliPHOSAnalyze.h"
48#include "AliPHOSClusterizerv1.h"
49#include "AliPHOSTrackSegmentMakerv1.h"
26d4b141 50#include "AliPHOSPIDv1.h"
6ad0bfa0 51#include "AliPHOSReconstructioner.h"
52#include "AliPHOSDigit.h"
53#include "AliPHOSTrackSegment.h"
54#include "AliPHOSRecParticle.h"
83974468 55#include "AliPHOSIndexToObject.h"
c3b9b3f9 56#include "AliPHOSCPV.h"
6ad0bfa0 57
58ClassImp(AliPHOSAnalyze)
59
60
61//____________________________________________________________________________
62 AliPHOSAnalyze::AliPHOSAnalyze()
63{
b2a60966 64 // default ctor (useless)
6ad0bfa0 65
66 fRootFile = 0 ;
67}
68
69//____________________________________________________________________________
70AliPHOSAnalyze::AliPHOSAnalyze(Text_t * name)
71{
83974468 72 // ctor: analyze events from root file "name"
6ad0bfa0 73
92862013 74 Bool_t ok = OpenRootFile(name) ;
75 if ( !ok ) {
6ad0bfa0 76 cout << " AliPHOSAnalyze > Error opening " << name << endl ;
77 }
78 else {
eecb6765 79 //========== Get AliRun object from file
80 gAlice = (AliRun*) fRootFile->Get("gAlice") ;
fc879520 81
eecb6765 82 //=========== Get the PHOS object and associated geometry from the file
788dcca4 83 fPHOS = (AliPHOSv1 *)gAlice->GetDetector("PHOS") ;
eecb6765 84 fGeom = AliPHOSGeometry::GetInstance( fPHOS->GetGeometry()->GetName(), fPHOS->GetGeometry()->GetTitle() );
fc879520 85
eecb6765 86 //========== Initializes the Index to Object converter
87 fObjGetter = AliPHOSIndexToObject::GetInstance(fPHOS) ;
88 //========== Current event number
89 fEvt = -999 ;
fc879520 90
6ad0bfa0 91 }
c3b9b3f9 92 fDebugLevel = 0;
46b146ca 93 ResetHistograms() ;
6ad0bfa0 94}
95
88714635 96//____________________________________________________________________________
97AliPHOSAnalyze::AliPHOSAnalyze(const AliPHOSAnalyze & ana)
98{
99 // copy ctor
100 ( (AliPHOSAnalyze &)ana ).Copy(*this) ;
101}
102
103//____________________________________________________________________________
191c1c41 104void AliPHOSAnalyze::Copy(TObject & obj)
88714635 105{
106 // copy an analysis into an other one
107 TObject::Copy(obj) ;
108 // I do nothing more because the copy is silly but the Code checkers requires one
109}
110
6ad0bfa0 111//____________________________________________________________________________
112AliPHOSAnalyze::~AliPHOSAnalyze()
113{
114 // dtor
115
83974468 116 if (fRootFile->IsOpen() )
117 fRootFile->Close() ;
c3b9b3f9 118 if(fRootFile) {delete fRootFile ; fRootFile=0 ;}
119 if(fPHOS) {delete fPHOS ; fPHOS =0 ;}
120 if(fClu) {delete fClu ; fClu =0 ;}
121 if(fPID) {delete fPID ; fPID =0 ;}
122 if(fRec) {delete fRec ; fRec =0 ;}
123 if(fTrs) {delete fTrs ; fTrs =0 ;}
6ad0bfa0 124
125}
126
127//____________________________________________________________________________
46b146ca 128void AliPHOSAnalyze::ActivePPSD(Int_t Nevents=1){
6ad0bfa0 129
46b146ca 130 fhEnergyCorrelations = new TH2F("hEnergyCorrelations","hEnergyCorrelations",40, 0., 0.15, 30, 0., 3.e-5);
131 //========== Create the Clusterizer
132 fClu = new AliPHOSClusterizerv1() ;
133 fClu->SetEmcEnergyThreshold(0.05) ;
134 fClu->SetEmcClusteringThreshold(0.20) ;
135 fClu->SetPpsdEnergyThreshold (0.0000002) ;
136 fClu->SetPpsdClusteringThreshold(0.0000001) ;
137 fClu->SetLocalMaxCut(0.03) ;
138 fClu->SetCalibrationParameters(0., 0.00000001) ;
b2a60966 139
46b146ca 140 Int_t ievent;
83974468 141
46b146ca 142 for ( ievent=0; ievent<Nevents; ievent++)
143 {
144
145 //========== Event Number>
146 if ( ( log10((Float_t)(ievent+1)) - (Int_t)(log10((Float_t)(ievent+1))) ) == 0. )
147 cout << "AnalyzeResolutions > " << "Event is " << ievent << endl ;
148
149 //=========== Connects the various Tree's for evt
150 gAlice->GetEvent(ievent);
b2a60966 151
46b146ca 152 //=========== Gets the Kine TTree
153 gAlice->TreeK()->GetEvent(0) ;
c3b9b3f9 154
46b146ca 155 //=========== Get the Digit Tree
156 gAlice->TreeD()->GetEvent(0) ;
157
158 //========== Creating branches ===================================
159 AliPHOSRecPoint::RecPointsList ** EmcRecPoints = fPHOS->EmcRecPoints() ;
160 gAlice->TreeR()->SetBranchAddress( "PHOSEmcRP", EmcRecPoints ) ;
161
162 AliPHOSRecPoint::RecPointsList ** PpsdRecPoints = fPHOS->PpsdRecPoints() ;
163 gAlice->TreeR()->SetBranchAddress( "PHOSPpsdRP", PpsdRecPoints ) ;
164
165 AliPHOSTrackSegment::TrackSegmentsList ** TrackSegmentsList = fPHOS->TrackSegments() ;
166 if( (*TrackSegmentsList) )
167 (*TrackSegmentsList)->Clear() ;
168 gAlice->TreeR()->SetBranchAddress( "PHOSTS", TrackSegmentsList ) ;
169
170 AliPHOSRecParticle::RecParticlesList ** RecParticleList = fPHOS->RecParticles() ;
171 if( (*RecParticleList) )
172 (*RecParticleList)->Clear() ;
173 gAlice->TreeR()->SetBranchAddress( "PHOSRP", RecParticleList ) ;
174
175
176 //=========== Gets the Reconstraction TTree
177 gAlice->TreeR()->GetEvent(0) ;
178
179 AliPHOSPpsdRecPoint * RecPoint ;
180 Int_t relid[4] ;
181 TIter nextRP(*fPHOS->PpsdRecPoints() ) ;
182 while( ( RecPoint = (AliPHOSPpsdRecPoint *)nextRP() ) )
183 {
184 if(!(RecPoint->GetUp()) ) {
185 AliPHOSDigit *digit ;
186 Int_t iDigit ;
187 for(iDigit = 0; iDigit < fPHOS->Digits()->GetEntries(); iDigit++)
188 {
189 digit = (AliPHOSDigit *) fPHOS->Digits()->At(iDigit) ;
190 fGeom->AbsToRelNumbering(digit->GetId(), relid) ;
191 if((relid[2]==1)&&(relid[3]==1)&&(relid[0]==RecPoint->GetPHOSMod())){
192 Float_t ConvertorEnergy = fClu->Calibrate(digit->GetAmp()) ;
193 fhEnergyCorrelations->Fill(ConvertorEnergy,RecPoint->GetTotalEnergy() );
194 break ;
195 }
196 }
197 break ;
198 }
199 }
200 }
201 SaveHistograms() ;
202 fhEnergyCorrelations->Draw("BOX") ;
6ad0bfa0 203}
204
92862013 205//____________________________________________________________________________
c3b9b3f9 206void AliPHOSAnalyze::AnalyzeManyEvents(Int_t Nevents, Int_t module)
92862013 207{
b2a60966 208 // analyzes Nevents events in a single PHOS module
46b146ca 209 // Events should be reconstructed by Reconstruct()
92862013 210
211 if ( fRootFile == 0 )
212 cout << "AnalyzeManyEvents > " << "Root File not openned" << endl ;
213 else
214 {
92862013 215 //========== Booking Histograms
216 cout << "AnalyzeManyEvents > " << "Booking Histograms" << endl ;
217 BookingHistograms();
46b146ca 218
92862013 219 Int_t ievent;
6a3f1304 220 Int_t relid[4] ;
92862013 221 AliPHOSDigit * digit ;
46b146ca 222 AliPHOSEmcRecPoint * emc ;
223 AliPHOSPpsdRecPoint * ppsd ;
38f8ae6d 224 // AliPHOSTrackSegment * tracksegment ;
225 AliPHOSRecParticle * recparticle;
46b146ca 226
92862013 227 for ( ievent=0; ievent<Nevents; ievent++)
228 {
de9ec31b 229 //========== Event Number>
46b146ca 230 if ( ( log10((Float_t)(ievent+1)) - (Int_t)(log10((Float_t)(ievent+1))) ) == 0. )
231 cout << "AnalyzeManyEvents > " << "Event is " << ievent << endl ;
232
92862013 233 //=========== Connects the various Tree's for evt
234 gAlice->GetEvent(ievent);
46b146ca 235
92862013 236 //=========== Gets the Digit TTree
237 gAlice->TreeD()->GetEvent(0) ;
46b146ca 238
92862013 239 //=========== Gets the number of entries in the Digits array
240 TIter nextdigit(fPHOS->Digits()) ;
241 while( ( digit = (AliPHOSDigit *)nextdigit() ) )
242 {
243 fGeom->AbsToRelNumbering(digit->GetId(), relid) ;
46b146ca 244 if (fClu->IsInEmc(digit)) fhEmcDigit->Fill(fClu->Calibrate(digit->GetAmp())) ;
245 else
c3b9b3f9 246 {
247 if (relid[1]<17) fhVetoDigit->Fill(fClu->Calibrate(digit->GetAmp()));
248 if (relid[1]>16) fhConvertorDigit->Fill(fClu->Calibrate(digit->GetAmp()));
249 }
92862013 250 }
11e586af 251
46b146ca 252
253 //=========== Cluster in module
254 TIter nextEmc(*fPHOS->EmcRecPoints() ) ;
255 while((emc = (AliPHOSEmcRecPoint *)nextEmc()))
256 {
257 if ( emc->GetPHOSMod() == module )
258 {
259 fhEmcCluster->Fill( emc->GetTotalEnergy() );
260 TIter nextPpsd( *fPHOS->PpsdRecPoints()) ;
261 while((ppsd = (AliPHOSPpsdRecPoint *)nextPpsd()))
262 {
263 if ( ppsd->GetPHOSMod() == module )
264 {
265 if (!ppsd->GetUp()) fhConvertorEmc->Fill(emc->GetTotalEnergy(),ppsd->GetTotalEnergy()) ;
266 }
267 }
268 }
269 }
270
271 //=========== Cluster in module PPSD Down
272 TIter nextPpsd(*fPHOS->PpsdRecPoints() ) ;
273 while((ppsd = (AliPHOSPpsdRecPoint *)nextPpsd()))
274 {
275 if ( ppsd->GetPHOSMod() == module )
276 {
277 if (!ppsd->GetUp()) fhConvertorCluster->Fill(ppsd->GetTotalEnergy()) ;
278 if (ppsd->GetUp()) fhVetoCluster ->Fill(ppsd->GetTotalEnergy()) ;
279 }
280 }
281
92862013 282 //========== TRackSegments in the event
fc879520 283 TIter nextRecParticle(*fPHOS->RecParticles() ) ;
38f8ae6d 284 while((recparticle = (AliPHOSRecParticle *)nextRecParticle()))
92862013 285 {
38f8ae6d 286 if ( recparticle->GetPHOSTrackSegment()->GetPHOSMod() == module )
92862013 287 {
b2a60966 288 cout << "Particle type is " << recparticle->GetType() << endl ;
289 Int_t numberofprimaries = 0 ;
290 Int_t * listofprimaries = recparticle->GetPrimaries(numberofprimaries) ;
291 cout << "Number of primaries = " << numberofprimaries << endl ;
292 Int_t index ;
293 for ( index = 0 ; index < numberofprimaries ; index++)
294 cout << " primary # " << index << " = " << listofprimaries[index] << endl ;
134ce69a 295 }
296 }
134ce69a 297 } // endfor
46b146ca 298 SaveHistograms();
134ce69a 299 } // endif
300} // endfunction
301
302//____________________________________________________________________________
c3b9b3f9 303void AliPHOSAnalyze::ReconstructCPV(Int_t Nevents )
304{
305
306 // Perform reconstruction of EMC and CPV (GPS2 or IHEP) for <Nevents> events
307 // Yuri Kharlov. 19 October 2000
308
309 Int_t ievent ;
310 for ( ievent=0; ievent<Nevents; ievent++) {
311 if (ievent==0) {
312 cout << "Analyze > Starting Reconstructing " << endl ;
313 //========== Create the Clusterizer
314 fClu = new AliPHOSClusterizerv1() ;
315 fClu->SetEmcEnergyThreshold(0.05) ;
316 fClu->SetEmcClusteringThreshold(0.20) ;
317 fClu->SetLocalMaxCut(0.03) ;
318 if (strcmp(fGeom->GetName(),"GPS2") == 0) {
319 fClu->SetPpsdEnergyThreshold (0.0000002) ;
320 fClu->SetPpsdClusteringThreshold(0.0000001) ;
321 }
322 else if (strcmp(fGeom->GetName(),"IHEP") == 0) {
323 fClu->SetLocalMaxCutCPV(0.03) ;
324 fClu->SetLogWeightCutCPV(4.0) ;
325 fClu->SetPpsdEnergyThreshold (0.09) ;
326 }
327 fClu->SetCalibrationParameters(0., 0.00000001) ;
328
329 //========== Creates the track segment maker
330 fTrs = new AliPHOSTrackSegmentMakerv1() ;
331
332 //========== Creates the particle identifier for GPS2 only
333 if (strcmp(fGeom->GetName(),"GPS2") == 0) {
334 fPID = new AliPHOSPIDv1() ;
335 fPID->SetShowerProfileCuts(0.3, 1.8, 0.3, 1.8 ) ;
336 }
337
338 //========== Creates the Reconstructioner
339 fRec = new AliPHOSReconstructioner(fClu, fTrs, fPID) ;
340 if (fDebugLevel != 0) fRec -> SetDebugReconstruction(kTRUE);
341 }
342
343 if (fDebugLevel != 0 ||
344 (ievent+1) % (Int_t)TMath::Power( 10, (Int_t)TMath::Log10(ievent+1) ) == 0)
345 cout << "======= Analyze ======> Event " << ievent+1 << endl ;
346
347 //=========== Connects the various Tree's for evt
348 gAlice->GetEvent(ievent);
349
350 //=========== Gets the Digit TTree
351 gAlice->TreeD()->GetEvent(0) ;
352
353 //=========== Do the reconstruction
354 fPHOS->Reconstruction(fRec);
355 }
356
357 if(fClu) {delete fClu ; fClu =0 ;}
358 if(fPID) {delete fPID ; fPID =0 ;}
359 if(fRec) {delete fRec ; fRec =0 ;}
360 if(fTrs) {delete fTrs ; fTrs =0 ;}
361
362}
363//-------------------------------------------------------------------------------------
364
365void AliPHOSAnalyze::Reconstruct(Int_t Nevents )
46b146ca 366{
fc879520 367 Int_t ievent ;
368 for ( ievent=0; ievent<Nevents; ievent++)
369 {
370 if (ievent==0)
371 {
46b146ca 372 cout << "Analyze > Starting Reconstructing " << endl ;
fc879520 373 //========== Create the Clusterizer
374 fClu = new AliPHOSClusterizerv1() ;
375 fClu->SetEmcEnergyThreshold(0.05) ;
376 fClu->SetEmcClusteringThreshold(0.20) ;
377 fClu->SetPpsdEnergyThreshold (0.0000002) ;
378 fClu->SetPpsdClusteringThreshold(0.0000001) ;
379 fClu->SetLocalMaxCut(0.03) ;
380 fClu->SetCalibrationParameters(0., 0.00000001) ;
381
382 //========== Creates the track segment maker
383 fTrs = new AliPHOSTrackSegmentMakerv1() ;
384 // fTrs->UnsetUnfoldFlag() ;
385
386 //========== Creates the particle identifier
387 fPID = new AliPHOSPIDv1() ;
388 fPID->SetShowerProfileCuts(0.3, 1.8, 0.3, 1.8 ) ;
389
390 //========== Creates the Reconstructioner
391 fRec = new AliPHOSReconstructioner(fClu, fTrs, fPID) ;
c3b9b3f9 392// fRec -> SetDebugReconstruction(kTRUE);
134ce69a 393
fc879520 394 }
395
396 //========== Event Number>
c3b9b3f9 397 if ((ievent+1) % (Int_t)TMath::Power( 10, (Int_t)TMath::Log10(ievent+1) ) == 0)
398 cout << "======= Analyze ======> Event " << ievent+1 << endl ;
fc879520 399
400 //=========== Connects the various Tree's for evt
401 gAlice->GetEvent(ievent);
134ce69a 402
fc879520 403 //=========== Gets the Digit TTree
404 gAlice->TreeD()->GetEvent(0) ;
134ce69a 405
fc879520 406 //=========== Do the reconstruction
407 fPHOS->Reconstruction(fRec);
408 }
134ce69a 409
c3b9b3f9 410 if(fClu) {delete fClu ; fClu =0 ;}
411 if(fPID) {delete fPID ; fPID =0 ;}
412 if(fRec) {delete fRec ; fRec =0 ;}
413 if(fTrs) {delete fTrs ; fTrs =0 ;}
414
fc879520 415}
416//-------------------------------------------------------------------------------------
417
418// TClonesArray AllDigitArray = TClonesArray("AliPHOSDigit",1000) ;
419// TClonesArray * PhotonsList ;
420// TClonesArray * FalsDigitsList ;
421// TClonesArray AllPrimary = TClonesArray("TParticle",5000) ;
422// TFile * file2 = new TFile("ph100.root") ; // file with added photons
423// gAlice = (AliRun*) file2->Get("gAlice") ;
424// Int_t ievent;
425// Int_t NDigits[Nevents+1] ;
426// NDigits[0]=0 ;
427// Int_t NAllDigits = 0;
428// Int_t NprimPerEvent = 20 ;
429// for (ievent=0; ievent <Nevents; ievent++)
430// {
431// PhotonsList = gAlice->Particles(); //Primary
432// FalsDigitsList = ((AliPHOSv1 *)gAlice->GetDetector("PHOS"))->Digits(); //Digits
433// gAlice->GetEvent(ievent) ;
434// gAlice->TreeD()->GetEvent(0) ;
435// gAlice->TreeK()->GetEvent(0) ;
436// //Copy Primary
437// Int_t Nprim ;
438// for(Nprim = 0 ;Nprim < NprimPerEvent ; Nprim++)
439// new (AllPrimary[Nprim+ievent*NprimPerEvent]) TParticle(*((TParticle *) PhotonsList->At(Nprim))) ;
440
441// //Copy Digits
442// TIter nextDigit(FalsDigitsList) ;
443// AliPHOSDigit * FalseDigit ;
444// NDigits[ievent+1] = NDigits[ievent]+ FalsDigitsList->GetEntriesFast() ;
445// while( (FalseDigit = (AliPHOSDigit *) nextDigit()))
446// {
447// new (AllDigitArray[NAllDigits]) AliPHOSDigit(FalseDigit->GetPrimary(1),FalseDigit->GetId(),FalseDigit->GetAmp()) ;
448// NAllDigits++ ;
449// }
450// }
451// file2->Close() ;
452
453
454
455// //Add primary particles
456// cout << "# of Primaries before add " << PrimaryList->GetEntriesFast() << endl;
457// Int_t NTruePrimary = 0 ; //PrimaryList->GetEntriesFast() ;
458// Int_t Nprim ;
459// for(Nprim = 0; Nprim < NprimPerEvent; Nprim++)
460// new ((*PrimaryList)[NTruePrimary+Nprim]) TParticle(*((TParticle *) AllPrimary.At(Nprim+ievent*NprimPerEvent))) ;
461
462// cout << "# of Primaries after add " << PrimaryList->GetEntriesFast() <<endl;
463
464// cout << "Digits before add " << DigitsList->GetEntries() << endl ;
465// cout << "Digits to add " << NDigits[ievent+1]- NDigits[ievent]<< endl ;
466
467 //=========== Add fals digits ==============================
468// TIter nextDigit(DigitsList) ;
469// AliPHOSDigit * FalseDigit ;
470// AliPHOSDigit * RealDigit ;
471// Int_t NTrueDigits = DigitsList->GetEntriesFast() ;
472// Int_t Ndigit ;
473// for(Ndigit=NDigits[ievent];Ndigit<NDigits[ievent+1];Ndigit++)
474// {
475// FalseDigit = (AliPHOSDigit*) AllDigitArray.At(Ndigit) ;
476// Bool_t Add = kTRUE ;
477// AliPHOSDigit tmpDigit=AliPHOSDigit(FalseDigit->GetPrimary(1)+NTruePrimary,FalseDigit->GetId(),FalseDigit->GetAmp()) ;
478
479// while( (RealDigit = (AliPHOSDigit *) nextDigit()) && Add)
480// {
481// if((*RealDigit) == (tmpDigit))
482// {
483// *RealDigit=*RealDigit+tmpDigit ;
484// Add = kFALSE ;
485// }
486// }
487// if(Add)
488// {
489// new ((*DigitsList)[NTrueDigits]) AliPHOSDigit(FalseDigit->GetPrimary(1)+NTruePrimary,FalseDigit->GetId(),FalseDigit->GetAmp()) ;
490// ((AliPHOSDigit *)DigitsList->At(NTrueDigits))->SetIndexInList(NTrueDigits) ;
491// NTrueDigits++ ;
492// }
493// }
494// cout << "Digits after add " << DigitsList->GetEntries() << endl ;
134ce69a 495
134ce69a 496
c3b9b3f9 497//____________________________________________________________________________
498void AliPHOSAnalyze::ReadAndPrintCPV(Int_t Nevents)
499{
500 //
501 // Read and print generated and reconstructed hits in CPV
502 // Author: Yuri Kharlov
503 // 12 October 2000
504 //
505
506 cout << "Start CPV Analysis"<< endl ;
507 for ( Int_t ievent=0; ievent<Nevents; ievent++) {
508
509 //========== Event Number>
510 cout << endl << "==== ReadAndPrintCPV ====> Event is " << ievent+1 << endl ;
511
512 //=========== Connects the various Tree's for evt
513 gAlice->GetEvent(ievent);
514
515 //=========== Get the Hits Tree
516 gAlice->ResetHits();
517 gAlice->TreeH()->GetEvent(0);
518
519 //========== Creating branches ===================================
520 AliPHOSRecPoint::RecPointsList ** EmcRecPoints = fPHOS->EmcRecPoints() ;
521 gAlice->TreeR()->SetBranchAddress( "PHOSEmcRP" , EmcRecPoints ) ;
522
523 AliPHOSRecPoint::RecPointsList ** CpvRecPoints = fPHOS->PpsdRecPoints() ;
524 gAlice->TreeR()->SetBranchAddress( "PHOSPpsdRP", CpvRecPoints ) ;
525
526 //=========== Gets the Reconstruction TTree
527 gAlice->TreeR()->GetEvent(0) ;
528
529 // Read and print CPV hits
530
531 TClonesArray *CPVhits;
532 for (Int_t iModule=0; iModule < fGeom->GetNModules(); iModule++) {
533 CPVModule cpvModule = fPHOS->GetCPVModule(iModule);
534 CPVhits = cpvModule.Hits();
535 Int_t nCPVhits = CPVhits->GetEntriesFast();
536 for (Int_t ihit=0; ihit<nCPVhits; ihit++) {
537 CPVHit *cpvHit = (CPVHit*)CPVhits->UncheckedAt(ihit);
538 TLorentzVector p = cpvHit->GetMomentum();
539 Float_t xgen = cpvHit->GetX();
540 Float_t zgen = cpvHit->GetY();
541 Int_t ipart = cpvHit->GetIpart();
542 printf("CPV hit in module %d: ",iModule+1);
543 printf(" p = (%f, %f, %f, %f) GeV,\n",
544 p.Px(),p.Py(),p.Pz(),p.Energy());
545 printf(" xy = (%8.4f, %8.4f) cm, ipart = %d\n",
546 xgen,zgen,ipart);
547 }
548 }
549
550 // Read and print CPV reconstructed points
551
552 TIter nextRP(*fPHOS->PpsdRecPoints() ) ;
553 AliPHOSPpsdRecPoint *cpvRecPoint ;
554 while( ( cpvRecPoint = (AliPHOSPpsdRecPoint *)nextRP() ) ) {
555 TVector3 locpos;
556 cpvRecPoint->GetLocalPosition(locpos);
557 Int_t PHOSModule = cpvRecPoint->GetPHOSMod();
558 printf("CPV recpoint in module %d: (X,Y,Z) = (%f,%f,%f) cm\n",
559 PHOSModule,locpos.X(),locpos.Y(),locpos.Z());
560 }
561 }
562}
563
564//____________________________________________________________________________
565void AliPHOSAnalyze::AnalyzeCPV(Int_t Nevents)
566{
567 //
568 // Analyzes CPV characteristics
569 // Author: Yuri Kharlov
570 // 9 October 2000
571 //
572
573 // Book histograms
574
575 TH1F *hDx = new TH1F("hDx" ,"CPV x-resolution@reconstruction",100,-5. , 5.);
576 TH1F *hDz = new TH1F("hDz" ,"CPV z-resolution@reconstruction",100,-5. , 5.);
577 TH1S *hNrp = new TH1S("hNrp" ,"CPV rec.point multiplicity", 21,-0.5,20.5);
578
579 cout << "Start CPV Analysis"<< endl ;
580 for ( Int_t ievent=0; ievent<Nevents; ievent++) {
581
582 //========== Event Number>
583 if ( (ievent+1) % (Int_t)TMath::Power( 10, (Int_t)TMath::Log10(ievent+1) ) == 0)
584 cout << endl << "==== AnalyzeCPV ====> Event is " << ievent+1 << endl ;
585
586 //=========== Connects the various Tree's for evt
587 gAlice->GetEvent(ievent);
588
589 //=========== Get the Hits Tree
590 gAlice->ResetHits();
591 gAlice->TreeH()->GetEvent(0);
592
593 //========== Creating branches ===================================
594 AliPHOSRecPoint::RecPointsList ** EmcRecPoints = fPHOS->EmcRecPoints() ;
595 gAlice->TreeR()->SetBranchAddress( "PHOSEmcRP" , EmcRecPoints ) ;
596
597 AliPHOSRecPoint::RecPointsList ** CpvRecPoints = fPHOS->PpsdRecPoints() ;
598 gAlice->TreeR()->SetBranchAddress( "PHOSPpsdRP", CpvRecPoints ) ;
599
600 //=========== Gets the Reconstruction TTree
601 gAlice->TreeR()->GetEvent(0) ;
602
603 TIter nextRP(*fPHOS->PpsdRecPoints() ) ;
604 AliPHOSEmcRecPoint *cpvRecPoint ;
605 CPVModule cpvModule;
606 TClonesArray *CPVhits;
607 while( ( cpvRecPoint = (AliPHOSEmcRecPoint *)nextRP() ) ) {
608 TVector3 locpos;
609 cpvRecPoint->GetLocalPosition(locpos);
610 Int_t PHOSModule = cpvRecPoint->GetPHOSMod();
611 Int_t rpMult = cpvRecPoint->GetDigitsMultiplicity();
612 Float_t xrec = locpos.X();
613 Float_t zrec = locpos.Z();
614 Float_t dxmin = 1.e+10;
615 Float_t dzmin = 1.e+10;
616
617 cpvModule = fPHOS->GetCPVModule(PHOSModule-1);
618 CPVhits = cpvModule.Hits();
619 Int_t nCPVhits = CPVhits->GetEntriesFast();
620 for (Int_t ihit=0; ihit<nCPVhits; ihit++) {
621 CPVHit *cpvHit = (CPVHit*)CPVhits->UncheckedAt(ihit);
622 Float_t xgen = cpvHit->GetX();
623 Float_t zgen = cpvHit->GetY();
624 if ( TMath::Abs(xgen-xrec) < TMath::Abs(dxmin) ) dxmin = xgen-xrec;
625 if ( TMath::Abs(zgen-zrec) < TMath::Abs(dzmin) ) dzmin = zgen-zrec;
626 }
627 cpvModule.Clear();
628 hDx ->Fill(dxmin);
629 hDz ->Fill(dzmin);
630 hNrp ->Fill(rpMult);
631 }
632
633 }
634
635 // Save histograms
636
637 Text_t outputname[80] ;
638 sprintf(outputname,"%s.analyzed",fRootFile->GetName());
639 TFile output(outputname,"RECREATE");
640 output.cd();
641
642 hDx ->Write() ;
643 hDz ->Write() ;
644 hNrp ->Write() ;
645
646 // Plot histograms
647
648 TCanvas *CPVcanvas = new TCanvas("CPV","CPV analysis",20,20,300,900);
649 gStyle->SetOptStat(111111);
650 gStyle->SetOptFit(1);
651 gStyle->SetOptDate(1);
652 CPVcanvas->Divide(3,1);
653
654 CPVcanvas->cd(1);
655 gPad->SetFillColor(10);
656 hNrp->SetFillColor(16);
657 hNrp->Draw();
658
659 CPVcanvas->cd(2);
660 gPad->SetFillColor(10);
661 hDx->SetFillColor(16);
662 hDx->Fit("gaus");
663 hDx->Draw();
664
665 CPVcanvas->cd(3);
666 gPad->SetFillColor(10);
667 hDz->SetFillColor(16);
668 hDz->Fit("gaus");
669 hDz->Draw();
670
671 CPVcanvas->Print("CPV.ps");
672
673}
674
fc879520 675//____________________________________________________________________________
676 void AliPHOSAnalyze::AnalyzeResolutions(Int_t Nevents )
677{
678 // analyzes Nevents events and calculate Energy and Position resolution as well as
679 // probaility of correct indentifiing of the incident particle
134ce69a 680
fc879520 681 //========== Booking Histograms
682 cout << "AnalyzeResolutions > " << "Booking Histograms" << endl ;
683 BookResolutionHistograms();
134ce69a 684
fc879520 685 Int_t Counter[9][5] ;
686 Int_t i1,i2,TotalInd = 0 ;
687 for(i1 = 0; i1<9; i1++)
688 for(i2 = 0; i2<5; i2++)
689 Counter[i1][i2] = 0 ;
690
691 Int_t TotalPrimary = 0 ;
692 Int_t TotalRecPart = 0 ;
693 Int_t TotalRPwithPrim = 0 ;
694 Int_t ievent;
695
696 cout << "Start Analysing"<< endl ;
697 for ( ievent=0; ievent<Nevents; ievent++)
698 {
699
700 //========== Event Number>
701 if ( ( log10((Float_t)(ievent+1)) - (Int_t)(log10((Float_t)(ievent+1))) ) == 0. )
702 cout << "AnalyzeResolutions > " << "Event is " << ievent << endl ;
703
704 //=========== Connects the various Tree's for evt
705 gAlice->GetEvent(ievent);
134ce69a 706
fc879520 707
134ce69a 708
fc879520 709 //=========== Gets the Kine TTree
710 gAlice->TreeK()->GetEvent(0) ;
711
712 //=========== Gets the list of Primari Particles
713 TClonesArray * PrimaryList = gAlice->Particles();
714
715 TParticle * Primary ;
716 Int_t iPrimary ;
717 for ( iPrimary = 0 ; iPrimary < PrimaryList->GetEntries() ; iPrimary++)
134ce69a 718 {
fc879520 719 Primary = (TParticle*)PrimaryList->UncheckedAt(iPrimary) ;
720 Int_t PrimaryType = Primary->GetPdgCode() ;
721 if( PrimaryType == 22 )
722 fhPrimary->Fill(Primary->Energy()) ;
723 }
724
725 //=========== Get the Digit Tree
726 gAlice->TreeD()->GetEvent(0) ;
134ce69a 727
fc879520 728 //========== Creating branches ===================================
729 AliPHOSRecPoint::RecPointsList ** EmcRecPoints = fPHOS->EmcRecPoints() ;
730 gAlice->TreeR()->SetBranchAddress( "PHOSEmcRP", EmcRecPoints ) ;
134ce69a 731
fc879520 732 AliPHOSRecPoint::RecPointsList ** PpsdRecPoints = fPHOS->PpsdRecPoints() ;
733 gAlice->TreeR()->SetBranchAddress( "PHOSPpsdRP", PpsdRecPoints ) ;
734
735 AliPHOSTrackSegment::TrackSegmentsList ** TrackSegmentsList = fPHOS->TrackSegments() ;
736 if( (*TrackSegmentsList) )
737 (*TrackSegmentsList)->Clear() ;
738 gAlice->TreeR()->SetBranchAddress( "PHOSTS", TrackSegmentsList ) ;
739
740 AliPHOSRecParticle::RecParticlesList ** RecParticleList = fPHOS->RecParticles() ;
741 if( (*RecParticleList) )
742 (*RecParticleList)->Clear() ;
743 gAlice->TreeR()->SetBranchAddress( "PHOSRP", RecParticleList ) ;
744
134ce69a 745
fc879520 746 //=========== Gets the Reconstraction TTree
747 gAlice->TreeR()->GetEvent(0) ;
748
46b146ca 749 cout << ievent << " " << (*EmcRecPoints) << " " <<(*PpsdRecPoints) <<fPHOS->Digits()<< endl ;
750 cout << " " << " " << (*EmcRecPoints)->GetEntries() << " " <<(*PpsdRecPoints)->GetEntries() <<fPHOS->Digits()->GetEntries()<< endl ;
751
fc879520 752 AliPHOSRecParticle * RecParticle ;
753 Int_t iRecParticle ;
754 for(iRecParticle = 0; iRecParticle < (*RecParticleList)->GetEntries() ;iRecParticle++ )
755 {
756 RecParticle = (AliPHOSRecParticle *) (*RecParticleList)->At(iRecParticle) ;
757
758 Int_t ModuleNumberRec ;
759 Double_t RecX, RecZ ;
760 fGeom->ImpactOnEmc(RecParticle->Theta(), RecParticle->Phi(), ModuleNumberRec, RecX, RecZ) ;
761
762 Double_t MinDistance = 10000 ;
763 Int_t ClosestPrimary = -1 ;
764
765 Int_t numberofprimaries ;
766 Int_t * listofprimaries = RecParticle->GetPrimaries(numberofprimaries) ;
767 Int_t index ;
768 TParticle * Primary ;
769 Double_t Distance = MinDistance ;
770 for ( index = 0 ; index < numberofprimaries ; index++)
771 {
772 Primary = (TParticle*)PrimaryList->UncheckedAt(listofprimaries[index]) ;
773 Int_t ModuleNumber ;
774 Double_t PrimX, PrimZ ;
775 fGeom->ImpactOnEmc(Primary->Theta(), Primary->Phi(), ModuleNumber, PrimX, PrimZ) ;
776 if(ModuleNumberRec == ModuleNumber)
777 Distance = TMath::Sqrt((RecX-PrimX)*(RecX-PrimX)+(RecZ-PrimZ)*(RecZ-PrimZ) ) ;
778 if(MinDistance > Distance)
134ce69a 779 {
fc879520 780 MinDistance = Distance ;
781 ClosestPrimary = listofprimaries[index] ;
782 }
783 }
784 TotalRecPart++ ;
785 if(ClosestPrimary >=0 )
786 {
787 fhPhotonAllEnergy->Fill(((TParticle *) PrimaryList->At(ClosestPrimary))->Energy(), RecParticle->Energy() ) ;
788 fhPhotonAllPosition->Fill(((TParticle *) PrimaryList->At(ClosestPrimary))->Energy(),Distance) ;
789 TotalRPwithPrim++;
790 Int_t PrimaryType = ((TParticle *)PrimaryList->At(ClosestPrimary))->GetPdgCode() ;
791 TParticlePDG* PDGparticle = ((TParticle *)PrimaryList->At(ClosestPrimary))->GetPDG();
792 Double_t charge = PDGparticle->Charge() ;
793 Int_t PrimaryCode ;
794 switch(PrimaryType)
795 {
796 case 22:
797 PrimaryCode = 0; //Photon
798 break;
799 case 11 :
800 PrimaryCode = 1; //Electron
801 break;
802 case -11 :
803 PrimaryCode = 1; //positron
804 break;
805 case 321 :
806 PrimaryCode = 4; //K+
807 break;
808 case -321 :
809 PrimaryCode = 4; //K-
810 break;
811 case 310 :
812 PrimaryCode = 4; //K0s
813 break;
814 case 130 :
815 PrimaryCode = 4; //K0l
816 break;
817 default:
818 if(charge)
819 PrimaryCode = 2; //Charged hadron
820 else
821 PrimaryCode = 3; //Neutral hadron
822 break;
134ce69a 823 }
824
fc879520 825 switch(RecParticle->GetType())
134ce69a 826 {
fc879520 827 case AliPHOSFastRecParticle::kGAMMA:
828 if(PrimaryType == 22){
829 fhPhotonEnergy->Fill(((TParticle *) PrimaryList->At(ClosestPrimary))->Energy(), RecParticle->Energy() ) ;
830 fhPhotonPosition->Fill(((TParticle *) PrimaryList->At(ClosestPrimary))->Energy(),Distance) ;
831 fhPhotonReg->Fill(RecParticle->Energy() ) ;
832 fhPhotonEM->Fill(RecParticle->Energy() ) ;
833 fhPhotPhot->Fill(RecParticle->Energy() ) ;
834 }
835 if(PrimaryType == 2112){ //neutron
836 fhNReg->Fill(RecParticle->Energy() ) ;
837 fhNEM->Fill(RecParticle->Energy() ) ;
838 }
839
840 if(PrimaryType == -2112){ //neutron ~
841 fhNBarReg->Fill(RecParticle->Energy() ) ;
842 fhNBarEM->Fill(RecParticle->Energy() ) ;
843
844 }
845 if(PrimaryCode == 2){
846 fhChargedReg->Fill(RecParticle->Energy() ) ;
847 fhChargedEM->Fill(RecParticle->Energy() ) ;
848 }
849
850 fhAllReg->Fill(RecParticle->Energy() ) ;
851 fhAllEM->Fill(RecParticle->Energy() ) ;
852 Counter[0][PrimaryCode]++;
853 break;
854 case AliPHOSFastRecParticle::kELECTRON:
855 if(PrimaryType == 11 || PrimaryType == -11){
856 fhElectronEnergy->Fill(((TParticle *)PrimaryList->At(ClosestPrimary))->Energy(), RecParticle->Energy() ) ;
857 fhElectronPosition->Fill(((TParticle *)PrimaryList->At(ClosestPrimary))->Energy(),Distance ) ;
858 }
859 if(PrimaryType == 22)
860 fhPhotElec->Fill(RecParticle->Energy() ) ;
861
862 Counter[1][PrimaryCode]++;
863 break;
864 case AliPHOSFastRecParticle::kNEUTRALHA:
865 if(PrimaryType == 22)
866 fhPhotNeuH->Fill(RecParticle->Energy() ) ;
867
868 fhNeutralHadronEnergy->Fill( ((TParticle *)PrimaryList->At(ClosestPrimary))->Energy(), RecParticle->Energy()) ;
869 fhNeutralHadronPosition->Fill(((TParticle *)PrimaryList->At(ClosestPrimary))->Energy() ,Distance ) ;
870 Counter[2][PrimaryCode]++;
871 break ;
872 case AliPHOSFastRecParticle::kNEUTRALEM:
873 if(PrimaryType == 22 || PrimaryType == 11 || PrimaryType == -11){
874 fhNeutralEMEnergy->Fill(((TParticle *)PrimaryList->At(ClosestPrimary))->Energy(),RecParticle->Energy() ) ;
875 fhNeutralEMPosition->Fill(((TParticle *)PrimaryList->At(ClosestPrimary))->Energy(),Distance ) ;
876 }
877
878 if(PrimaryType == 22){ //photon
879 fhPhotNuEM->Fill(RecParticle->Energy() ) ;
880 fhPhotonEM->Fill(RecParticle->Energy() ) ;
881 }
882 if(PrimaryType == 2112) //neutron
883 fhNEM->Fill(RecParticle->Energy() ) ;
884
885 if(PrimaryType == -2112) //neutron ~
886 fhNBarEM->Fill(RecParticle->Energy() ) ;
887
888 if(PrimaryCode == 2)
889 fhChargedEM->Fill(RecParticle->Energy() ) ;
890
891 fhAllEM->Fill(RecParticle->Energy() ) ;
892
893 Counter[3][PrimaryCode]++;
894 break ;
895 case AliPHOSFastRecParticle::kCHARGEDHA:
896 if(PrimaryType == 22) //photon
897 fhPhotChHa->Fill(RecParticle->Energy() ) ;
898
899 fhChargedHadronEnergy->Fill(((TParticle *)PrimaryList->At(ClosestPrimary))->Energy(),RecParticle->Energy() ) ;
900 fhChargedHadronPosition->Fill(((TParticle *)PrimaryList->At(ClosestPrimary))->Energy(),Distance ) ;
901 Counter[4][PrimaryCode]++ ;
902 break ;
903 case AliPHOSFastRecParticle::kGAMMAHA:
904 if(PrimaryType == 22) //photon
905 fhPhotGaHa->Fill(RecParticle->Energy() ) ;
906 fhPhotonHadronEnergy->Fill(((TParticle *)PrimaryList->At(ClosestPrimary))->Energy(), RecParticle->Energy()) ;
907 fhPhotonHadronPosition->Fill(((TParticle *)PrimaryList->At(ClosestPrimary))->Energy(),Distance ) ;
908 Counter[5][PrimaryCode]++ ;
909 break ;
910 case AliPHOSFastRecParticle::kABSURDEM:
911 Counter[6][PrimaryCode]++ ;
912 break;
913 case AliPHOSFastRecParticle::kABSURDHA:
914 Counter[7][PrimaryCode]++ ;
915 break;
916 default:
917 Counter[8][PrimaryCode]++ ;
918 break;
92862013 919 }
920 }
fc879520 921 }
922 } // endfor
46b146ca 923 SaveHistograms();
fc879520 924 cout << "Resolutions: Analyzed " << Nevents << " event(s)" << endl ;
925 cout << "Resolutions: Total primary " << TotalPrimary << endl ;
926 cout << "Resoluitons: Total reconstracted " << TotalRecPart << endl ;
927 cout << "TotalReconstructed with Primarie " << TotalRPwithPrim << endl ;
928 cout << " Primary: Photon Electron Ch. Hadr. Neutr. Hadr Kaons" << endl ;
929 cout << " Detected as photon " << Counter[0][0] << " " << Counter[0][1] << " " << Counter[0][2] << " " <<Counter[0][3] << " " << Counter[0][4] << endl ;
930 cout << " Detected as electron " << Counter[1][0] << " " << Counter[1][1] << " " << Counter[1][2] << " " <<Counter[1][3] << " " << Counter[1][4] << endl ;
931 cout << " Detected as neutral hadron " << Counter[2][0] << " " << Counter[2][1] << " " << Counter[2][2] << " " <<Counter[2][3] << " " << Counter[2][4] << endl ;
932 cout << " Detected as neutral EM " << Counter[3][0] << " " << Counter[3][1] << " " << Counter[3][2] << " " <<Counter[3][3] << " " << Counter[3][4] << endl ;
933 cout << " Detected as charged hadron " << Counter[4][0] << " " << Counter[4][1] << " " << Counter[4][2] << " " <<Counter[4][3] << " " << Counter[4][4] << endl ;
934 cout << " Detected as gamma-hadron " << Counter[5][0] << " " << Counter[5][1] << " " << Counter[5][2] << " " <<Counter[5][3] << " " << Counter[5][4] << endl ;
935 cout << " Detected as Absurd EM " << Counter[6][0] << " " << Counter[6][1] << " " << Counter[6][2] << " " <<Counter[6][3] << " " << Counter[6][4] << endl ;
936 cout << " Detected as absurd hadron " << Counter[7][0] << " " << Counter[7][1] << " " << Counter[7][2] << " " <<Counter[7][3] << " " << Counter[7][4] << endl ;
937 cout << " Detected as undefined " << Counter[8][0] << " " << Counter[8][1] << " " << Counter[8][2] << " " <<Counter[8][3] << " " << Counter[8][4] << endl ;
938
939 for(i1 = 0; i1<9; i1++)
940 for(i2 = 0; i2<5; i2++)
941 TotalInd+=Counter[i1][i2] ;
942 cout << "Indentified particles " << TotalInd << endl ;
943
92862013 944} // endfunction
945
946
947//____________________________________________________________________________
948void AliPHOSAnalyze::BookingHistograms()
949{
b2a60966 950 // Books the histograms where the results of the analysis are stored (to be changed)
951
eecb6765 952 delete fhEmcDigit ;
953 delete fhVetoDigit ;
954 delete fhConvertorDigit ;
955 delete fhEmcCluster ;
956 delete fhVetoCluster ;
957 delete fhConvertorCluster ;
958 delete fhConvertorEmc ;
959
46b146ca 960 fhEmcDigit = new TH1F("hEmcDigit", "hEmcDigit", 1000, 0. , 25.);
961 fhVetoDigit = new TH1F("hVetoDigit", "hVetoDigit", 500, 0. , 3.e-5);
962 fhConvertorDigit = new TH1F("hConvertorDigit","hConvertorDigit", 500, 0. , 3.e-5);
963 fhEmcCluster = new TH1F("hEmcCluster", "hEmcCluster", 1000, 0. , 30.);
964 fhVetoCluster = new TH1F("hVetoCluster", "hVetoCluster", 500, 0. , 3.e-5);
965 fhConvertorCluster = new TH1F("hConvertorCluster","hConvertorCluster",500, 0. , 3.e-5);
966 fhConvertorEmc = new TH2F("hConvertorEmc", "hConvertorEmc", 200, 1. , 3., 200, 0., 3.e-5);
92862013 967
134ce69a 968}
969//____________________________________________________________________________
970void AliPHOSAnalyze::BookResolutionHistograms()
971{
972 // Books the histograms where the results of the Resolution analysis are stored
973
974 if(fhPhotonEnergy)
975 delete fhPhotonEnergy ;
fc879520 976 if(fhPhotonAllEnergy)
977 delete fhPhotonAllEnergy ;
134ce69a 978 if(fhElectronEnergy)
979 delete fhElectronEnergy ;
fc879520 980 if(fhElectronAllEnergy)
981 delete fhElectronAllEnergy ;
134ce69a 982 if(fhNeutralHadronEnergy)
983 delete fhNeutralHadronEnergy ;
984 if(fhNeutralEMEnergy)
985 delete fhNeutralEMEnergy ;
fc879520 986 if(fhNeutralEMAllEnergy)
987 delete fhNeutralEMAllEnergy ;
134ce69a 988 if(fhChargedHadronEnergy)
989 delete fhChargedHadronEnergy ;
990 if(fhPhotonHadronEnergy)
991 delete fhPhotonHadronEnergy ;
992 if(fhPhotonPosition)
993 delete fhPhotonPosition ;
fc879520 994 if(fhPhotonAllPosition)
995 delete fhPhotonAllPosition ;
134ce69a 996 if(fhElectronPosition)
997 delete fhElectronPosition ;
fc879520 998 if(fhElectronAllPosition)
999 delete fhElectronAllPosition ;
134ce69a 1000 if(fhNeutralHadronPosition)
1001 delete fhNeutralHadronPosition ;
1002 if(fhNeutralEMPosition)
1003 delete fhNeutralEMPosition ;
fc879520 1004 if(fhNeutralEMAllPosition)
1005 delete fhNeutralEMAllPosition ;
134ce69a 1006 if(fhChargedHadronPosition)
1007 delete fhChargedHadronPosition ;
1008 if(fhPhotonHadronPosition)
1009 delete fhPhotonHadronPosition ;
1010
1011 fhPhotonEnergy = new TH2F("hPhotonEnergy", "hPhotonEnergy", 100, 0., 5., 100, 0., 5.);
fc879520 1012 fhPhotonAllEnergy = new TH2F("hPhotonAllEnergy", "hPhotonAllEnergy", 100, 0., 5., 100, 0., 5.);
134ce69a 1013 fhElectronEnergy = new TH2F("hElectronEnergy","hElectronEnergy", 100, 0., 5., 100, 0., 5.);
fc879520 1014 fhElectronAllEnergy = new TH2F("hElectronAllEnergy","hElectronAllEnergy", 100, 0., 5., 100, 0., 5.);
134ce69a 1015 fhNeutralHadronEnergy = new TH2F("hNeutralHadronEnergy", "hNeutralHadronEnergy", 100, 0., 5., 100, 0., 5.);
1016 fhNeutralEMEnergy = new TH2F("hNeutralEMEnergy", "hNeutralEMEnergy", 100, 0., 5., 100, 0., 5.);
fc879520 1017 fhNeutralEMAllEnergy = new TH2F("hNeutralEMAllEnergy", "hNeutralEMAllEnergy", 100, 0., 5., 100, 0., 5.);
134ce69a 1018 fhChargedHadronEnergy = new TH2F("hChargedHadronEnergy", "hChargedHadronEnergy", 100, 0., 5., 100, 0., 5.);
1019 fhPhotonHadronEnergy = new TH2F("hPhotonHadronEnergy","hPhotonHadronEnergy", 100, 0., 5., 100, 0., 5.);
fc879520 1020 fhPhotonPosition = new TH2F("hPhotonPosition","hPhotonPosition", 20, 0., 5., 100, 0., 5.);
1021 fhPhotonAllPosition = new TH2F("hPhotonAllPosition","hPhotonAllPosition", 20, 0., 5., 100, 0., 5.);
1022 fhElectronPosition = new TH2F("hElectronPosition","hElectronPosition", 20, 0., 5., 100, 0., 5.);
1023 fhElectronAllPosition = new TH2F("hElectronAllPosition","hElectronAllPosition", 20, 0., 5., 100, 0., 5.);
1024 fhNeutralHadronPosition = new TH2F("hNeutralHadronPosition","hNeutralHadronPosition", 20, 0., 5., 100, 0., 5.);
1025 fhNeutralEMPosition = new TH2F("hNeutralEMPosition","hNeutralEMPosition", 20, 0., 5., 100, 0., 5.);
1026 fhNeutralEMAllPosition = new TH2F("hNeutralEMAllPosition","hNeutralEMAllPosition", 20, 0., 5., 100, 0., 5.);
1027 fhChargedHadronPosition = new TH2F("hChargedHadronPosition","hChargedHadronPosition", 20, 0., 5., 100, 0., 5.);
1028 fhPhotonHadronPosition = new TH2F("hPhotonHadronPosition","hPhotonHadronPosition", 20, 0., 5., 100, 0., 5.);
46b146ca 1029
1030 if(fhPhotonReg)
1031 delete fhPhotonReg ;
1032 if(fhAllReg)
1033 delete fhAllReg ;
1034 if(fhNReg)
1035 delete fhNReg ;
1036 if(fhNReg)
1037 delete fhNReg ;
1038 if(fhNReg)
1039 delete fhNReg ;
1040
fc879520 1041 fhPhotonReg = new TH1F("hPhotonReg","hPhotonReg", 20, 0., 5.);
1042 fhAllReg = new TH1F("hAllReg", "hAllReg", 20, 0., 5.);
1043 fhNReg = new TH1F("hNReg", "hNReg", 20, 0., 5.);
1044 fhNBarReg = new TH1F("hNBarReg", "hNBarReg", 20, 0., 5.);
1045 fhChargedReg= new TH1F("hChargedReg", "hChargedReg", 20, 0., 5.);
46b146ca 1046
1047 if(fhPhotonEM)
1048 delete fhPhotonEM ;
1049 if(fhAllEM)
1050 delete fhAllEM ;
1051 if(fhNEM)
1052 delete fhNEM ;
1053 if(fhNBarEM)
1054 delete fhNBarEM ;
1055 if(fhChargedEM)
1056 delete fhChargedEM ;
1057
fc879520 1058 fhPhotonEM = new TH1F("hPhotonEM","hPhotonEM", 20, 0., 5.);
1059 fhAllEM = new TH1F("hAllEM", "hAllEM", 20, 0., 5.);
1060 fhNEM = new TH1F("hNEM", "hNEM", 20, 0., 5.);
1061 fhNBarEM = new TH1F("hNBarEM", "hNBarEM", 20, 0., 5.);
1062 fhChargedEM= new TH1F("hChargedEM", "hChargedEM", 20, 0., 5.);
46b146ca 1063
1064 if(fhPrimary)
1065 delete fhPrimary ;
fc879520 1066 fhPrimary= new TH1F("hPrimary", "hPrimary", 20, 0., 5.);
1067
46b146ca 1068 if(fhPhotPhot)
1069 delete fhPhotPhot ;
1070 if(fhPhotElec)
1071 delete fhPhotElec ;
1072 if(fhPhotNeuH)
1073 delete fhPhotNeuH ;
1074 if(fhPhotNuEM)
1075 delete fhPhotNuEM ;
1076 if(fhPhotChHa)
1077 delete fhPhotChHa ;
1078 if(fhPhotGaHa)
1079 delete fhPhotGaHa ;
1080
fc879520 1081 fhPhotPhot = new TH1F("hPhotPhot","hPhotPhot", 20, 0., 5.); //Photon registered as photon
1082 fhPhotElec = new TH1F("hPhotElec","hPhotElec", 20, 0., 5.); //Photon registered as Electron
1083 fhPhotNeuH = new TH1F("hPhotNeuH","hPhotNeuH", 20, 0., 5.); //Photon registered as Neutral Hadron
1084 fhPhotNuEM = new TH1F("hPhotNuEM","hPhotNuEM", 20, 0., 5.); //Photon registered as Neutral EM
1085 fhPhotChHa = new TH1F("hPhotChHa","hPhotChHa", 20, 0., 5.); //Photon registered as Charged Hadron
1086 fhPhotGaHa = new TH1F("hPhotGaHa","hPhotGaHa", 20, 0., 5.); //Photon registered as Gamma-Hadron
1087
3fc07a80 1088
92862013 1089}
6ad0bfa0 1090//____________________________________________________________________________
1091Bool_t AliPHOSAnalyze::Init(Int_t evt)
1092{
b2a60966 1093 // Do a few initializations: open the root file
1094 // get the AliRun object
1095 // defines the clusterizer, tracksegment maker and particle identifier
1096 // sets the associated parameters
6ad0bfa0 1097
92862013 1098 Bool_t ok = kTRUE ;
6ad0bfa0 1099
1100 //========== Open galice root file
1101
1102 if ( fRootFile == 0 ) {
1103 Text_t * name = new Text_t[80] ;
1104 cout << "AnalyzeOneEvent > Enter file root file name : " ;
1105 cin >> name ;
92862013 1106 Bool_t ok = OpenRootFile(name) ;
1107 if ( !ok )
6ad0bfa0 1108 cout << " AliPHOSAnalyze > Error opening " << name << endl ;
1109 else {
1110 //========== Get AliRun object from file
1111
1112 gAlice = (AliRun*) fRootFile->Get("gAlice") ;
1113
1114 //=========== Get the PHOS object and associated geometry from the file
1115
788dcca4 1116 fPHOS = (AliPHOSv1 *)gAlice->GetDetector("PHOS") ;
aafe457d 1117 fGeom = fPHOS->GetGeometry();
1118 // fGeom = AliPHOSGeometry::GetInstance( fPHOS->GetGeometry()->GetName(), fPHOS->GetGeometry()->GetTitle() );
83974468 1119
92862013 1120 } // else !ok
6ad0bfa0 1121 } // if fRootFile
1122
92862013 1123 if ( ok ) {
6ad0bfa0 1124
1125 //========== Create the Clusterizer
1126
774e78c2 1127 fClu = new AliPHOSClusterizerv1() ;
1128 fClu->SetEmcEnergyThreshold(0.030) ;
aafe457d 1129 fClu->SetEmcClusteringThreshold(0.20) ;
92862013 1130 fClu->SetPpsdEnergyThreshold (0.0000002) ;
1131 fClu->SetPpsdClusteringThreshold(0.0000001) ;
6ad0bfa0 1132 fClu->SetLocalMaxCut(0.03) ;
92862013 1133 fClu->SetCalibrationParameters(0., 0.00000001) ;
6ad0bfa0 1134 cout << "AnalyzeOneEvent > using clusterizer " << fClu->GetName() << endl ;
1135 fClu->PrintParameters() ;
1136
1137 //========== Creates the track segment maker
1138
1139 fTrs = new AliPHOSTrackSegmentMakerv1() ;
1140 cout << "AnalyzeOneEvent > using tack segment maker " << fTrs->GetName() << endl ;
134ce69a 1141 // fTrs->UnsetUnfoldFlag() ;
6ad0bfa0 1142
26d4b141 1143 //========== Creates the particle identifier
6ad0bfa0 1144
26d4b141 1145 fPID = new AliPHOSPIDv1() ;
1146 cout << "AnalyzeOneEvent > using particle identifier " << fPID->GetName() << endl ;
2aad621e 1147 //fPID->SetShowerProfileCuts(Float_t l1m, Float_t l1M, Float_t l2m, Float_t l2M) ;
1148 fPID->SetShowerProfileCuts(0.7, 2.0 , 0.6 , 1.5) ;
1149
6ad0bfa0 1150 //========== Creates the Reconstructioner
1151
2aad621e 1152 fRec = new AliPHOSReconstructioner(fClu, fTrs, fPID) ;
c3b9b3f9 1153// fRec -> SetDebugReconstruction(kFALSE);
1154 fRec -> SetDebugReconstruction(kTRUE);
6ad0bfa0 1155
1156 //=========== Connect the various Tree's for evt
1157
1158 if ( evt == -999 ) {
1159 cout << "AnalyzeOneEvent > Enter event number : " ;
1160 cin >> evt ;
1161 cout << evt << endl ;
1162 }
1163 fEvt = evt ;
1164 gAlice->GetEvent(evt);
1165
1166 //=========== Get the Digit TTree
1167
1168 gAlice->TreeD()->GetEvent(0) ;
1169
92862013 1170 } // ok
6ad0bfa0 1171
92862013 1172 return ok ;
6ad0bfa0 1173}
1174
92862013 1175
6ad0bfa0 1176//____________________________________________________________________________
92862013 1177void AliPHOSAnalyze::DisplayKineEvent(Int_t evt)
6ad0bfa0 1178{
b2a60966 1179 // Display particles from the Kine Tree in global Alice (theta, phi) coordinates.
1180 // One PHOS module at the time.
1181 // The particle type can be selected.
1182
6ad0bfa0 1183 if (evt == -999)
1184 evt = fEvt ;
1185
92862013 1186 Int_t module ;
6ad0bfa0 1187 cout << "DisplayKineEvent > which module (1-5, -1: all) ? " ;
92862013 1188 cin >> module ; cout << module << endl ;
6ad0bfa0 1189
92862013 1190 Int_t testparticle ;
6ad0bfa0 1191 cout << " 22 : PHOTON " << endl
1192 << " (-)11 : (POSITRON)ELECTRON " << endl
1193 << " (-)2112 : (ANTI)NEUTRON " << endl
1194 << " -999 : Everything else " << endl ;
1195 cout << "DisplayKineEvent > enter PDG particle code to display " ;
92862013 1196 cin >> testparticle ; cout << testparticle << endl ;
6ad0bfa0 1197
92862013 1198 Text_t histoname[80] ;
1199 sprintf(histoname,"Event %d: Incident particles in module %d", evt, module) ;
6ad0bfa0 1200
92862013 1201 Double_t tm, tM, pm, pM ; // min and Max theta and phi covered by module
cf0c2bc1 1202 fGeom->EmcModuleCoverage(module, tm, tM, pm, pM, AliPHOSGeometry::Degre() ) ;
6ad0bfa0 1203
1204 Double_t theta, phi ;
cf0c2bc1 1205 fGeom->EmcXtalCoverage(theta, phi, AliPHOSGeometry::Degre() ) ;
6ad0bfa0 1206
1207 Int_t tdim = (Int_t)( (tM - tm) / theta ) ;
1208 Int_t pdim = (Int_t)( (pM - pm) / phi ) ;
1209
1210 tm -= theta ;
1211 tM += theta ;
1212 pm -= phi ;
1213 pM += phi ;
1214
92862013 1215 TH2F * histoparticle = new TH2F("histoparticle", histoname,
6ad0bfa0 1216 pdim, pm, pM, tdim, tm, tM) ;
92862013 1217 histoparticle->SetStats(kFALSE) ;
6ad0bfa0 1218
1219 // Get pointers to Alice Particle TClonesArray
1220
92862013 1221 TParticle * particle;
1222 TClonesArray * particlearray = gAlice->Particles();
6ad0bfa0 1223
92862013 1224 Text_t canvasname[80];
1225 sprintf(canvasname,"Particles incident in PHOS/EMC module # %d",module) ;
1226 TCanvas * kinecanvas = new TCanvas("kinecanvas", canvasname, 650, 500) ;
6ad0bfa0 1227
1228 // get the KINE Tree
1229
92862013 1230 TTree * kine = gAlice->TreeK() ;
1231 Stat_t nParticles = kine->GetEntries() ;
1232 cout << "DisplayKineEvent > events in kine " << nParticles << endl ;
6ad0bfa0 1233
1234 // loop over particles
1235
92862013 1236 Double_t kRADDEG = 180. / TMath::Pi() ;
6ad0bfa0 1237 Int_t index1 ;
1238 Int_t nparticlein = 0 ;
92862013 1239 for (index1 = 0 ; index1 < nParticles ; index1++){
1240 Int_t nparticle = particlearray->GetEntriesFast() ;
6ad0bfa0 1241 Int_t index2 ;
1242 for( index2 = 0 ; index2 < nparticle ; index2++) {
92862013 1243 particle = (TParticle*)particlearray->UncheckedAt(index2) ;
1244 Int_t particletype = particle->GetPdgCode() ;
1245 if (testparticle == -999 || testparticle == particletype) {
1246 Double_t phi = particle->Phi() ;
1247 Double_t theta = particle->Theta() ;
6ad0bfa0 1248 Int_t mod ;
1249 Double_t x, z ;
92862013 1250 fGeom->ImpactOnEmc(theta, phi, mod, z, x) ;
1251 if ( mod == module ) {
6ad0bfa0 1252 nparticlein++ ;
b2a60966 1253 if (particle->Energy() > fClu->GetEmcClusteringThreshold() )
1254 histoparticle->Fill(phi*kRADDEG, theta*kRADDEG, particle->Energy() ) ;
6ad0bfa0 1255 }
1256 }
1257 }
1258 }
92862013 1259 kinecanvas->Draw() ;
1260 histoparticle->Draw("color") ;
1261 TPaveText * pavetext = new TPaveText(294, 100, 300, 101);
6ad0bfa0 1262 Text_t text[40] ;
1263 sprintf(text, "Particles: %d ", nparticlein) ;
92862013 1264 pavetext->AddText(text) ;
1265 pavetext->Draw() ;
1266 kinecanvas->Update();
6ad0bfa0 1267
1268}
1269//____________________________________________________________________________
1270void AliPHOSAnalyze::DisplayRecParticles()
1271{
b2a60966 1272 // Display reconstructed particles in global Alice(theta, phi) coordinates.
1273 // One PHOS module at the time.
1274 // Click on symbols indicate the reconstructed particle type.
1275
6ad0bfa0 1276 if (fEvt == -999) {
15605d3c 1277 cout << "DisplayRecParticles > Analyze an event first ... (y/n) " ;
6ad0bfa0 1278 Text_t answer[1] ;
1279 cin >> answer ; cout << answer ;
46b146ca 1280// if ( answer == "y" )
1281// AnalyzeOneEvent() ;
6ad0bfa0 1282 }
1283 if (fEvt != -999) {
1284
92862013 1285 Int_t module ;
15605d3c 1286 cout << "DisplayRecParticles > which module (1-5, -1: all) ? " ;
92862013 1287 cin >> module ; cout << module << endl ;
1288 Text_t histoname[80] ;
1289 sprintf(histoname,"Event %d: Reconstructed particles in module %d", fEvt, module) ;
1290 Double_t tm, tM, pm, pM ; // min and Max theta and phi covered by module
cf0c2bc1 1291 fGeom->EmcModuleCoverage(module, tm, tM, pm, pM, AliPHOSGeometry::Degre() ) ;
6ad0bfa0 1292 Double_t theta, phi ;
cf0c2bc1 1293 fGeom->EmcXtalCoverage(theta, phi, AliPHOSGeometry::Degre() ) ;
6ad0bfa0 1294 Int_t tdim = (Int_t)( (tM - tm) / theta ) ;
1295 Int_t pdim = (Int_t)( (pM - pm) / phi ) ;
1296 tm -= theta ;
1297 tM += theta ;
1298 pm -= phi ;
92862013 1299 TH2F * histoRparticle = new TH2F("histoRparticle", histoname,
6ad0bfa0 1300 pdim, pm, pM, tdim, tm, tM) ;
92862013 1301 histoRparticle->SetStats(kFALSE) ;
1302 Text_t canvasname[80] ;
1303 sprintf(canvasname, "Reconstructed particles in PHOSmodule # %d", module) ;
1304 TCanvas * rparticlecanvas = new TCanvas("RparticleCanvas", canvasname, 650, 500) ;
fc879520 1305 AliPHOSRecParticle::RecParticlesList * rpl = *fPHOS->RecParticles() ;
92862013 1306 Int_t nRecParticles = rpl->GetEntries() ;
1307 Int_t nRecParticlesInModule = 0 ;
6ad0bfa0 1308 TIter nextRecPart(rpl) ;
1309 AliPHOSRecParticle * rp ;
92862013 1310 cout << "DisplayRecParticles > " << nRecParticles << " reconstructed particles " << endl ;
1311 Double_t kRADDEG = 180. / TMath::Pi() ;
6ad0bfa0 1312 while ( (rp = (AliPHOSRecParticle *)nextRecPart() ) ) {
1313 AliPHOSTrackSegment * ts = rp->GetPHOSTrackSegment() ;
b2a60966 1314 if ( ts->GetPHOSMod() == module ) {
1315 Int_t numberofprimaries = 0 ;
134ce69a 1316 Int_t * listofprimaries = 0;
1317 rp->GetPrimaries(numberofprimaries) ;
b2a60966 1318 cout << "Number of primaries = " << numberofprimaries << endl ;
1319 Int_t index ;
1320 for ( index = 0 ; index < numberofprimaries ; index++)
1321 cout << " primary # " << index << " = " << listofprimaries[index] << endl ;
1322
92862013 1323 nRecParticlesInModule++ ;
1324 Double_t theta = rp->Theta() * kRADDEG ;
1325 Double_t phi = rp->Phi() * kRADDEG ;
6ad0bfa0 1326 Double_t energy = rp->Energy() ;
92862013 1327 histoRparticle->Fill(phi, theta, energy) ;
6ad0bfa0 1328 }
1329 }
92862013 1330 histoRparticle->Draw("color") ;
15605d3c 1331
1332 nextRecPart.Reset() ;
1333 while ( (rp = (AliPHOSRecParticle *)nextRecPart() ) ) {
1334 AliPHOSTrackSegment * ts = rp->GetPHOSTrackSegment() ;
1335 if ( ts->GetPHOSMod() == module )
1336 rp->Draw("P") ;
1337 }
1338
6ad0bfa0 1339 Text_t text[80] ;
92862013 1340 sprintf(text, "reconstructed particles: %d", nRecParticlesInModule) ;
1341 TPaveText * pavetext = new TPaveText(292, 100, 300, 101);
1342 pavetext->AddText(text) ;
1343 pavetext->Draw() ;
1344 rparticlecanvas->Update() ;
6ad0bfa0 1345 }
1346}
1347
1348//____________________________________________________________________________
1349void AliPHOSAnalyze::DisplayRecPoints()
1350{
b2a60966 1351 // Display reconstructed points in local PHOS-module (x, z) coordinates.
1352 // One PHOS module at the time.
1353 // Click on symbols displays the EMC cluster, or PPSD information.
1354
6ad0bfa0 1355 if (fEvt == -999) {
1356 cout << "DisplayRecPoints > Analyze an event first ... (y/n) " ;
1357 Text_t answer[1] ;
1358 cin >> answer ; cout << answer ;
46b146ca 1359// if ( answer == "y" )
1360// AnalyzeOneEvent() ;
6ad0bfa0 1361 }
1362 if (fEvt != -999) {
1363
92862013 1364 Int_t module ;
6ad0bfa0 1365 cout << "DisplayRecPoints > which module (1-5, -1: all) ? " ;
92862013 1366 cin >> module ; cout << module << endl ;
6ad0bfa0 1367
92862013 1368 Text_t canvasname[80];
1369 sprintf(canvasname,"Digits in PHOS/EMC module # %d",module) ;
1370 TCanvas * modulecanvas = new TCanvas("module", canvasname, 650, 500) ;
1371 modulecanvas->Draw() ;
6ad0bfa0 1372
92862013 1373 //=========== Creating 2d-histogram of the PHOS module
6ad0bfa0 1374 // a little bit junkie but is used to test Geom functinalities
1375
92862013 1376 Double_t tm, tM, pm, pM ; // min and Max theta and phi covered by module
6ad0bfa0 1377
92862013 1378 fGeom->EmcModuleCoverage(module, tm, tM, pm, pM);
6ad0bfa0 1379 // convert angles into coordinates local to the EMC module of interest
1380
92862013 1381 Int_t emcModuleNumber ;
1382 Double_t emcModulexm, emcModulezm ; // minimum local coordinate in a given EMCA module
1383 Double_t emcModulexM, emcModulezM ; // maximum local coordinate in a given EMCA module
1384 fGeom->ImpactOnEmc(tm, pm, emcModuleNumber, emcModulezm, emcModulexm) ;
1385 fGeom->ImpactOnEmc(tM, pM, emcModuleNumber, emcModulezM, emcModulexM) ;
1386 Int_t xdim = (Int_t)( ( emcModulexM - emcModulexm ) / fGeom->GetCrystalSize(0) ) ;
1387 Int_t zdim = (Int_t)( ( emcModulezM - emcModulezm ) / fGeom->GetCrystalSize(2) ) ;
1388 Float_t xmin = emcModulexm - fGeom->GetCrystalSize(0) ;
1389 Float_t xMax = emcModulexM + fGeom->GetCrystalSize(0) ;
1390 Float_t zmin = emcModulezm - fGeom->GetCrystalSize(2) ;
1391 Float_t zMax = emcModulezM + fGeom->GetCrystalSize(2) ;
1392 Text_t histoname[80];
1393 sprintf(histoname,"Event %d: Digits and RecPoints in module %d", fEvt, module) ;
1394 TH2F * hModule = new TH2F("HistoReconstructed", histoname,
6ad0bfa0 1395 xdim, xmin, xMax, zdim, zmin, zMax) ;
1396 hModule->SetMaximum(2.0);
1397 hModule->SetMinimum(0.0);
1398 hModule->SetStats(kFALSE);
1399
1400 TIter next(fPHOS->Digits()) ;
92862013 1401 Float_t energy, y, z;
1402 Float_t etot=0.;
1403 Int_t relid[4]; Int_t nDigits = 0 ;
6ad0bfa0 1404 AliPHOSDigit * digit ;
78ccc411 1405
1406 // Making 2D histogram of the EMC module
6ad0bfa0 1407 while((digit = (AliPHOSDigit *)next()))
1408 {
92862013 1409 fGeom->AbsToRelNumbering(digit->GetId(), relid) ;
78ccc411 1410 if (relid[0] == module && relid[1] == 0)
5959e315 1411 {
92862013 1412 energy = fClu->Calibrate(digit->GetAmp()) ;
774e78c2 1413 cout << "Energy is " << energy << " and threshold is " << fClu->GetEmcEnergyThreshold() << endl;
78ccc411 1414 if (energy > fClu->GetEmcEnergyThreshold() ){
5959e315 1415 nDigits++ ;
1416 etot += energy ;
1417 fGeom->RelPosInModule(relid,y,z) ;
1418 hModule->Fill(y, z, energy) ;
3fc07a80 1419 }
6ad0bfa0 1420 }
1421 }
92862013 1422 cout <<"DrawRecPoints > Found in module "
1423 << module << " " << nDigits << " digits with total energy " << etot << endl ;
6ad0bfa0 1424 hModule->Draw("col2") ;
1425
92862013 1426 //=========== Cluster in module
6ad0bfa0 1427
83974468 1428 // TClonesArray * emcRP = fPHOS->EmcClusters() ;
fc879520 1429 TObjArray * emcRP = *(fPHOS->EmcRecPoints()) ;
83974468 1430
92862013 1431 etot = 0.;
1432 Int_t totalnClusters = 0 ;
1433 Int_t nClusters = 0 ;
1434 TIter nextemc(emcRP) ;
6ad0bfa0 1435 AliPHOSEmcRecPoint * emc ;
1436 while((emc = (AliPHOSEmcRecPoint *)nextemc()))
1437 {
3fc07a80 1438 // Int_t numberofprimaries ;
1439 // Int_t * primariesarray = new Int_t[10] ;
1440 // emc->GetPrimaries(numberofprimaries, primariesarray) ;
92862013 1441 totalnClusters++ ;
1442 if ( emc->GetPHOSMod() == module )
6ad0bfa0 1443 {
92862013 1444 nClusters++ ;
1445 energy = emc->GetTotalEnergy() ;
1446 etot+= energy ;
26d4b141 1447 emc->Draw("M") ;
6ad0bfa0 1448 }
1449 }
92862013 1450 cout << "DrawRecPoints > Found " << totalnClusters << " EMC Clusters in PHOS" << endl ;
1451 cout << "DrawRecPoints > Found in module " << module << " " << nClusters << " EMC Clusters " << endl ;
1452 cout << "DrawRecPoints > total energy " << etot << endl ;
6ad0bfa0 1453
92862013 1454 TPaveText * pavetext = new TPaveText(22, 80, 83, 90);
6ad0bfa0 1455 Text_t text[40] ;
92862013 1456 sprintf(text, "digits: %d; clusters: %d", nDigits, nClusters) ;
1457 pavetext->AddText(text) ;
1458 pavetext->Draw() ;
1459 modulecanvas->Update();
6ad0bfa0 1460
92862013 1461 //=========== Cluster in module PPSD Down
6ad0bfa0 1462
83974468 1463 // TClonesArray * ppsdRP = fPHOS->PpsdClusters() ;
fc879520 1464 TObjArray * ppsdRP = *(fPHOS->PpsdRecPoints() );
83974468 1465
92862013 1466 etot = 0.;
1467 TIter nextPpsd(ppsdRP) ;
1468 AliPHOSPpsdRecPoint * ppsd ;
1469 while((ppsd = (AliPHOSPpsdRecPoint *)nextPpsd()))
6ad0bfa0 1470 {
92862013 1471 totalnClusters++ ;
1472 if ( ppsd->GetPHOSMod() == module )
6ad0bfa0 1473 {
92862013 1474 nClusters++ ;
1475 energy = ppsd->GetEnergy() ;
1476 etot+=energy ;
1477 if (!ppsd->GetUp()) ppsd->Draw("P") ;
6ad0bfa0 1478 }
1479 }
92862013 1480 cout << "DrawRecPoints > Found " << totalnClusters << " Ppsd Down Clusters in PHOS" << endl ;
1481 cout << "DrawRecPoints > Found in module " << module << " " << nClusters << " Ppsd Down Clusters " << endl ;
1482 cout << "DrawRecPoints > total energy " << etot << endl ;
6ad0bfa0 1483
92862013 1484 //=========== Cluster in module PPSD Up
6ad0bfa0 1485
fc879520 1486 ppsdRP = *(fPHOS->PpsdRecPoints()) ;
83974468 1487
92862013 1488 etot = 0.;
1489 TIter nextPpsdUp(ppsdRP) ;
1490 while((ppsd = (AliPHOSPpsdRecPoint *)nextPpsdUp()))
6ad0bfa0 1491 {
92862013 1492 totalnClusters++ ;
1493 if ( ppsd->GetPHOSMod() == module )
6ad0bfa0 1494 {
92862013 1495 nClusters++ ;
1496 energy = ppsd->GetEnergy() ;
1497 etot+=energy ;
1498 if (ppsd->GetUp()) ppsd->Draw("P") ;
6ad0bfa0 1499 }
1500 }
92862013 1501 cout << "DrawRecPoints > Found " << totalnClusters << " Ppsd Up Clusters in PHOS" << endl ;
1502 cout << "DrawRecPoints > Found in module " << module << " " << nClusters << " Ppsd Up Clusters " << endl ;
1503 cout << "DrawRecPoints > total energy " << etot << endl ;
6ad0bfa0 1504
1505 } // if !-999
1506}
1507
1508//____________________________________________________________________________
1509void AliPHOSAnalyze::DisplayTrackSegments()
1510{
b2a60966 1511 // Display track segments in local PHOS-module (x, z) coordinates.
1512 // One PHOS module at the time.
1513 // One symbol per PHOS subsystem: EMC, upper PPSD, lower PPSD.
1514
6ad0bfa0 1515 if (fEvt == -999) {
1516 cout << "DisplayTrackSegments > Analyze an event first ... (y/n) " ;
1517 Text_t answer[1] ;
1518 cin >> answer ; cout << answer ;
46b146ca 1519// if ( answer == "y" )
1520// AnalyzeOneEvent() ;
6ad0bfa0 1521 }
1522 if (fEvt != -999) {
1523
92862013 1524 Int_t module ;
6ad0bfa0 1525 cout << "DisplayTrackSegments > which module (1-5, -1: all) ? " ;
92862013 1526 cin >> module ; cout << module << endl ;
1527 //=========== Creating 2d-histogram of the PHOS module
6ad0bfa0 1528 // a little bit junkie but is used to test Geom functinalities
1529
92862013 1530 Double_t tm, tM, pm, pM ; // min and Max theta and phi covered by module
6ad0bfa0 1531
92862013 1532 fGeom->EmcModuleCoverage(module, tm, tM, pm, pM);
6ad0bfa0 1533 // convert angles into coordinates local to the EMC module of interest
1534
92862013 1535 Int_t emcModuleNumber ;
1536 Double_t emcModulexm, emcModulezm ; // minimum local coordinate in a given EMCA module
1537 Double_t emcModulexM, emcModulezM ; // maximum local coordinate in a given EMCA module
1538 fGeom->ImpactOnEmc(tm, pm, emcModuleNumber, emcModulezm, emcModulexm) ;
1539 fGeom->ImpactOnEmc(tM, pM, emcModuleNumber, emcModulezM, emcModulexM) ;
1540 Int_t xdim = (Int_t)( ( emcModulexM - emcModulexm ) / fGeom->GetCrystalSize(0) ) ;
1541 Int_t zdim = (Int_t)( ( emcModulezM - emcModulezm ) / fGeom->GetCrystalSize(2) ) ;
1542 Float_t xmin = emcModulexm - fGeom->GetCrystalSize(0) ;
1543 Float_t xMax = emcModulexM + fGeom->GetCrystalSize(0) ;
1544 Float_t zmin = emcModulezm - fGeom->GetCrystalSize(2) ;
1545 Float_t zMax = emcModulezM + fGeom->GetCrystalSize(2) ;
1546 Text_t histoname[80];
1547 sprintf(histoname,"Event %d: Track Segments in module %d", fEvt, module) ;
1548 TH2F * histotrack = new TH2F("histotrack", histoname,
6ad0bfa0 1549 xdim, xmin, xMax, zdim, zmin, zMax) ;
92862013 1550 histotrack->SetStats(kFALSE);
1551 Text_t canvasname[80];
1552 sprintf(canvasname,"Track segments in PHOS/EMC-PPSD module # %d", module) ;
1553 TCanvas * trackcanvas = new TCanvas("TrackSegmentCanvas", canvasname, 650, 500) ;
1554 histotrack->Draw() ;
6ad0bfa0 1555
fc879520 1556 AliPHOSTrackSegment::TrackSegmentsList * trsegl = *(fPHOS->TrackSegments()) ;
6ad0bfa0 1557 AliPHOSTrackSegment * trseg ;
1558
92862013 1559 Int_t nTrackSegments = trsegl->GetEntries() ;
6ad0bfa0 1560 Int_t index ;
92862013 1561 Float_t etot = 0 ;
1562 Int_t nTrackSegmentsInModule = 0 ;
1563 for(index = 0; index < nTrackSegments ; index++){
6ad0bfa0 1564 trseg = (AliPHOSTrackSegment * )trsegl->At(index) ;
92862013 1565 etot+= trseg->GetEnergy() ;
1566 if ( trseg->GetPHOSMod() == module ) {
1567 nTrackSegmentsInModule++ ;
6ad0bfa0 1568 trseg->Draw("P");
1569 }
1570 }
1571 Text_t text[80] ;
92862013 1572 sprintf(text, "track segments: %d", nTrackSegmentsInModule) ;
1573 TPaveText * pavetext = new TPaveText(22, 80, 83, 90);
1574 pavetext->AddText(text) ;
1575 pavetext->Draw() ;
1576 trackcanvas->Update() ;
1577 cout << "DisplayTrackSegments > Found " << trsegl->GetEntries() << " Track segments with total energy "<< etot << endl ;
6ad0bfa0 1578
1579 }
1580}
1581//____________________________________________________________________________
1582Bool_t AliPHOSAnalyze::OpenRootFile(Text_t * name)
1583{
b2a60966 1584 // Open the root file named "name"
1585
1586 fRootFile = new TFile(name, "update") ;
6ad0bfa0 1587 return fRootFile->IsOpen() ;
1588}
92862013 1589//____________________________________________________________________________
46b146ca 1590// void AliPHOSAnalyze::SavingHistograms()
1591// {
1592// // Saves the histograms in a root file named "name.analyzed"
1593
1594// Text_t outputname[80] ;
1595// sprintf(outputname,"%s.analyzed",fRootFile->GetName());
1596// TFile output(outputname,"RECREATE");
1597// output.cd();
134ce69a 1598// if (fhEmcDigit )
1599// fhEmcDigit->Write() ;
1600// if (fhVetoDigit )
1601// fhVetoDigit->Write() ;
1602// if (fhConvertorDigit )
1603// fhConvertorDigit->Write() ;
1604// if (fhEmcCluster )
1605// fhEmcCluster->Write() ;
1606// if (fhVetoCluster )
1607// fhVetoCluster->Write() ;
1608// if (fhConvertorCluster )
1609// fhConvertorCluster->Write() ;
1610// if (fhConvertorEmc )
1611// fhConvertorEmc->Write() ;
1612// if (fhPhotonEnergy)
1613// fhPhotonEnergy->Write() ;
1614// if (fhPhotonPositionX)
1615// fhPhotonPositionX->Write() ;
1616// if (fhPhotonPositionY)
1617// fhPhotonPositionX->Write() ;
1618// if (fhElectronEnergy)
1619// fhElectronEnergy->Write() ;
1620// if (fhElectronPositionX)
1621// fhElectronPositionX->Write() ;
1622// if (fhElectronPositionY)
1623// fhElectronPositionX->Write() ;
1624// if (fhNeutralHadronEnergy)
1625// fhNeutralHadronEnergy->Write() ;
1626// if (fhNeutralHadronPositionX)
1627// fhNeutralHadronPositionX->Write() ;
1628// if (fhNeutralHadronPositionY)
1629// fhNeutralHadronPositionX->Write() ;
1630// if (fhNeutralEMEnergy)
1631// fhNeutralEMEnergy->Write() ;
1632// if (fhNeutralEMPositionX)
1633// fhNeutralEMPositionX->Write() ;
1634// if (fhNeutralEMPositionY)
1635// fhNeutralEMPositionX->Write() ;
1636// if (fhChargedHadronEnergy)
1637// fhChargedHadronEnergy->Write() ;
1638// if (fhChargedHadronPositionX)
1639// fhChargedHadronPositionX->Write() ;
1640// if (fhChargedHadronPositionY)
1641// fhChargedHadronPositionX->Write() ;
1642// if (fhPhotonHadronEnergy)
1643// fhPhotonHadronEnergy->Write() ;
1644// if (fhPhotonHadronPositionX)
1645// fhPhotonHadronPositionX->Write() ;
1646// if (fhPhotonHadronPositionY)
1647// fhPhotonHadronPositionX->Write() ;
1648
46b146ca 1649// output.Write();
1650// output.Close();
1651// }
134ce69a 1652//____________________________________________________________________________
46b146ca 1653void AliPHOSAnalyze::SaveHistograms()
134ce69a 1654{
1655 // Saves the histograms in a root file named "name.analyzed"
1656
1657 Text_t outputname[80] ;
1658 sprintf(outputname,"%s.analyzed",fRootFile->GetName());
1659 TFile output(outputname,"RECREATE");
1660 output.cd();
46b146ca 1661
b9bbdad1 1662 if (fhPhotonEnergy)
1663 fhPhotonEnergy->Write() ;
fc879520 1664 if (fhPhotonAllEnergy)
1665 fhPhotonAllEnergy->Write() ;
134ce69a 1666 if (fhPhotonPosition)
1667 fhPhotonPosition->Write() ;
fc879520 1668 if (fhPhotonAllPosition)
1669 fhPhotonAllPosition->Write() ;
b9bbdad1 1670 if (fhElectronEnergy)
1671 fhElectronEnergy->Write() ;
fc879520 1672 if (fhElectronAllEnergy)
1673 fhElectronAllEnergy->Write() ;
134ce69a 1674 if (fhElectronPosition)
1675 fhElectronPosition->Write() ;
fc879520 1676 if (fhElectronAllPosition)
1677 fhElectronAllPosition->Write() ;
b9bbdad1 1678 if (fhNeutralHadronEnergy)
1679 fhNeutralHadronEnergy->Write() ;
134ce69a 1680 if (fhNeutralHadronPosition)
1681 fhNeutralHadronPosition->Write() ;
b9bbdad1 1682 if (fhNeutralEMEnergy)
1683 fhNeutralEMEnergy->Write() ;
fc879520 1684 if (fhNeutralEMAllEnergy)
1685 fhNeutralEMAllEnergy->Write() ;
134ce69a 1686 if (fhNeutralEMPosition)
1687 fhNeutralEMPosition->Write() ;
fc879520 1688 if (fhNeutralEMAllPosition)
1689 fhNeutralEMAllPosition->Write() ;
b9bbdad1 1690 if (fhChargedHadronEnergy)
1691 fhChargedHadronEnergy->Write() ;
134ce69a 1692 if (fhChargedHadronPosition)
1693 fhChargedHadronPosition->Write() ;
c1d256cb 1694 if (fhPhotonHadronEnergy)
1695 fhPhotonHadronEnergy->Write() ;
134ce69a 1696 if (fhPhotonHadronPosition)
1697 fhPhotonHadronPosition->Write() ;
fc879520 1698 if (fhPhotonReg)
1699 fhPhotonReg->Write() ;
1700 if (fhAllReg)
1701 fhAllReg->Write() ;
1702 if(fhNReg)
1703 fhNReg->Write() ;
1704 if(fhNBarReg)
1705 fhNBarReg->Write() ;
1706 if(fhChargedReg)
1707 fhChargedReg->Write() ;
1708 if (fhPhotonEM)
1709 fhPhotonEM->Write() ;
1710 if (fhAllEM)
1711 fhAllEM->Write() ;
1712 if(fhNEM)
1713 fhNEM->Write() ;
1714 if(fhNBarEM)
1715 fhNBarEM->Write() ;
1716 if(fhChargedEM)
1717 fhChargedEM->Write() ;
1718 if(fhPrimary)
1719 fhPrimary->Write() ;
1720 if(fhPhotPhot)
1721 fhPhotPhot->Write() ;
1722 if(fhPhotElec)
1723 fhPhotElec->Write() ;
1724 if(fhPhotNeuH)
1725 fhPhotNeuH->Write() ;
1726 if(fhPhotNuEM)
1727 fhPhotNuEM->Write() ;
1728 if(fhPhotNuEM)
1729 fhPhotNuEM->Write() ;
1730 if(fhPhotChHa)
1731 fhPhotChHa->Write() ;
1732 if(fhPhotGaHa)
1733 fhPhotGaHa->Write() ;
46b146ca 1734 if(fhEnergyCorrelations)
1735 fhEnergyCorrelations->Write() ;
1736
92862013 1737 output.Write();
1738 output.Close();
1739}
46b146ca 1740//____________________________________________________________________________
1741void AliPHOSAnalyze::ResetHistograms()
1742{
1743 fhEnergyCorrelations = 0 ; //Energy correlations between Eloss in Convertor and PPSD(2)
1744
1745 fhEmcDigit = 0 ; // Histo of digit energies in the Emc
1746 fhVetoDigit = 0 ; // Histo of digit energies in the Veto
1747 fhConvertorDigit = 0 ; // Histo of digit energies in the Convertor
1748 fhEmcCluster = 0 ; // Histo of Cluster energies in Emc
1749 fhVetoCluster = 0 ; // Histo of Cluster energies in Veto
1750 fhConvertorCluster = 0 ; // Histo of Cluster energies in Convertor
1751 fhConvertorEmc = 0 ; // 2d Convertor versus Emc energies
1752
1753 fhPhotonEnergy = 0 ; // Spectrum of detected photons with photon primary
1754 fhPhotonAllEnergy = 0 ; // Total spectrum of detected photons
1755 fhElectronEnergy = 0 ; // Spectrum of detected electrons with electron primary
1756 fhElectronAllEnergy = 0 ; // Total spectrum of detected electrons
1757 fhNeutralHadronEnergy = 0 ; // Spectrum of detected neutral hadron
1758 fhNeutralEMEnergy = 0 ; // Spectrum of detected neutral EM with EM primary
1759 fhNeutralEMAllEnergy = 0 ; // Spectrum of detected neutral EM
1760 fhChargedHadronEnergy = 0 ; // Spectrum of detected charged
1761 fhPhotonHadronEnergy = 0 ; // Spectrum of detected Photon-Hadron
1762 fhPhotonPosition = 0 ; // Position Resolution of photons with photon primary
1763 fhPhotonAllPosition = 0 ; // Position Resolution of photons
1764 fhElectronPosition = 0 ; // Position Resolution of electrons with electron primary
1765 fhElectronAllPosition = 0 ; // Position Resolution of electrons
1766 fhNeutralHadronPosition = 0 ; // Position Resolution of neutral hadron
1767 fhNeutralEMPosition = 0 ; // Position Resolution of neutral EM with EM primary
1768 fhNeutralEMAllPosition = 0 ; // Position Resolution of neutral EM
1769 fhChargedHadronPosition = 0 ; // Position Resolution of charged
1770 fhPhotonHadronPosition = 0 ; // Position Resolution of Photon-Hadron
1771 fhPhotonPositionY = 0 ; // Y distribution of detected photons
1772 fhElectronPositionY = 0 ; // Y distribution of detected electrons
1773 fhNeutralHadronPositionY = 0 ; // Y distribution of detected neutral hadron
1774 fhNeutralEMPositionY = 0 ; // Y distribution of detected neutral EM
1775 fhChargedHadronPositionY = 0 ; // Y distribution of detected charged
1776 fhPhotonHadronPositionY = 0 ; // Y distribution of detected Photon-Hadron
1777 fhPhotonReg = 0 ;
1778 fhAllReg = 0 ;
1779 fhNReg = 0 ;
1780 fhNBarReg = 0 ;
1781 fhChargedReg = 0 ;
1782 fhPhotonEM = 0 ;
1783 fhAllEM = 0 ;
1784 fhNEM = 0 ;
1785 fhNBarEM = 0 ;
1786 fhChargedEM = 0 ;
1787 fhPrimary = 0 ;
1788
1789 fhPhotPhot = 0 ;
1790 fhPhotElec = 0 ;
1791 fhPhotNeuH = 0 ;
1792 fhPhotNuEM = 0 ;
1793 fhPhotChHa = 0 ;
1794 fhPhotGaHa = 0 ;
1795
134ce69a 1796
46b146ca 1797}
134ce69a 1798
1799