]>
Commit | Line | Data |
---|---|---|
6ad0bfa0 | 1 | /************************************************************************** |
2 | * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. * | |
3 | * * | |
4 | * Author: The ALICE Off-line Project. * | |
5 | * Contributors are mentioned in the code where appropriate. * | |
6 | * * | |
7 | * Permission to use, copy, modify and distribute this software and its * | |
8 | * documentation strictly for non-commercial purposes is hereby granted * | |
9 | * without fee, provided that the above copyright notice appears in all * | |
10 | * copies and that both the copyright notice and this permission notice * | |
11 | * appear in the supporting documentation. The authors make no claims * | |
12 | * about the suitability of this software for any purpose. It is * | |
13 | * provided "as is" without express or implied warranty. * | |
14 | **************************************************************************/ | |
15 | ||
b2a60966 | 16 | /* $Id$ */ |
17 | ||
6ad0bfa0 | 18 | //_________________________________________________________________________ |
5f20d3fb | 19 | // Algorythm class to analyze PHOSv1 events: |
b2a60966 | 20 | // Construct histograms and displays them. |
21 | // Use the macro EditorBar.C for best access to the functionnalities | |
22 | // | |
2aad621e | 23 | //*-- Author: Y. Schutz (SUBATECH) & Gines Martinez (SUBATECH) |
6ad0bfa0 | 24 | ////////////////////////////////////////////////////////////////////////////// |
25 | ||
26 | // --- ROOT system --- | |
27 | ||
92862013 | 28 | #include "TFile.h" |
29 | #include "TH1.h" | |
30 | #include "TPad.h" | |
6ad0bfa0 | 31 | #include "TH2.h" |
32 | #include "TParticle.h" | |
33 | #include "TClonesArray.h" | |
34 | #include "TTree.h" | |
35 | #include "TMath.h" | |
36 | #include "TCanvas.h" | |
37 | ||
38 | // --- Standard library --- | |
39 | ||
de9ec31b | 40 | #include <iostream.h> |
41 | #include <stdio.h> | |
92862013 | 42 | |
6ad0bfa0 | 43 | // --- AliRoot header files --- |
44 | ||
45 | #include "AliRun.h" | |
46 | #include "AliPHOSAnalyze.h" | |
47 | #include "AliPHOSClusterizerv1.h" | |
48 | #include "AliPHOSTrackSegmentMakerv1.h" | |
26d4b141 | 49 | #include "AliPHOSPIDv1.h" |
6ad0bfa0 | 50 | #include "AliPHOSReconstructioner.h" |
51 | #include "AliPHOSDigit.h" | |
52 | #include "AliPHOSTrackSegment.h" | |
53 | #include "AliPHOSRecParticle.h" | |
83974468 | 54 | #include "AliPHOSIndexToObject.h" |
6ad0bfa0 | 55 | |
56 | ClassImp(AliPHOSAnalyze) | |
57 | ||
58 | ||
59 | //____________________________________________________________________________ | |
60 | AliPHOSAnalyze::AliPHOSAnalyze() | |
61 | { | |
b2a60966 | 62 | // default ctor (useless) |
6ad0bfa0 | 63 | |
64 | fRootFile = 0 ; | |
65 | } | |
66 | ||
67 | //____________________________________________________________________________ | |
68 | AliPHOSAnalyze::AliPHOSAnalyze(Text_t * name) | |
69 | { | |
83974468 | 70 | // ctor: analyze events from root file "name" |
6ad0bfa0 | 71 | |
92862013 | 72 | Bool_t ok = OpenRootFile(name) ; |
73 | if ( !ok ) { | |
6ad0bfa0 | 74 | cout << " AliPHOSAnalyze > Error opening " << name << endl ; |
75 | } | |
76 | else { | |
77 | gAlice = (AliRun*) fRootFile->Get("gAlice"); | |
5f20d3fb | 78 | fPHOS = (AliPHOSv1 *)gAlice->GetDetector("PHOS") ; |
134ce69a | 79 | fGeom = AliPHOSGeometry::GetInstance( fPHOS->GetGeometry()->GetName(), fPHOS->GetGeometry()->GetTitle() ) ; |
aafe457d | 80 | |
6ad0bfa0 | 81 | fEvt = -999 ; |
82 | } | |
83 | } | |
84 | ||
88714635 | 85 | //____________________________________________________________________________ |
86 | AliPHOSAnalyze::AliPHOSAnalyze(const AliPHOSAnalyze & ana) | |
87 | { | |
88 | // copy ctor | |
89 | ( (AliPHOSAnalyze &)ana ).Copy(*this) ; | |
90 | } | |
91 | ||
92 | //____________________________________________________________________________ | |
191c1c41 | 93 | void AliPHOSAnalyze::Copy(TObject & obj) |
88714635 | 94 | { |
95 | // copy an analysis into an other one | |
96 | TObject::Copy(obj) ; | |
97 | // I do nothing more because the copy is silly but the Code checkers requires one | |
98 | } | |
99 | ||
6ad0bfa0 | 100 | //____________________________________________________________________________ |
101 | AliPHOSAnalyze::~AliPHOSAnalyze() | |
102 | { | |
103 | // dtor | |
104 | ||
83974468 | 105 | if (fRootFile->IsOpen() ) |
106 | fRootFile->Close() ; | |
6ad0bfa0 | 107 | delete fRootFile ; |
108 | fRootFile = 0 ; | |
109 | ||
110 | delete fPHOS ; | |
111 | fPHOS = 0 ; | |
112 | ||
113 | delete fClu ; | |
114 | fClu = 0 ; | |
115 | ||
26d4b141 | 116 | delete fPID ; |
117 | fPID = 0 ; | |
6ad0bfa0 | 118 | |
119 | delete fRec ; | |
120 | fRec = 0 ; | |
121 | ||
122 | delete fTrs ; | |
123 | fTrs = 0 ; | |
124 | ||
125 | } | |
126 | ||
127 | //____________________________________________________________________________ | |
128 | void AliPHOSAnalyze::AnalyzeOneEvent(Int_t evt) | |
129 | { | |
b2a60966 | 130 | // analyze one single event with id=evt |
131 | ||
132 | TStopwatch ts ; | |
133 | ts.Start() ; | |
92862013 | 134 | Bool_t ok = Init(evt) ; |
6ad0bfa0 | 135 | |
92862013 | 136 | if ( ok ) { |
6ad0bfa0 | 137 | //=========== Get the number of entries in the Digits array |
138 | ||
139 | Int_t nId = fPHOS->Digits()->GetEntries(); | |
140 | printf("AnalyzeOneEvent > Number of entries in the Digit array is %d \n",nId); | |
141 | ||
142 | //=========== Do the reconstruction | |
143 | ||
144 | cout << "AnalyzeOneEvent > Found " << nId << " digits in PHOS" << endl ; | |
145 | ||
146 | fPHOS->Reconstruction(fRec); | |
134ce69a | 147 | |
6ad0bfa0 | 148 | // =========== End of reconstruction |
b2a60966 | 149 | |
cfaa8d98 | 150 | // Deleting fClu, fTrs, fPID et fRec |
151 | fClu->Delete(); | |
152 | fTrs->Delete(); | |
153 | fPID->Delete(); | |
154 | fRec->Delete(); | |
155 | ||
b2a60966 | 156 | // =========== Write the root file |
157 | ||
83974468 | 158 | |
b2a60966 | 159 | // =========== Finish |
160 | ||
6ad0bfa0 | 161 | cout << "AnalyzeOneEvent > event # " << fEvt << " processed" << endl ; |
92862013 | 162 | } // ok |
6ad0bfa0 | 163 | else |
164 | cout << "AnalyzeOneEvent > filed to process event # " << evt << endl ; | |
b2a60966 | 165 | |
166 | ts.Stop() ; cout << "CPU time = " << ts.CpuTime() << endl ; | |
167 | cout << "Real time = " << ts.RealTime() << endl ; | |
6ad0bfa0 | 168 | } |
169 | ||
92862013 | 170 | //____________________________________________________________________________ |
b2a60966 | 171 | void AliPHOSAnalyze::AnalyzeManyEvents(Int_t Nevents, Int_t module) |
92862013 | 172 | { |
b2a60966 | 173 | // analyzes Nevents events in a single PHOS module |
92862013 | 174 | |
175 | if ( fRootFile == 0 ) | |
176 | cout << "AnalyzeManyEvents > " << "Root File not openned" << endl ; | |
177 | else | |
178 | { | |
179 | //========== Get AliRun object from file | |
180 | gAlice = (AliRun*) fRootFile->Get("gAlice") ; | |
181 | //=========== Get the PHOS object and associated geometry from the file | |
5f20d3fb | 182 | fPHOS = (AliPHOSv1 *)gAlice->GetDetector("PHOS") ; |
92862013 | 183 | fGeom = AliPHOSGeometry::GetInstance( fPHOS->GetGeometry()->GetName(), fPHOS->GetGeometry()->GetTitle() ); |
184 | //========== Booking Histograms | |
185 | cout << "AnalyzeManyEvents > " << "Booking Histograms" << endl ; | |
186 | BookingHistograms(); | |
187 | Int_t ievent; | |
6a3f1304 | 188 | Int_t relid[4] ; |
92862013 | 189 | AliPHOSDigit * digit ; |
134ce69a | 190 | // AliPHOSEmcRecPoint * emc ; |
191 | // AliPHOSPpsdRecPoint * ppsd ; | |
38f8ae6d | 192 | // AliPHOSTrackSegment * tracksegment ; |
193 | AliPHOSRecParticle * recparticle; | |
92862013 | 194 | for ( ievent=0; ievent<Nevents; ievent++) |
195 | { | |
196 | if (ievent==0) cout << "AnalyzeManyEvents > " << "Starting Analyzing " << endl ; | |
197 | //========== Create the Clusterizer | |
198 | fClu = new AliPHOSClusterizerv1() ; | |
2aad621e | 199 | fClu->SetEmcEnergyThreshold(0.05) ; |
3fc07a80 | 200 | fClu->SetEmcClusteringThreshold(0.50) ; |
92862013 | 201 | fClu->SetPpsdEnergyThreshold (0.0000002) ; |
202 | fClu->SetPpsdClusteringThreshold(0.0000001) ; | |
203 | fClu->SetLocalMaxCut(0.03) ; | |
204 | fClu->SetCalibrationParameters(0., 0.00000001) ; | |
205 | //========== Creates the track segment maker | |
6727be7e | 206 | fTrs = new AliPHOSTrackSegmentMakerv1() ; |
134ce69a | 207 | // fTrs->UnsetUnfoldFlag() ; |
26d4b141 | 208 | //========== Creates the particle identifier |
209 | fPID = new AliPHOSPIDv1() ; | |
6a3f1304 | 210 | fPID->SetShowerProfileCuts(0.3, 1.8, 0.3, 1.8 ) ; |
09fc14a0 | 211 | fPID->Print() ; |
92862013 | 212 | //========== Creates the Reconstructioner |
26d4b141 | 213 | fRec = new AliPHOSReconstructioner(fClu, fTrs, fPID) ; |
de9ec31b | 214 | //========== Event Number> |
215 | if ( ( log10((Float_t)(ievent+1)) - (Int_t)(log10((Float_t)(ievent+1))) ) == 0. ) | |
216 | cout << "AnalyzeManyEvents > " << "Event is " << ievent << endl ; | |
92862013 | 217 | //=========== Connects the various Tree's for evt |
218 | gAlice->GetEvent(ievent); | |
219 | //=========== Gets the Digit TTree | |
220 | gAlice->TreeD()->GetEvent(0) ; | |
221 | //=========== Gets the number of entries in the Digits array | |
222 | TIter nextdigit(fPHOS->Digits()) ; | |
223 | while( ( digit = (AliPHOSDigit *)nextdigit() ) ) | |
224 | { | |
225 | fGeom->AbsToRelNumbering(digit->GetId(), relid) ; | |
226 | if (fClu->IsInEmc(digit)) fhEmcDigit->Fill(fClu->Calibrate(digit->GetAmp())) ; | |
227 | else | |
228 | { | |
229 | if (relid[1]<17) fhVetoDigit->Fill(fClu->Calibrate(digit->GetAmp())); | |
230 | if (relid[1]>16) fhConvertorDigit->Fill(fClu->Calibrate(digit->GetAmp())); | |
231 | } | |
232 | } | |
233 | //=========== Do the reconstruction | |
234 | fPHOS->Reconstruction(fRec); | |
134ce69a | 235 | |
236 | // //=========== Cluster in module | |
237 | // TIter nextEmc(fPHOS->EmcRecPoints() ) ; | |
238 | // while((emc = (AliPHOSEmcRecPoint *)nextEmc())) | |
239 | // { | |
240 | // if ( emc->GetPHOSMod() == module ) | |
241 | // { | |
242 | // fhEmcCluster->Fill( emc->GetTotalEnergy() ); | |
243 | // TIter nextPpsd( fPHOS->PpsdRecPoints()) ; | |
244 | // while((ppsd = (AliPHOSPpsdRecPoint *)nextPpsd())) | |
245 | // { | |
246 | // if ( ppsd->GetPHOSMod() == module ) | |
247 | // { | |
248 | // if (!ppsd->GetUp()) fhConvertorEmc->Fill(emc->GetTotalEnergy(),ppsd->GetTotalEnergy()) ; | |
249 | // } | |
250 | // } | |
251 | // } | |
252 | // } | |
253 | // //=========== Cluster in module PPSD Down | |
254 | // TIter nextPpsd(fPHOS->PpsdRecPoints() ) ; | |
255 | // while((ppsd = (AliPHOSPpsdRecPoint *)nextPpsd())) | |
256 | // { | |
257 | // if ( ppsd->GetPHOSMod() == module ) | |
258 | // { | |
259 | // if (!ppsd->GetUp()) fhConvertorCluster->Fill(ppsd->GetTotalEnergy()) ; | |
260 | // if (ppsd->GetUp()) fhVetoCluster ->Fill(ppsd->GetTotalEnergy()) ; | |
261 | // } | |
262 | // } | |
92862013 | 263 | //========== TRackSegments in the event |
38f8ae6d | 264 | TIter nextRecParticle(fPHOS->RecParticles() ) ; |
265 | while((recparticle = (AliPHOSRecParticle *)nextRecParticle())) | |
92862013 | 266 | { |
38f8ae6d | 267 | if ( recparticle->GetPHOSTrackSegment()->GetPHOSMod() == module ) |
92862013 | 268 | { |
b2a60966 | 269 | cout << "Particle type is " << recparticle->GetType() << endl ; |
270 | Int_t numberofprimaries = 0 ; | |
271 | Int_t * listofprimaries = recparticle->GetPrimaries(numberofprimaries) ; | |
272 | cout << "Number of primaries = " << numberofprimaries << endl ; | |
273 | Int_t index ; | |
274 | for ( index = 0 ; index < numberofprimaries ; index++) | |
275 | cout << " primary # " << index << " = " << listofprimaries[index] << endl ; | |
134ce69a | 276 | // switch(recparticle->GetType()) |
277 | // { | |
278 | // case AliPHOSFastRecParticle::kGAMMA: | |
279 | // fhPhotonEnergy->Fill(recparticle->Energy() ) ; | |
280 | // //fhPhotonPositionX->Fill(recpart. ) ; | |
281 | // //fhPhotonPositionY->Fill(recpart. ) ; | |
282 | // cout << "PHOTON" << endl; | |
283 | // break; | |
284 | // case AliPHOSFastRecParticle::kELECTRON: | |
285 | // fhElectronEnergy->Fill(recparticle->Energy() ) ; | |
286 | // //fhElectronPositionX->Fill(recpart. ) ; | |
287 | // //fhElectronPositionY->Fill(recpart. ) ; | |
288 | // cout << "ELECTRON" << endl; | |
289 | // break; | |
290 | // case AliPHOSFastRecParticle::kNEUTRALHA: | |
291 | // fhNeutralHadronEnergy->Fill(recparticle->Energy() ) ; | |
292 | // //fhNeutralHadronPositionX->Fill(recpart. ) ; | |
293 | // //fhNeutralHadronPositionY->Fill(recpart. ) ; | |
294 | // cout << "NEUTRAl HADRON" << endl; | |
295 | // break ; | |
296 | // case AliPHOSFastRecParticle::kNEUTRALEM: | |
297 | // fhNeutralEMEnergy->Fill(recparticle->Energy() ) ; | |
298 | // //fhNeutralEMPositionX->Fill(recpart. ) ; | |
299 | // //fhNeutralEMPositionY->Fill(recpart. ) ; | |
300 | // //cout << "NEUTRAL EM" << endl; | |
301 | // break ; | |
302 | // case AliPHOSFastRecParticle::kCHARGEDHA: | |
303 | // fhChargedHadronEnergy->Fill(recparticle->Energy() ) ; | |
304 | // //fhChargedHadronPositionX->Fill(recpart. ) ; | |
305 | // //fhChargedHadronPositionY->Fill(recpart. ) ; | |
306 | // cout << "CHARGED HADRON" << endl; | |
307 | // break ; | |
308 | // case AliPHOSFastRecParticle::kGAMMAHA: | |
309 | // fhPhotonHadronEnergy->Fill(recparticle->Energy() ) ; | |
310 | // //fhPhotonHadronPositionX->Fill(recpart. ) ; | |
311 | // //fhPhotonHadronPositionY->Fill(recpart. ) ; | |
312 | // cout << "PHOTON HADRON" << endl; | |
313 | // break ; | |
314 | // } | |
315 | } | |
316 | } | |
317 | // Deleting fClu, fTrs, fPID et fRec | |
318 | fClu->Delete(); | |
319 | fTrs->Delete(); | |
320 | fPID->Delete(); | |
321 | fRec->Delete(); | |
322 | ||
323 | } // endfor | |
324 | SavingHistograms(); | |
325 | } // endif | |
326 | } // endfunction | |
327 | ||
328 | //____________________________________________________________________________ | |
329 | void AliPHOSAnalyze::AnalyzeResolutions(Int_t Nevents ) | |
330 | { | |
331 | // analyzes Nevents events and calculate Energy and Position resolution as well as | |
332 | // probaility of correct indentifiing of the incident particle | |
333 | ||
334 | if ( fRootFile == 0 ) | |
335 | cout << "AnalyzeResolutions > " << "Root File not openned" << endl ; | |
336 | else | |
337 | { | |
338 | //========== Get AliRun object from file | |
339 | gAlice = (AliRun*) fRootFile->Get("gAlice") ; | |
340 | //=========== Get the PHOS object and associated geometry from the file | |
341 | fPHOS = (AliPHOSv1 *)gAlice->GetDetector("PHOS") ; | |
342 | fGeom = AliPHOSGeometry::GetInstance( fPHOS->GetGeometry()->GetName(), fPHOS->GetGeometry()->GetTitle() ); | |
343 | ||
344 | //========== Initializes the Index to Object converter | |
345 | fObjGetter = AliPHOSIndexToObject::GetInstance(fPHOS) ; | |
346 | ||
347 | //========== Booking Histograms | |
348 | cout << "AnalyzeResolutions > " << "Booking Histograms" << endl ; | |
349 | BookResolutionHistograms(); | |
350 | ||
351 | Int_t ievent; | |
352 | for ( ievent=0; ievent<Nevents; ievent++) | |
353 | { | |
354 | if (ievent==0) cout << "AnalyzeResolutions > " << "Starting Analyzing " << endl ; | |
355 | ||
356 | //========== Create the Clusterizer | |
357 | fClu = new AliPHOSClusterizerv1() ; | |
358 | fClu->SetEmcEnergyThreshold(0.05) ; | |
359 | fClu->SetEmcClusteringThreshold(0.50) ; | |
360 | fClu->SetPpsdEnergyThreshold (0.0000002) ; | |
361 | fClu->SetPpsdClusteringThreshold(0.0000001) ; | |
362 | fClu->SetLocalMaxCut(0.03) ; | |
363 | fClu->SetCalibrationParameters(0., 0.00000001) ; | |
364 | ||
365 | //========== Creates the track segment maker | |
366 | fTrs = new AliPHOSTrackSegmentMakerv1() ; | |
367 | // fTrs->UnsetUnfoldFlag() ; | |
368 | ||
369 | //========== Creates the particle identifier | |
370 | fPID = new AliPHOSPIDv1() ; | |
371 | fPID->SetShowerProfileCuts(0.3, 1.8, 0.3, 1.8 ) ; | |
372 | ||
373 | //========== Creates the Reconstructioner | |
374 | fRec = new AliPHOSReconstructioner(fClu, fTrs, fPID) ; | |
375 | ||
376 | //========== Event Number> | |
377 | if ( ( log10((Float_t)(ievent+1)) - (Int_t)(log10((Float_t)(ievent+1))) ) == 0. ) | |
378 | cout << "AnalyzeResolutions > " << "Event is " << ievent << endl ; | |
379 | ||
380 | //=========== Connects the various Tree's for evt | |
381 | gAlice->GetEvent(ievent); | |
382 | ||
383 | //=========== Get the Digit Tree | |
384 | gAlice->TreeD()->GetEvent(0) ; | |
385 | ||
386 | //=========== Do the reconstruction | |
387 | fPHOS->Reconstruction(fRec); | |
388 | ||
389 | //========== | |
390 | TClonesArray * RecParticleList = new TClonesArray("AliPHOSRecParticle", 2000) ; | |
391 | TBranch * RecBranch = gAlice->TreeR()->GetBranch("PHOSRP"); | |
392 | RecBranch-> SetAddress(&RecParticleList); | |
393 | ||
394 | //=========== Gets the Reconstraction TTree | |
395 | gAlice->TreeR()->GetEvent(0) ; | |
396 | ||
397 | //=========== Gets the Kine TTree | |
398 | gAlice->TreeK()->GetEvent(0) ; | |
399 | ||
400 | //=========== Gets the list of Primari Particles | |
401 | TClonesArray * PrimaryList = gAlice->Particles(); | |
402 | ||
403 | //========== TRackSegments in the event | |
404 | TIter nextRecParticle(RecParticleList) ; | |
405 | AliPHOSRecParticle * RecParticle ; | |
406 | ||
407 | while( (RecParticle = (AliPHOSRecParticle *) nextRecParticle())) | |
408 | { | |
409 | Int_t ModuleNumberRec ; | |
410 | Double_t RecX, RecZ ; | |
411 | fGeom->ImpactOnEmc(RecParticle->Theta(), RecParticle->Phi(), ModuleNumberRec, RecX, RecZ) ; | |
412 | ||
413 | Double_t MinDistance = 10000 ; | |
414 | Int_t ClosestPrimary = -1 ; | |
415 | ||
416 | Int_t numberofprimaries ; | |
417 | Int_t * listofprimaries = RecParticle->GetPrimaries(numberofprimaries) ; | |
418 | ||
419 | Int_t index ; | |
420 | TParticle * Primary ; | |
421 | Double_t Distance = MinDistance ; | |
422 | for ( index = 0 ; index < numberofprimaries ; index++) | |
423 | { | |
424 | Primary = (TParticle*)PrimaryList->UncheckedAt(listofprimaries[index]) ; | |
425 | ||
426 | Int_t ModuleNumber ; | |
427 | Double_t PrimX, PrimZ ; | |
428 | fGeom->ImpactOnEmc(Primary->Theta(), Primary->Phi(), ModuleNumber, PrimX, PrimZ) ; | |
429 | if(ModuleNumberRec == ModuleNumber) | |
430 | Distance = TMath::Sqrt((RecX-PrimX)*(RecX-PrimX)+(RecZ-PrimZ)*(RecZ-PrimZ) ) ; | |
431 | if(MinDistance > Distance) | |
432 | { | |
433 | MinDistance = Distance ; | |
434 | ClosestPrimary = listofprimaries[index] ; | |
435 | } | |
436 | } | |
437 | ||
438 | if(ClosestPrimary >=0 ) | |
439 | { | |
440 | switch(RecParticle->GetType()) | |
92862013 | 441 | { |
9ec91567 | 442 | case AliPHOSFastRecParticle::kGAMMA: |
134ce69a | 443 | fhPhotonEnergy->Fill(RecParticle->Energy(),((TParticle *) PrimaryList->At(ClosestPrimary))->Energy() ) ; |
444 | fhPhotonPosition->Fill(RecParticle->Energy(),Distance) ; | |
92862013 | 445 | break; |
9ec91567 | 446 | case AliPHOSFastRecParticle::kELECTRON: |
134ce69a | 447 | fhElectronEnergy->Fill(RecParticle->Energy(),((TParticle *)PrimaryList->At(ClosestPrimary))->Energy() ) ; |
448 | fhElectronPosition->Fill(RecParticle->Energy(),Distance ) ; | |
92862013 | 449 | break; |
9ec91567 | 450 | case AliPHOSFastRecParticle::kNEUTRALHA: |
134ce69a | 451 | fhNeutralHadronEnergy->Fill( RecParticle->Energy(),((TParticle *)PrimaryList->At(ClosestPrimary))->Energy()) ; |
452 | fhNeutralHadronPosition->Fill(RecParticle->Energy(),Distance ) ; | |
b9bbdad1 | 453 | break ; |
9ec91567 | 454 | case AliPHOSFastRecParticle::kNEUTRALEM: |
134ce69a | 455 | fhNeutralEMEnergy->Fill( RecParticle->Energy(),((TParticle *)PrimaryList->At(ClosestPrimary))->Energy() ) ; |
456 | fhNeutralEMPosition->Fill(RecParticle->Energy(),Distance ) ; | |
92862013 | 457 | break ; |
9ec91567 | 458 | case AliPHOSFastRecParticle::kCHARGEDHA: |
134ce69a | 459 | fhChargedHadronEnergy->Fill(RecParticle->Energy(),((TParticle *)PrimaryList->At(ClosestPrimary))->Energy() ) ; |
460 | fhChargedHadronPosition->Fill(RecParticle->Energy(),Distance ) ; | |
92862013 | 461 | break ; |
9ec91567 | 462 | case AliPHOSFastRecParticle::kGAMMAHA: |
134ce69a | 463 | fhPhotonHadronEnergy->Fill( RecParticle->Energy(),((TParticle *)PrimaryList->At(ClosestPrimary))->Energy()) ; |
464 | fhPhotonHadronPosition->Fill(RecParticle->Energy(),Distance ) ; | |
c1d256cb | 465 | break ; |
92862013 | 466 | } |
467 | } | |
468 | } | |
134ce69a | 469 | |
26d4b141 | 470 | // Deleting fClu, fTrs, fPID et fRec |
92862013 | 471 | fClu->Delete(); |
472 | fTrs->Delete(); | |
26d4b141 | 473 | fPID->Delete(); |
92862013 | 474 | fRec->Delete(); |
134ce69a | 475 | PrimaryList->Delete() ; |
92862013 | 476 | |
477 | } // endfor | |
134ce69a | 478 | SaveResolutionHistograms(); |
92862013 | 479 | } // endif |
480 | } // endfunction | |
481 | ||
482 | ||
483 | //____________________________________________________________________________ | |
484 | void AliPHOSAnalyze::BookingHistograms() | |
485 | { | |
b2a60966 | 486 | // Books the histograms where the results of the analysis are stored (to be changed) |
487 | ||
b9bbdad1 | 488 | if (fhEmcDigit ) |
489 | delete fhEmcDigit ; | |
490 | if (fhVetoDigit ) | |
491 | delete fhVetoDigit ; | |
492 | if (fhConvertorDigit ) | |
493 | delete fhConvertorDigit ; | |
494 | if (fhEmcCluster ) | |
495 | delete fhEmcCluster ; | |
496 | if (fhVetoCluster ) | |
497 | delete fhVetoCluster ; | |
498 | if (fhConvertorCluster ) | |
499 | delete fhConvertorCluster ; | |
500 | if (fhConvertorEmc ) | |
501 | delete fhConvertorEmc ; | |
502 | ||
134ce69a | 503 | // fhEmcDigit = new TH1F("hEmcDigit", "hEmcDigit", 1000, 0. , 25.); |
504 | // fhVetoDigit = new TH1F("hVetoDigit", "hVetoDigit", 500, 0. , 3.e-5); | |
505 | // fhConvertorDigit = new TH1F("hConvertorDigit","hConvertorDigit", 500, 0. , 3.e-5); | |
506 | // fhEmcCluster = new TH1F("hEmcCluster", "hEmcCluster", 1000, 0. , 30.); | |
507 | // fhVetoCluster = new TH1F("hVetoCluster", "hVetoCluster", 500, 0. , 3.e-5); | |
508 | // fhConvertorCluster = new TH1F("hConvertorCluster","hConvertorCluster",500, 0. , 3.e-5); | |
509 | // fhConvertorEmc = new TH2F("hConvertorEmc", "hConvertorEmc", 200, 1. , 3., 200, 0., 3.e-5); | |
510 | // fhPhotonEnergy = new TH1F("hPhotonEnergy", "hPhotonEnergy", 1000, 0. , 30.); | |
511 | // fhElectronEnergy = new TH1F("hElectronEnergy","hElectronEnergy", 1000, 0. , 30.); | |
512 | // fhNeutralHadronEnergy = new TH1F("hNeutralHadronEnergy", "hNeutralHadronEnergy", 1000, 0. , 30.); | |
513 | // fhNeutralEMEnergy = new TH1F("hNeutralEMEnergy", "hNeutralEMEnergy", 1000, 0. , 30.); | |
514 | // fhChargedHadronEnergy = new TH1F("hChargedHadronEnergy", "hChargedHadronEnergy", 1000, 0. , 30.); | |
515 | // fhPhotonHadronEnergy = new TH1F("hPhotonHadronEnergy","hPhotonHadronEnergy",500,-80. , 80.); | |
516 | // fhPhotonPositionX = new TH1F("hPhotonPositionX","hPhotonPositionX", 500,-80. , 80.); | |
517 | // fhElectronPositionX = new TH1F("hElectronPositionX","hElectronPositionX",500,-80. , 80.); | |
518 | // fhNeutralHadronPositionX = new TH1F("hNeutralHadronPositionX","hNeutralHadronPositionX",500,-80. , 80.); | |
519 | // fhNeutralEMPositionX = new TH1F("hNeutralEMPositionX","hNeutralEMPositionX",500,-80. , 80.); | |
520 | // fhChargedHadronPositionX = new TH1F("hChargedHadronPositionX","hChargedHadronPositionX",500,-80. , 80.); | |
521 | // fhPhotonHadronPositionX = new TH1F("hPhotonHadronPositionX","hPhotonHadronPositionX",500,-80. , 80.); | |
522 | // fhPhotonPositionY = new TH1F("hPhotonPositionY","hPhotonPositionY", 500,-80. , 80.); | |
523 | // fhElectronPositionY = new TH1F("hElectronPositionY","hElectronPositionY",500,-80. , 80.); | |
524 | // fhNeutralHadronPositionY = new TH1F("hNeutralHadronPositionY","hNeutralHadronPositionY",500,-80. , 80.); | |
525 | // fhNeutralEMPositionY = new TH1F("hNeutralEMPositionY","hNeutralEMPositionY",500,-80. , 80.); | |
526 | // fhChargedHadronPositionY = new TH1F("hChargedHadronPositionY","hChargedHadronPositionY",500,-80. , 80.); | |
527 | // fhPhotonHadronPositionY = new TH1F("hPhotonHadronPositionY","hPhotonHadronPositionY",500,-80. , 80.); | |
528 | ||
92862013 | 529 | |
134ce69a | 530 | } |
531 | //____________________________________________________________________________ | |
532 | void AliPHOSAnalyze::BookResolutionHistograms() | |
533 | { | |
534 | // Books the histograms where the results of the Resolution analysis are stored | |
535 | ||
536 | if(fhPhotonEnergy) | |
537 | delete fhPhotonEnergy ; | |
538 | if(fhElectronEnergy) | |
539 | delete fhElectronEnergy ; | |
540 | if(fhNeutralHadronEnergy) | |
541 | delete fhNeutralHadronEnergy ; | |
542 | if(fhNeutralEMEnergy) | |
543 | delete fhNeutralEMEnergy ; | |
544 | if(fhChargedHadronEnergy) | |
545 | delete fhChargedHadronEnergy ; | |
546 | if(fhPhotonHadronEnergy) | |
547 | delete fhPhotonHadronEnergy ; | |
548 | if(fhPhotonPosition) | |
549 | delete fhPhotonPosition ; | |
550 | if(fhElectronPosition) | |
551 | delete fhElectronPosition ; | |
552 | if(fhNeutralHadronPosition) | |
553 | delete fhNeutralHadronPosition ; | |
554 | if(fhNeutralEMPosition) | |
555 | delete fhNeutralEMPosition ; | |
556 | if(fhChargedHadronPosition) | |
557 | delete fhChargedHadronPosition ; | |
558 | if(fhPhotonHadronPosition) | |
559 | delete fhPhotonHadronPosition ; | |
560 | ||
561 | fhPhotonEnergy = new TH2F("hPhotonEnergy", "hPhotonEnergy", 100, 0., 5., 100, 0., 5.); | |
562 | fhElectronEnergy = new TH2F("hElectronEnergy","hElectronEnergy", 100, 0., 5., 100, 0., 5.); | |
563 | fhNeutralHadronEnergy = new TH2F("hNeutralHadronEnergy", "hNeutralHadronEnergy", 100, 0., 5., 100, 0., 5.); | |
564 | fhNeutralEMEnergy = new TH2F("hNeutralEMEnergy", "hNeutralEMEnergy", 100, 0., 5., 100, 0., 5.); | |
565 | fhChargedHadronEnergy = new TH2F("hChargedHadronEnergy", "hChargedHadronEnergy", 100, 0., 5., 100, 0., 5.); | |
566 | fhPhotonHadronEnergy = new TH2F("hPhotonHadronEnergy","hPhotonHadronEnergy", 100, 0., 5., 100, 0., 5.); | |
567 | fhPhotonPosition = new TH2F("hPhotonPosition","hPhotonPosition", 100, 0., 5., 100, 0., 5.); | |
568 | fhElectronPosition = new TH2F("hElectronPosition","hElectronPosition", 100, 0., 5., 100, 0., 5.); | |
569 | fhNeutralHadronPosition = new TH2F("hNeutralHadronPosition","hNeutralHadronPosition", 100, 0., 5., 100, 0., 5.); | |
570 | fhNeutralEMPosition = new TH2F("hNeutralEMPosition","hNeutralEMPosition", 100, 0., 5., 100, 0., 5.); | |
571 | fhChargedHadronPosition = new TH2F("hChargedHadronPosition","hChargedHadronPosition", 100, 0., 5., 100, 0., 5.); | |
572 | fhPhotonHadronPosition = new TH2F("hPhotonHadronPosition","hPhotonHadronPosition", 100, 0., 5., 100, 0., 5.); | |
3fc07a80 | 573 | |
92862013 | 574 | } |
6ad0bfa0 | 575 | //____________________________________________________________________________ |
576 | Bool_t AliPHOSAnalyze::Init(Int_t evt) | |
577 | { | |
b2a60966 | 578 | // Do a few initializations: open the root file |
579 | // get the AliRun object | |
580 | // defines the clusterizer, tracksegment maker and particle identifier | |
581 | // sets the associated parameters | |
6ad0bfa0 | 582 | |
92862013 | 583 | Bool_t ok = kTRUE ; |
6ad0bfa0 | 584 | |
585 | //========== Open galice root file | |
586 | ||
587 | if ( fRootFile == 0 ) { | |
588 | Text_t * name = new Text_t[80] ; | |
589 | cout << "AnalyzeOneEvent > Enter file root file name : " ; | |
590 | cin >> name ; | |
92862013 | 591 | Bool_t ok = OpenRootFile(name) ; |
592 | if ( !ok ) | |
6ad0bfa0 | 593 | cout << " AliPHOSAnalyze > Error opening " << name << endl ; |
594 | else { | |
595 | //========== Get AliRun object from file | |
596 | ||
597 | gAlice = (AliRun*) fRootFile->Get("gAlice") ; | |
598 | ||
599 | //=========== Get the PHOS object and associated geometry from the file | |
600 | ||
5f20d3fb | 601 | fPHOS = (AliPHOSv1 *)gAlice->GetDetector("PHOS") ; |
aafe457d | 602 | fGeom = fPHOS->GetGeometry(); |
603 | // fGeom = AliPHOSGeometry::GetInstance( fPHOS->GetGeometry()->GetName(), fPHOS->GetGeometry()->GetTitle() ); | |
83974468 | 604 | |
92862013 | 605 | } // else !ok |
6ad0bfa0 | 606 | } // if fRootFile |
607 | ||
92862013 | 608 | if ( ok ) { |
6ad0bfa0 | 609 | |
83974468 | 610 | |
611 | //========== Initializes the Index to Object converter | |
612 | ||
613 | fObjGetter = AliPHOSIndexToObject::GetInstance(fPHOS) ; | |
614 | ||
6ad0bfa0 | 615 | //========== Create the Clusterizer |
616 | ||
774e78c2 | 617 | fClu = new AliPHOSClusterizerv1() ; |
618 | fClu->SetEmcEnergyThreshold(0.030) ; | |
aafe457d | 619 | fClu->SetEmcClusteringThreshold(0.20) ; |
92862013 | 620 | fClu->SetPpsdEnergyThreshold (0.0000002) ; |
621 | fClu->SetPpsdClusteringThreshold(0.0000001) ; | |
6ad0bfa0 | 622 | fClu->SetLocalMaxCut(0.03) ; |
92862013 | 623 | fClu->SetCalibrationParameters(0., 0.00000001) ; |
6ad0bfa0 | 624 | cout << "AnalyzeOneEvent > using clusterizer " << fClu->GetName() << endl ; |
625 | fClu->PrintParameters() ; | |
626 | ||
627 | //========== Creates the track segment maker | |
628 | ||
629 | fTrs = new AliPHOSTrackSegmentMakerv1() ; | |
630 | cout << "AnalyzeOneEvent > using tack segment maker " << fTrs->GetName() << endl ; | |
134ce69a | 631 | // fTrs->UnsetUnfoldFlag() ; |
6ad0bfa0 | 632 | |
26d4b141 | 633 | //========== Creates the particle identifier |
6ad0bfa0 | 634 | |
26d4b141 | 635 | fPID = new AliPHOSPIDv1() ; |
636 | cout << "AnalyzeOneEvent > using particle identifier " << fPID->GetName() << endl ; | |
2aad621e | 637 | //fPID->SetShowerProfileCuts(Float_t l1m, Float_t l1M, Float_t l2m, Float_t l2M) ; |
638 | fPID->SetShowerProfileCuts(0.7, 2.0 , 0.6 , 1.5) ; | |
639 | ||
6ad0bfa0 | 640 | //========== Creates the Reconstructioner |
641 | ||
2aad621e | 642 | fRec = new AliPHOSReconstructioner(fClu, fTrs, fPID) ; |
134ce69a | 643 | fRec -> SetDebugReconstruction(kFALSE); |
6ad0bfa0 | 644 | |
645 | //=========== Connect the various Tree's for evt | |
646 | ||
647 | if ( evt == -999 ) { | |
648 | cout << "AnalyzeOneEvent > Enter event number : " ; | |
649 | cin >> evt ; | |
650 | cout << evt << endl ; | |
651 | } | |
652 | fEvt = evt ; | |
653 | gAlice->GetEvent(evt); | |
654 | ||
655 | //=========== Get the Digit TTree | |
656 | ||
657 | gAlice->TreeD()->GetEvent(0) ; | |
658 | ||
92862013 | 659 | } // ok |
6ad0bfa0 | 660 | |
92862013 | 661 | return ok ; |
6ad0bfa0 | 662 | } |
663 | ||
92862013 | 664 | |
6ad0bfa0 | 665 | //____________________________________________________________________________ |
92862013 | 666 | void AliPHOSAnalyze::DisplayKineEvent(Int_t evt) |
6ad0bfa0 | 667 | { |
b2a60966 | 668 | // Display particles from the Kine Tree in global Alice (theta, phi) coordinates. |
669 | // One PHOS module at the time. | |
670 | // The particle type can be selected. | |
671 | ||
6ad0bfa0 | 672 | if (evt == -999) |
673 | evt = fEvt ; | |
674 | ||
92862013 | 675 | Int_t module ; |
6ad0bfa0 | 676 | cout << "DisplayKineEvent > which module (1-5, -1: all) ? " ; |
92862013 | 677 | cin >> module ; cout << module << endl ; |
6ad0bfa0 | 678 | |
92862013 | 679 | Int_t testparticle ; |
6ad0bfa0 | 680 | cout << " 22 : PHOTON " << endl |
681 | << " (-)11 : (POSITRON)ELECTRON " << endl | |
682 | << " (-)2112 : (ANTI)NEUTRON " << endl | |
683 | << " -999 : Everything else " << endl ; | |
684 | cout << "DisplayKineEvent > enter PDG particle code to display " ; | |
92862013 | 685 | cin >> testparticle ; cout << testparticle << endl ; |
6ad0bfa0 | 686 | |
92862013 | 687 | Text_t histoname[80] ; |
688 | sprintf(histoname,"Event %d: Incident particles in module %d", evt, module) ; | |
6ad0bfa0 | 689 | |
92862013 | 690 | Double_t tm, tM, pm, pM ; // min and Max theta and phi covered by module |
cf0c2bc1 | 691 | fGeom->EmcModuleCoverage(module, tm, tM, pm, pM, AliPHOSGeometry::Degre() ) ; |
6ad0bfa0 | 692 | |
693 | Double_t theta, phi ; | |
cf0c2bc1 | 694 | fGeom->EmcXtalCoverage(theta, phi, AliPHOSGeometry::Degre() ) ; |
6ad0bfa0 | 695 | |
696 | Int_t tdim = (Int_t)( (tM - tm) / theta ) ; | |
697 | Int_t pdim = (Int_t)( (pM - pm) / phi ) ; | |
698 | ||
699 | tm -= theta ; | |
700 | tM += theta ; | |
701 | pm -= phi ; | |
702 | pM += phi ; | |
703 | ||
92862013 | 704 | TH2F * histoparticle = new TH2F("histoparticle", histoname, |
6ad0bfa0 | 705 | pdim, pm, pM, tdim, tm, tM) ; |
92862013 | 706 | histoparticle->SetStats(kFALSE) ; |
6ad0bfa0 | 707 | |
708 | // Get pointers to Alice Particle TClonesArray | |
709 | ||
92862013 | 710 | TParticle * particle; |
711 | TClonesArray * particlearray = gAlice->Particles(); | |
6ad0bfa0 | 712 | |
92862013 | 713 | Text_t canvasname[80]; |
714 | sprintf(canvasname,"Particles incident in PHOS/EMC module # %d",module) ; | |
715 | TCanvas * kinecanvas = new TCanvas("kinecanvas", canvasname, 650, 500) ; | |
6ad0bfa0 | 716 | |
717 | // get the KINE Tree | |
718 | ||
92862013 | 719 | TTree * kine = gAlice->TreeK() ; |
720 | Stat_t nParticles = kine->GetEntries() ; | |
721 | cout << "DisplayKineEvent > events in kine " << nParticles << endl ; | |
6ad0bfa0 | 722 | |
723 | // loop over particles | |
724 | ||
92862013 | 725 | Double_t kRADDEG = 180. / TMath::Pi() ; |
6ad0bfa0 | 726 | Int_t index1 ; |
727 | Int_t nparticlein = 0 ; | |
92862013 | 728 | for (index1 = 0 ; index1 < nParticles ; index1++){ |
729 | Int_t nparticle = particlearray->GetEntriesFast() ; | |
6ad0bfa0 | 730 | Int_t index2 ; |
731 | for( index2 = 0 ; index2 < nparticle ; index2++) { | |
92862013 | 732 | particle = (TParticle*)particlearray->UncheckedAt(index2) ; |
733 | Int_t particletype = particle->GetPdgCode() ; | |
734 | if (testparticle == -999 || testparticle == particletype) { | |
735 | Double_t phi = particle->Phi() ; | |
736 | Double_t theta = particle->Theta() ; | |
6ad0bfa0 | 737 | Int_t mod ; |
738 | Double_t x, z ; | |
92862013 | 739 | fGeom->ImpactOnEmc(theta, phi, mod, z, x) ; |
740 | if ( mod == module ) { | |
6ad0bfa0 | 741 | nparticlein++ ; |
b2a60966 | 742 | if (particle->Energy() > fClu->GetEmcClusteringThreshold() ) |
743 | histoparticle->Fill(phi*kRADDEG, theta*kRADDEG, particle->Energy() ) ; | |
6ad0bfa0 | 744 | } |
745 | } | |
746 | } | |
747 | } | |
92862013 | 748 | kinecanvas->Draw() ; |
749 | histoparticle->Draw("color") ; | |
750 | TPaveText * pavetext = new TPaveText(294, 100, 300, 101); | |
6ad0bfa0 | 751 | Text_t text[40] ; |
752 | sprintf(text, "Particles: %d ", nparticlein) ; | |
92862013 | 753 | pavetext->AddText(text) ; |
754 | pavetext->Draw() ; | |
755 | kinecanvas->Update(); | |
6ad0bfa0 | 756 | |
757 | } | |
758 | //____________________________________________________________________________ | |
759 | void AliPHOSAnalyze::DisplayRecParticles() | |
760 | { | |
b2a60966 | 761 | // Display reconstructed particles in global Alice(theta, phi) coordinates. |
762 | // One PHOS module at the time. | |
763 | // Click on symbols indicate the reconstructed particle type. | |
764 | ||
6ad0bfa0 | 765 | if (fEvt == -999) { |
15605d3c | 766 | cout << "DisplayRecParticles > Analyze an event first ... (y/n) " ; |
6ad0bfa0 | 767 | Text_t answer[1] ; |
768 | cin >> answer ; cout << answer ; | |
769 | if ( answer == "y" ) | |
770 | AnalyzeOneEvent() ; | |
771 | } | |
772 | if (fEvt != -999) { | |
773 | ||
92862013 | 774 | Int_t module ; |
15605d3c | 775 | cout << "DisplayRecParticles > which module (1-5, -1: all) ? " ; |
92862013 | 776 | cin >> module ; cout << module << endl ; |
777 | Text_t histoname[80] ; | |
778 | sprintf(histoname,"Event %d: Reconstructed particles in module %d", fEvt, module) ; | |
779 | Double_t tm, tM, pm, pM ; // min and Max theta and phi covered by module | |
cf0c2bc1 | 780 | fGeom->EmcModuleCoverage(module, tm, tM, pm, pM, AliPHOSGeometry::Degre() ) ; |
6ad0bfa0 | 781 | Double_t theta, phi ; |
cf0c2bc1 | 782 | fGeom->EmcXtalCoverage(theta, phi, AliPHOSGeometry::Degre() ) ; |
6ad0bfa0 | 783 | Int_t tdim = (Int_t)( (tM - tm) / theta ) ; |
784 | Int_t pdim = (Int_t)( (pM - pm) / phi ) ; | |
785 | tm -= theta ; | |
786 | tM += theta ; | |
787 | pm -= phi ; | |
92862013 | 788 | TH2F * histoRparticle = new TH2F("histoRparticle", histoname, |
6ad0bfa0 | 789 | pdim, pm, pM, tdim, tm, tM) ; |
92862013 | 790 | histoRparticle->SetStats(kFALSE) ; |
791 | Text_t canvasname[80] ; | |
792 | sprintf(canvasname, "Reconstructed particles in PHOSmodule # %d", module) ; | |
793 | TCanvas * rparticlecanvas = new TCanvas("RparticleCanvas", canvasname, 650, 500) ; | |
88714635 | 794 | AliPHOSRecParticle::RecParticlesList * rpl = fPHOS->RecParticles() ; |
92862013 | 795 | Int_t nRecParticles = rpl->GetEntries() ; |
796 | Int_t nRecParticlesInModule = 0 ; | |
6ad0bfa0 | 797 | TIter nextRecPart(rpl) ; |
798 | AliPHOSRecParticle * rp ; | |
92862013 | 799 | cout << "DisplayRecParticles > " << nRecParticles << " reconstructed particles " << endl ; |
800 | Double_t kRADDEG = 180. / TMath::Pi() ; | |
6ad0bfa0 | 801 | while ( (rp = (AliPHOSRecParticle *)nextRecPart() ) ) { |
802 | AliPHOSTrackSegment * ts = rp->GetPHOSTrackSegment() ; | |
b2a60966 | 803 | if ( ts->GetPHOSMod() == module ) { |
804 | Int_t numberofprimaries = 0 ; | |
134ce69a | 805 | Int_t * listofprimaries = 0; |
806 | rp->GetPrimaries(numberofprimaries) ; | |
b2a60966 | 807 | cout << "Number of primaries = " << numberofprimaries << endl ; |
808 | Int_t index ; | |
809 | for ( index = 0 ; index < numberofprimaries ; index++) | |
810 | cout << " primary # " << index << " = " << listofprimaries[index] << endl ; | |
811 | ||
92862013 | 812 | nRecParticlesInModule++ ; |
813 | Double_t theta = rp->Theta() * kRADDEG ; | |
814 | Double_t phi = rp->Phi() * kRADDEG ; | |
6ad0bfa0 | 815 | Double_t energy = rp->Energy() ; |
92862013 | 816 | histoRparticle->Fill(phi, theta, energy) ; |
6ad0bfa0 | 817 | } |
818 | } | |
92862013 | 819 | histoRparticle->Draw("color") ; |
15605d3c | 820 | |
821 | nextRecPart.Reset() ; | |
822 | while ( (rp = (AliPHOSRecParticle *)nextRecPart() ) ) { | |
823 | AliPHOSTrackSegment * ts = rp->GetPHOSTrackSegment() ; | |
824 | if ( ts->GetPHOSMod() == module ) | |
825 | rp->Draw("P") ; | |
826 | } | |
827 | ||
6ad0bfa0 | 828 | Text_t text[80] ; |
92862013 | 829 | sprintf(text, "reconstructed particles: %d", nRecParticlesInModule) ; |
830 | TPaveText * pavetext = new TPaveText(292, 100, 300, 101); | |
831 | pavetext->AddText(text) ; | |
832 | pavetext->Draw() ; | |
833 | rparticlecanvas->Update() ; | |
6ad0bfa0 | 834 | } |
835 | } | |
836 | ||
837 | //____________________________________________________________________________ | |
838 | void AliPHOSAnalyze::DisplayRecPoints() | |
839 | { | |
b2a60966 | 840 | // Display reconstructed points in local PHOS-module (x, z) coordinates. |
841 | // One PHOS module at the time. | |
842 | // Click on symbols displays the EMC cluster, or PPSD information. | |
843 | ||
6ad0bfa0 | 844 | if (fEvt == -999) { |
845 | cout << "DisplayRecPoints > Analyze an event first ... (y/n) " ; | |
846 | Text_t answer[1] ; | |
847 | cin >> answer ; cout << answer ; | |
848 | if ( answer == "y" ) | |
849 | AnalyzeOneEvent() ; | |
850 | } | |
851 | if (fEvt != -999) { | |
852 | ||
92862013 | 853 | Int_t module ; |
6ad0bfa0 | 854 | cout << "DisplayRecPoints > which module (1-5, -1: all) ? " ; |
92862013 | 855 | cin >> module ; cout << module << endl ; |
6ad0bfa0 | 856 | |
92862013 | 857 | Text_t canvasname[80]; |
858 | sprintf(canvasname,"Digits in PHOS/EMC module # %d",module) ; | |
859 | TCanvas * modulecanvas = new TCanvas("module", canvasname, 650, 500) ; | |
860 | modulecanvas->Draw() ; | |
6ad0bfa0 | 861 | |
92862013 | 862 | //=========== Creating 2d-histogram of the PHOS module |
6ad0bfa0 | 863 | // a little bit junkie but is used to test Geom functinalities |
864 | ||
92862013 | 865 | Double_t tm, tM, pm, pM ; // min and Max theta and phi covered by module |
6ad0bfa0 | 866 | |
92862013 | 867 | fGeom->EmcModuleCoverage(module, tm, tM, pm, pM); |
6ad0bfa0 | 868 | // convert angles into coordinates local to the EMC module of interest |
869 | ||
92862013 | 870 | Int_t emcModuleNumber ; |
871 | Double_t emcModulexm, emcModulezm ; // minimum local coordinate in a given EMCA module | |
872 | Double_t emcModulexM, emcModulezM ; // maximum local coordinate in a given EMCA module | |
873 | fGeom->ImpactOnEmc(tm, pm, emcModuleNumber, emcModulezm, emcModulexm) ; | |
874 | fGeom->ImpactOnEmc(tM, pM, emcModuleNumber, emcModulezM, emcModulexM) ; | |
875 | Int_t xdim = (Int_t)( ( emcModulexM - emcModulexm ) / fGeom->GetCrystalSize(0) ) ; | |
876 | Int_t zdim = (Int_t)( ( emcModulezM - emcModulezm ) / fGeom->GetCrystalSize(2) ) ; | |
877 | Float_t xmin = emcModulexm - fGeom->GetCrystalSize(0) ; | |
878 | Float_t xMax = emcModulexM + fGeom->GetCrystalSize(0) ; | |
879 | Float_t zmin = emcModulezm - fGeom->GetCrystalSize(2) ; | |
880 | Float_t zMax = emcModulezM + fGeom->GetCrystalSize(2) ; | |
881 | Text_t histoname[80]; | |
882 | sprintf(histoname,"Event %d: Digits and RecPoints in module %d", fEvt, module) ; | |
883 | TH2F * hModule = new TH2F("HistoReconstructed", histoname, | |
6ad0bfa0 | 884 | xdim, xmin, xMax, zdim, zmin, zMax) ; |
885 | hModule->SetMaximum(2.0); | |
886 | hModule->SetMinimum(0.0); | |
887 | hModule->SetStats(kFALSE); | |
888 | ||
889 | TIter next(fPHOS->Digits()) ; | |
92862013 | 890 | Float_t energy, y, z; |
891 | Float_t etot=0.; | |
892 | Int_t relid[4]; Int_t nDigits = 0 ; | |
6ad0bfa0 | 893 | AliPHOSDigit * digit ; |
78ccc411 | 894 | |
895 | // Making 2D histogram of the EMC module | |
6ad0bfa0 | 896 | while((digit = (AliPHOSDigit *)next())) |
897 | { | |
92862013 | 898 | fGeom->AbsToRelNumbering(digit->GetId(), relid) ; |
78ccc411 | 899 | if (relid[0] == module && relid[1] == 0) |
5959e315 | 900 | { |
92862013 | 901 | energy = fClu->Calibrate(digit->GetAmp()) ; |
774e78c2 | 902 | cout << "Energy is " << energy << " and threshold is " << fClu->GetEmcEnergyThreshold() << endl; |
78ccc411 | 903 | if (energy > fClu->GetEmcEnergyThreshold() ){ |
5959e315 | 904 | nDigits++ ; |
905 | etot += energy ; | |
906 | fGeom->RelPosInModule(relid,y,z) ; | |
907 | hModule->Fill(y, z, energy) ; | |
3fc07a80 | 908 | } |
6ad0bfa0 | 909 | } |
910 | } | |
92862013 | 911 | cout <<"DrawRecPoints > Found in module " |
912 | << module << " " << nDigits << " digits with total energy " << etot << endl ; | |
6ad0bfa0 | 913 | hModule->Draw("col2") ; |
914 | ||
92862013 | 915 | //=========== Cluster in module |
6ad0bfa0 | 916 | |
83974468 | 917 | // TClonesArray * emcRP = fPHOS->EmcClusters() ; |
918 | TObjArray * emcRP = fPHOS->EmcRecPoints() ; | |
919 | ||
92862013 | 920 | etot = 0.; |
921 | Int_t totalnClusters = 0 ; | |
922 | Int_t nClusters = 0 ; | |
923 | TIter nextemc(emcRP) ; | |
6ad0bfa0 | 924 | AliPHOSEmcRecPoint * emc ; |
925 | while((emc = (AliPHOSEmcRecPoint *)nextemc())) | |
926 | { | |
3fc07a80 | 927 | // Int_t numberofprimaries ; |
928 | // Int_t * primariesarray = new Int_t[10] ; | |
929 | // emc->GetPrimaries(numberofprimaries, primariesarray) ; | |
92862013 | 930 | totalnClusters++ ; |
931 | if ( emc->GetPHOSMod() == module ) | |
6ad0bfa0 | 932 | { |
92862013 | 933 | nClusters++ ; |
934 | energy = emc->GetTotalEnergy() ; | |
935 | etot+= energy ; | |
26d4b141 | 936 | emc->Draw("M") ; |
6ad0bfa0 | 937 | } |
938 | } | |
92862013 | 939 | cout << "DrawRecPoints > Found " << totalnClusters << " EMC Clusters in PHOS" << endl ; |
940 | cout << "DrawRecPoints > Found in module " << module << " " << nClusters << " EMC Clusters " << endl ; | |
941 | cout << "DrawRecPoints > total energy " << etot << endl ; | |
6ad0bfa0 | 942 | |
92862013 | 943 | TPaveText * pavetext = new TPaveText(22, 80, 83, 90); |
6ad0bfa0 | 944 | Text_t text[40] ; |
92862013 | 945 | sprintf(text, "digits: %d; clusters: %d", nDigits, nClusters) ; |
946 | pavetext->AddText(text) ; | |
947 | pavetext->Draw() ; | |
948 | modulecanvas->Update(); | |
6ad0bfa0 | 949 | |
92862013 | 950 | //=========== Cluster in module PPSD Down |
6ad0bfa0 | 951 | |
83974468 | 952 | // TClonesArray * ppsdRP = fPHOS->PpsdClusters() ; |
953 | TObjArray * ppsdRP = fPHOS->PpsdRecPoints() ; | |
954 | ||
92862013 | 955 | etot = 0.; |
956 | TIter nextPpsd(ppsdRP) ; | |
957 | AliPHOSPpsdRecPoint * ppsd ; | |
958 | while((ppsd = (AliPHOSPpsdRecPoint *)nextPpsd())) | |
6ad0bfa0 | 959 | { |
92862013 | 960 | totalnClusters++ ; |
961 | if ( ppsd->GetPHOSMod() == module ) | |
6ad0bfa0 | 962 | { |
92862013 | 963 | nClusters++ ; |
964 | energy = ppsd->GetEnergy() ; | |
965 | etot+=energy ; | |
966 | if (!ppsd->GetUp()) ppsd->Draw("P") ; | |
6ad0bfa0 | 967 | } |
968 | } | |
92862013 | 969 | cout << "DrawRecPoints > Found " << totalnClusters << " Ppsd Down Clusters in PHOS" << endl ; |
970 | cout << "DrawRecPoints > Found in module " << module << " " << nClusters << " Ppsd Down Clusters " << endl ; | |
971 | cout << "DrawRecPoints > total energy " << etot << endl ; | |
6ad0bfa0 | 972 | |
92862013 | 973 | //=========== Cluster in module PPSD Up |
6ad0bfa0 | 974 | |
83974468 | 975 | ppsdRP = fPHOS->PpsdRecPoints() ; |
976 | ||
92862013 | 977 | etot = 0.; |
978 | TIter nextPpsdUp(ppsdRP) ; | |
979 | while((ppsd = (AliPHOSPpsdRecPoint *)nextPpsdUp())) | |
6ad0bfa0 | 980 | { |
92862013 | 981 | totalnClusters++ ; |
982 | if ( ppsd->GetPHOSMod() == module ) | |
6ad0bfa0 | 983 | { |
92862013 | 984 | nClusters++ ; |
985 | energy = ppsd->GetEnergy() ; | |
986 | etot+=energy ; | |
987 | if (ppsd->GetUp()) ppsd->Draw("P") ; | |
6ad0bfa0 | 988 | } |
989 | } | |
92862013 | 990 | cout << "DrawRecPoints > Found " << totalnClusters << " Ppsd Up Clusters in PHOS" << endl ; |
991 | cout << "DrawRecPoints > Found in module " << module << " " << nClusters << " Ppsd Up Clusters " << endl ; | |
992 | cout << "DrawRecPoints > total energy " << etot << endl ; | |
6ad0bfa0 | 993 | |
994 | } // if !-999 | |
995 | } | |
996 | ||
997 | //____________________________________________________________________________ | |
998 | void AliPHOSAnalyze::DisplayTrackSegments() | |
999 | { | |
b2a60966 | 1000 | // Display track segments in local PHOS-module (x, z) coordinates. |
1001 | // One PHOS module at the time. | |
1002 | // One symbol per PHOS subsystem: EMC, upper PPSD, lower PPSD. | |
1003 | ||
6ad0bfa0 | 1004 | if (fEvt == -999) { |
1005 | cout << "DisplayTrackSegments > Analyze an event first ... (y/n) " ; | |
1006 | Text_t answer[1] ; | |
1007 | cin >> answer ; cout << answer ; | |
1008 | if ( answer == "y" ) | |
1009 | AnalyzeOneEvent() ; | |
1010 | } | |
1011 | if (fEvt != -999) { | |
1012 | ||
92862013 | 1013 | Int_t module ; |
6ad0bfa0 | 1014 | cout << "DisplayTrackSegments > which module (1-5, -1: all) ? " ; |
92862013 | 1015 | cin >> module ; cout << module << endl ; |
1016 | //=========== Creating 2d-histogram of the PHOS module | |
6ad0bfa0 | 1017 | // a little bit junkie but is used to test Geom functinalities |
1018 | ||
92862013 | 1019 | Double_t tm, tM, pm, pM ; // min and Max theta and phi covered by module |
6ad0bfa0 | 1020 | |
92862013 | 1021 | fGeom->EmcModuleCoverage(module, tm, tM, pm, pM); |
6ad0bfa0 | 1022 | // convert angles into coordinates local to the EMC module of interest |
1023 | ||
92862013 | 1024 | Int_t emcModuleNumber ; |
1025 | Double_t emcModulexm, emcModulezm ; // minimum local coordinate in a given EMCA module | |
1026 | Double_t emcModulexM, emcModulezM ; // maximum local coordinate in a given EMCA module | |
1027 | fGeom->ImpactOnEmc(tm, pm, emcModuleNumber, emcModulezm, emcModulexm) ; | |
1028 | fGeom->ImpactOnEmc(tM, pM, emcModuleNumber, emcModulezM, emcModulexM) ; | |
1029 | Int_t xdim = (Int_t)( ( emcModulexM - emcModulexm ) / fGeom->GetCrystalSize(0) ) ; | |
1030 | Int_t zdim = (Int_t)( ( emcModulezM - emcModulezm ) / fGeom->GetCrystalSize(2) ) ; | |
1031 | Float_t xmin = emcModulexm - fGeom->GetCrystalSize(0) ; | |
1032 | Float_t xMax = emcModulexM + fGeom->GetCrystalSize(0) ; | |
1033 | Float_t zmin = emcModulezm - fGeom->GetCrystalSize(2) ; | |
1034 | Float_t zMax = emcModulezM + fGeom->GetCrystalSize(2) ; | |
1035 | Text_t histoname[80]; | |
1036 | sprintf(histoname,"Event %d: Track Segments in module %d", fEvt, module) ; | |
1037 | TH2F * histotrack = new TH2F("histotrack", histoname, | |
6ad0bfa0 | 1038 | xdim, xmin, xMax, zdim, zmin, zMax) ; |
92862013 | 1039 | histotrack->SetStats(kFALSE); |
1040 | Text_t canvasname[80]; | |
1041 | sprintf(canvasname,"Track segments in PHOS/EMC-PPSD module # %d", module) ; | |
1042 | TCanvas * trackcanvas = new TCanvas("TrackSegmentCanvas", canvasname, 650, 500) ; | |
1043 | histotrack->Draw() ; | |
6ad0bfa0 | 1044 | |
88714635 | 1045 | AliPHOSTrackSegment::TrackSegmentsList * trsegl = fPHOS->TrackSegments() ; |
6ad0bfa0 | 1046 | AliPHOSTrackSegment * trseg ; |
1047 | ||
92862013 | 1048 | Int_t nTrackSegments = trsegl->GetEntries() ; |
6ad0bfa0 | 1049 | Int_t index ; |
92862013 | 1050 | Float_t etot = 0 ; |
1051 | Int_t nTrackSegmentsInModule = 0 ; | |
1052 | for(index = 0; index < nTrackSegments ; index++){ | |
6ad0bfa0 | 1053 | trseg = (AliPHOSTrackSegment * )trsegl->At(index) ; |
92862013 | 1054 | etot+= trseg->GetEnergy() ; |
1055 | if ( trseg->GetPHOSMod() == module ) { | |
1056 | nTrackSegmentsInModule++ ; | |
6ad0bfa0 | 1057 | trseg->Draw("P"); |
1058 | } | |
1059 | } | |
1060 | Text_t text[80] ; | |
92862013 | 1061 | sprintf(text, "track segments: %d", nTrackSegmentsInModule) ; |
1062 | TPaveText * pavetext = new TPaveText(22, 80, 83, 90); | |
1063 | pavetext->AddText(text) ; | |
1064 | pavetext->Draw() ; | |
1065 | trackcanvas->Update() ; | |
1066 | cout << "DisplayTrackSegments > Found " << trsegl->GetEntries() << " Track segments with total energy "<< etot << endl ; | |
6ad0bfa0 | 1067 | |
1068 | } | |
1069 | } | |
1070 | //____________________________________________________________________________ | |
1071 | Bool_t AliPHOSAnalyze::OpenRootFile(Text_t * name) | |
1072 | { | |
b2a60966 | 1073 | // Open the root file named "name" |
1074 | ||
1075 | fRootFile = new TFile(name, "update") ; | |
6ad0bfa0 | 1076 | return fRootFile->IsOpen() ; |
1077 | } | |
92862013 | 1078 | //____________________________________________________________________________ |
1079 | void AliPHOSAnalyze::SavingHistograms() | |
1080 | { | |
b2a60966 | 1081 | // Saves the histograms in a root file named "name.analyzed" |
1082 | ||
1083 | Text_t outputname[80] ; | |
3862b681 | 1084 | sprintf(outputname,"%s.analyzed",fRootFile->GetName()); |
1085 | TFile output(outputname,"RECREATE"); | |
92862013 | 1086 | output.cd(); |
134ce69a | 1087 | // if (fhEmcDigit ) |
1088 | // fhEmcDigit->Write() ; | |
1089 | // if (fhVetoDigit ) | |
1090 | // fhVetoDigit->Write() ; | |
1091 | // if (fhConvertorDigit ) | |
1092 | // fhConvertorDigit->Write() ; | |
1093 | // if (fhEmcCluster ) | |
1094 | // fhEmcCluster->Write() ; | |
1095 | // if (fhVetoCluster ) | |
1096 | // fhVetoCluster->Write() ; | |
1097 | // if (fhConvertorCluster ) | |
1098 | // fhConvertorCluster->Write() ; | |
1099 | // if (fhConvertorEmc ) | |
1100 | // fhConvertorEmc->Write() ; | |
1101 | // if (fhPhotonEnergy) | |
1102 | // fhPhotonEnergy->Write() ; | |
1103 | // if (fhPhotonPositionX) | |
1104 | // fhPhotonPositionX->Write() ; | |
1105 | // if (fhPhotonPositionY) | |
1106 | // fhPhotonPositionX->Write() ; | |
1107 | // if (fhElectronEnergy) | |
1108 | // fhElectronEnergy->Write() ; | |
1109 | // if (fhElectronPositionX) | |
1110 | // fhElectronPositionX->Write() ; | |
1111 | // if (fhElectronPositionY) | |
1112 | // fhElectronPositionX->Write() ; | |
1113 | // if (fhNeutralHadronEnergy) | |
1114 | // fhNeutralHadronEnergy->Write() ; | |
1115 | // if (fhNeutralHadronPositionX) | |
1116 | // fhNeutralHadronPositionX->Write() ; | |
1117 | // if (fhNeutralHadronPositionY) | |
1118 | // fhNeutralHadronPositionX->Write() ; | |
1119 | // if (fhNeutralEMEnergy) | |
1120 | // fhNeutralEMEnergy->Write() ; | |
1121 | // if (fhNeutralEMPositionX) | |
1122 | // fhNeutralEMPositionX->Write() ; | |
1123 | // if (fhNeutralEMPositionY) | |
1124 | // fhNeutralEMPositionX->Write() ; | |
1125 | // if (fhChargedHadronEnergy) | |
1126 | // fhChargedHadronEnergy->Write() ; | |
1127 | // if (fhChargedHadronPositionX) | |
1128 | // fhChargedHadronPositionX->Write() ; | |
1129 | // if (fhChargedHadronPositionY) | |
1130 | // fhChargedHadronPositionX->Write() ; | |
1131 | // if (fhPhotonHadronEnergy) | |
1132 | // fhPhotonHadronEnergy->Write() ; | |
1133 | // if (fhPhotonHadronPositionX) | |
1134 | // fhPhotonHadronPositionX->Write() ; | |
1135 | // if (fhPhotonHadronPositionY) | |
1136 | // fhPhotonHadronPositionX->Write() ; | |
1137 | ||
1138 | output.Write(); | |
1139 | output.Close(); | |
1140 | } | |
1141 | //____________________________________________________________________________ | |
1142 | void AliPHOSAnalyze::SaveResolutionHistograms() | |
1143 | { | |
1144 | // Saves the histograms in a root file named "name.analyzed" | |
1145 | ||
1146 | Text_t outputname[80] ; | |
1147 | sprintf(outputname,"%s.analyzed",fRootFile->GetName()); | |
1148 | TFile output(outputname,"RECREATE"); | |
1149 | output.cd(); | |
b9bbdad1 | 1150 | if (fhPhotonEnergy) |
1151 | fhPhotonEnergy->Write() ; | |
134ce69a | 1152 | if (fhPhotonPosition) |
1153 | fhPhotonPosition->Write() ; | |
b9bbdad1 | 1154 | if (fhElectronEnergy) |
1155 | fhElectronEnergy->Write() ; | |
134ce69a | 1156 | if (fhElectronPosition) |
1157 | fhElectronPosition->Write() ; | |
b9bbdad1 | 1158 | if (fhNeutralHadronEnergy) |
1159 | fhNeutralHadronEnergy->Write() ; | |
134ce69a | 1160 | if (fhNeutralHadronPosition) |
1161 | fhNeutralHadronPosition->Write() ; | |
b9bbdad1 | 1162 | if (fhNeutralEMEnergy) |
1163 | fhNeutralEMEnergy->Write() ; | |
134ce69a | 1164 | if (fhNeutralEMPosition) |
1165 | fhNeutralEMPosition->Write() ; | |
b9bbdad1 | 1166 | if (fhChargedHadronEnergy) |
1167 | fhChargedHadronEnergy->Write() ; | |
134ce69a | 1168 | if (fhChargedHadronPosition) |
1169 | fhChargedHadronPosition->Write() ; | |
c1d256cb | 1170 | if (fhPhotonHadronEnergy) |
1171 | fhPhotonHadronEnergy->Write() ; | |
134ce69a | 1172 | if (fhPhotonHadronPosition) |
1173 | fhPhotonHadronPosition->Write() ; | |
92862013 | 1174 | |
1175 | output.Write(); | |
1176 | output.Close(); | |
1177 | } | |
134ce69a | 1178 | |
1179 | ||
1180 |