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