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