]> git.uio.no Git - u/mrichter/AliRoot.git/blame - PHOS/testsimglobal.C
Updated PID with pi0 identification
[u/mrichter/AliRoot.git] / PHOS / testsimglobal.C
CommitLineData
493ef90d 1// Root
cede0770 2#include "TH2.h"
3#include "TH1.h"
4#include "TFile.h"
5#include "TTree.h"
6#include "TBranch.h"
7#include "TClonesArray.h"
8#include "TCanvas.h"
9#include "TSystem.h"
493ef90d 10#include "TString.h"
cede0770 11#include "TF1.h"
12#include "TFormula.h"
13#include "TFolder.h"
14#include "TStopwatch.h"
15#include "TObjArray.h"
493ef90d 16
17//AliRoot
18#include "AliPHOSHit.h"
19#include "AliPHOSGetter.h"
cede0770 20#include "AliPHOSGeometry.h"
21#include "AliPHOSDigit.h"
22#include "AliPHOSSDigitizer.h"
23#include "AliPHOSDigitizer.h"
24#include "AliPHOSClusterizer.h"
25#include "AliPHOSClusterizerv1.h"
26#include "AliPHOSTrackSegmentMaker.h"
27#include "AliPHOSTrackSegmentMakerv1.h"
28#include "AliPHOSPID.h"
29#include "AliPHOSPIDv1.h"
30
493ef90d 31Bool_t sim_exam() ;
32Bool_t sdigit() ;
33Bool_t digit() ;
34Bool_t recpoint() ;
35Bool_t track() ;
36Bool_t particle() ;
37void write_info(TString) ;
cede0770 38
493ef90d 39
40//____________________________________________________________________________
41void testsimglobal()
cede0770 42{
cede0770 43
94de8339 44 gSystem->SetIncludePath("-I$ROOTSYS/include -I$ALICE_ROOT/include -I.");
493ef90d 45 Bool_t error = kFALSE ;
46 TString mess("") ;
47
48 if ( sim_exam() ){
49 mess = "Examination ended successfully." ;
50 write_info(mess) ;
51 }
52 else
53 error = kTRUE ;
54
55 if (!error){
56 mess = "Beginning of the PHOS reconstruction." ;
57 write_info(mess) ;
58 if(sdigit()){
59 mess = "SDigits process ended successfully." ;
60 write_info(mess) ;
61 }
62 else
63 error = kTRUE ;
64 }
65
66 if (!error){
67 if (digit()){
68 mess = "Digits process ended successfully." ;
69 write_info(mess) ;
70 }
71 else
72 error = kTRUE ;
73 }
74
75 if (!error){
76 if (recpoint()){
77 mess = "Cluster process ended successfully." ;
78 write_info(mess) ;
79 }
80 else
81 error = kTRUE ;
82 }
83
84 if (!error){
85 if (track()){
86 mess = "TrackSegments process ended successfully." ;
87 write_info(mess) ;
88 }
89 else
90 error = kTRUE ;
91 }
92
93 if (!error){
94 if (particle()){
95 mess = "RecParticles process ended successfully." ;
96 write_info(mess) ;
97 }
98 else
99 error = kTRUE ;
100 }
101
102 if(!error){
103 mess = "reconstruction ended successfully." ;
104 write_info(mess) ;
105 }
106 else {
69c46d16 107 gSystem->Exec("uuencode $ALICE_ROOT/PHOS/testPHOS.root testPHOS.root | mail -s 'PHOS INSTALLATION ERROR' schutz@in2p3.fr") ;
493ef90d 108 mess = "Error message to PHOS director has been sent, please wait for his answer before launching the whole simulation again." ;
109 write_info(mess) ;
110 }
111}
cede0770 112
cede0770 113
493ef90d 114//____________________________________________________________________________
115Bool_t sim_exam()
116{
117
118 // Definition of the alarm bounds for the examination
119 const Float_t maxAlaHitsM = 12.79 ; // total multiplicity
94de8339 120 const Float_t maxAlaTotEn = 19.34 ; // total energy
121 const Float_t maxAlaTotEnB = 18.35 ; // per block multiplicity
493ef90d 122 const Float_t maxAlaHitsMB = 11.1 ; // per block energy
123
124 TString mess("") ;
cede0770 125
493ef90d 126 mess = "Beginning of the simulation examination." ;
127 write_info(mess) ;
128
129 // define the array in which the events that have not reach th EMCA will be put.
130 Bool_t error = kFALSE ;
131
132 TStopwatch timer ;
133 timer.Start() ;
134
135 // Get the number of events generated in the simulation
136 AliPHOSGetter * gime = AliPHOSGetter::GetInstance("testPHOS.root") ;
137 Int_t maxevent = gime->MaxEvent() ;
138 gime->Event(0,"Q") ;
139
cede0770 140 // Examine the alarms
607fc376 141 TObjArray * alahm = dynamic_cast<TObjArray*>(dynamic_cast<TFolder*>(gime->Alarms())->FindObject("HitsM")) ;
493ef90d 142 Float_t ratiohm = 100.0*static_cast<Float_t>(alahm->GetEntries())/static_cast<Float_t>(maxevent) ;
cede0770 143
493ef90d 144 TObjArray * alaet = dynamic_cast<TObjArray*>(gime->Alarms()->FindObject("TotEn")) ;
145 Float_t ratioet = 100.0*static_cast<Float_t>(alaet->GetEntries())/static_cast<Float_t>(maxevent) ;
cede0770 146
493ef90d 147 // Define the alarms per block and examine them
148 char namemul[80], namen[80] ;
149 TObjArray* alahmb[5] ;
150 TObjArray* alaenb[5] ;
151 Float_t ratiohmb[5], ratioenb[5] ;
152
153 Int_t i = 0 ;
154 for (i = 0 ; i < 5 ; i++){
155 sprintf(namemul,"HitsMB%d",i+1) ;
156 sprintf(namen,"TotEnB%d",i+1) ;
157 alahmb[i] = (TObjArray*) (gime->Alarms()->FindObject(namemul)) ;
158 ratiohmb[i] = 100.0*(Float_t)alahmb[i]->GetEntries()/(Float_t)maxevent ;
159 alaenb[i] = (TObjArray*)(gime->Alarms()->FindObject(namen)) ;
160 ratioenb[i] = 100.0*(Float_t)alaenb[i]->GetEntries()/(Float_t)maxevent ;
161
162 if (ratiohmb[i] > maxAlaHitsMB){
163 error = kTRUE ;
69c46d16 164 mess = "Examination detected an error in " + TString(namemul) ;
493ef90d 165 write_info(mess) ;
166 }
167
168 if (ratioenb[i] > maxAlaTotEnB) {
169 error = kTRUE ;
69c46d16 170 mess = "Examination detected an error in " + TString(namen) ;
493ef90d 171 write_info(mess) ;
172 }
cede0770 173 }
493ef90d 174
175
176 timer.Stop() ;
177 timer.Print() ;
178
179
180 if (ratiohm > maxAlaHitsM){
181 error = kTRUE ;
607fc376 182 mess = "Examination detected an error in HitsM: ";
183 mess += ratiohm ;
184 mess += " > " ;
185 mess += maxAlaHitsM ;
493ef90d 186 write_info(mess) ;
187 }
188
cede0770 189 if (ratioet>maxAlaTotEn){
493ef90d 190 error = kTRUE ;
607fc376 191 sprintf(mess.Data(), "Examination detected an error in TotEn: %f > %f", ratioet, maxAlaTotEn) ;
493ef90d 192 write_info(mess) ;
cede0770 193 }
493ef90d 194
cede0770 195 // Condition that will launch the general loop that builds the histograms in order to allow a further analysis.
196
493ef90d 197 if (!error)
198 return kTRUE ;
199 else {
200 mess = "Examination sets up the file that will be sent to PHOS director (30s)." ;
201 write_info(mess) ;
cede0770 202
493ef90d 203 Int_t index = 0 ;
204 Int_t nhits = 0 ;
cede0770 205
493ef90d 206 TH1F * his = new TH1F("Total Multiplicity", "Total Multiplicity in PHOS", 200, 0., 200.) ;
207 TH1F * hisnrg = new TH1F("Total Energy", "Total energy in PHOS", 200, 0., 12.) ;
208
209 // name and define the different histograms per block involved in the analysis
210
211 TClonesArray hisba("TH1F") ;
212 TClonesArray hisbanrg("TH1F") ;
213 Int_t i = 0 ;
214
215 char name[80], title[80] ;
216 for (i = 0 ; i < 6 ; i++) {
217 sprintf(name,"multiplicity for block %d",i) ;
218 sprintf(title,"multiplicity per blocks, block %d",i) ;
219 new(hisba[i]) TH1F(name, title, 100, 0., 100.) ;
cede0770 220
493ef90d 221 sprintf(name,"total energy for block %d",i) ;
222 sprintf(title,"total energy per block, block %d",i) ;
223 new(hisbanrg[i]) TH1F(name, title, 200, 0., 12.) ;
224 }
225
226 // define the global histograms, the hit pointer and give the means to get the actual block reached by the photon
227
228 AliPHOSHit * hit ;
229 TH1F * hist ;
230 TH1F * histnrg ;
231 TH2F * hbiz = new TH2F ("hbiz","hbiz", 200, 0., 200., 200, 0., 12.) ;
232 const AliPHOSGeometry * geom = gime->PHOSGeometry() ;
233
234 // the very general loop
235
236 for (index = 0 ; index < maxevent ; index++) {
237 //get the number of the event
238 gime->Event(index) ;
239 // get the number of cells reached during this event and fill the total multiplicity histogram
240 Int_t n = gime->Hits()->GetEntries() ;
241 nhits += n ;
242 his->Fill(n) ;
243 // Get the data per block
75236460 244 const TClonesArray * hita = static_cast<const TClonesArray *>(gime -> Hits()) ;
493ef90d 245 TIter next(hita) ;
246 Float_t Et = 0. ;
247 Int_t id = 0, block = 0 ;
248 Int_t nhit[6], rid[4] ;
249 Float_t etblock[6] ;
250 Int_t i = 0 ;
cede0770 251
493ef90d 252 for ( i = 0; i < 6 ; i++) {
253 nhit[i] = 0 ;
254 etblock[i] = 0 ;
255 }
256
257 while ( (hit = static_cast<AliPHOSHit *>(next())) ) {
258 id = hit->GetId() ;
259 if (geom->IsInEMC(id) ) {
260 Et += ( hit -> GetEnergy()) ;
261 geom->AbsToRelNumbering(id,rid) ;
262 block = rid[0] ;
263 nhit[block]++ ;
264 etblock[block] += ( hit -> GetEnergy()) ;
265 }
266 }
267
268 //Fill all the histograms but total multiplicity, already done
269 hist = static_cast<TH1F*>(hisba.At(block)) ;
270 hist->Fill(nhit[block]) ;
271 histnrg = static_cast<TH1F*>(hisbanrg.At(block)) ;
272 histnrg->Fill(etblock[block]) ;
273 hisnrg -> Fill(Et) ;
274 hbiz->Fill(n,Et) ;
cede0770 275 }
493ef90d 276
277 TFile * file = gROOT -> GetFile("testPHOS.root") ;
278 file -> Write() ;
279 file->Close() ;
280 return kFALSE ;
281 }
282}
cede0770 283
284
493ef90d 285//____________________________________________________________________________
286Bool_t sdigit()
287{
288 //SDigits process
607fc376 289 AliPHOSGetter * gime = AliPHOSGetter::GetInstance("testPHOS.root") ;
cede0770 290 const Float_t maxSDigits = 62.89 ;
291 const Float_t widSDigits = TMath::Sqrt(maxSDigits) ;
493ef90d 292
293 TString mess("") ;
294 TString reconame = "test suite" ;
295
296 AliPHOSSDigitizer *sd = new AliPHOSSDigitizer("testPHOS.root",reconame.Data()) ;
493ef90d 297
298 sd->ExecuteTask("deb") ;
299
607fc376 300 Float_t nSDigits = static_cast<Float_t>((dynamic_cast<AliPHOSSDigitizer*>(gime->SDigitizer()))->GetSDigitsInRun()) / static_cast<Float_t>(gime->MaxEvent()) ;
493ef90d 301 if ( nSDigits < maxSDigits-widSDigits ||
302 nSDigits > maxSDigits+widSDigits ) {
303 mess = "Error detected in the SDigits process. Sending error file to PHOS director." ;
25b01074 304 cout << "sdigit() : nsDigits = " << nSDigits
305 << " maxSDigits,widSDigits= " << maxSDigits << "," << widSDigits << endl ;
493ef90d 306 write_info(mess) ;
307 return kFALSE ;
308 }
309 else
310 return kTRUE ;
311}
312
313
314//____________________________________________________________________________
315Bool_t digit()
316{
317
318 //Digits process
319 AliPHOSGetter * gime = AliPHOSGetter::GetInstance("testPHOS.root") ;
320 TString reconame = "test suite" ;
607fc376 321 const Float_t maxDigits = 40. ;
cede0770 322 const Float_t widDigits = TMath::Sqrt(maxDigits) ;
493ef90d 323
324 TString mess("") ;
325
326 AliPHOSDigitizer *d = new AliPHOSDigitizer("testPHOS.root",reconame.Data()) ;
327
328 d->ExecuteTask("deb") ;
329
607fc376 330 Float_t nDigits = static_cast<Float_t>((dynamic_cast<AliPHOSDigitizer*>(gime->Digitizer()))->GetDigitsInRun()) / static_cast<Float_t>(gime->MaxEvent()) ;
493ef90d 331
332 if ( nDigits < maxDigits-widDigits || nDigits > maxDigits+widDigits ) {
25b01074 333 cout << "digit() : nDigits = " << nDigits
334 << " maxDigits,widDigits= " << maxDigits << "," << widDigits << endl ;
493ef90d 335 mess = "Error detected in the Digits process. Sending error file to PHOS director." ;
336 write_info(mess) ;
337 return kFALSE ;
cede0770 338 }
493ef90d 339 else
340 return kTRUE ;
341}
cede0770 342
cede0770 343
493ef90d 344//____________________________________________________________________________
345Bool_t recpoint()
346{
347
348 //RecPoints process
349 AliPHOSGetter * gime = AliPHOSGetter::GetInstance("testPHOS.root") ;
350 TString reconame = "test suite" ;
351
94de8339 352 const Float_t maxRecPoints = 1.0 ;
493ef90d 353 const Float_t widRecPoints = TMath::Sqrt(maxRecPoints) ;
354
355 TString mess("") ;
cede0770 356
493ef90d 357 AliPHOSClusterizer * cluster = new AliPHOSClusterizerv1("testPHOS.root", reconame.Data()) ;
358
359 cluster->ExecuteTask("deb") ;
360
361 Float_t nRecPoints = static_cast<Float_t>(gime->Clusterizer(reconame.Data())->GetRecPointsInRun()) /
362 static_cast<Float_t>(gime->MaxEvent()) ;
363
364 if ( nRecPoints < maxRecPoints-widRecPoints
365 || nRecPoints > maxRecPoints+widRecPoints ) {
25b01074 366 cout << "recpoint() : nRecPoints = " << nRecPoints
367 << " maxRecPoints,widRecPoints= " << maxRecPoints << "," << widRecPoints << endl ;
493ef90d 368 mess = "Error detected in the Clusterizing process. Sending error file to PHOS director." ;
369 write_info(mess) ;
370 return kFALSE ;
371 }
372 else
373 return kTRUE ;
374}
cede0770 375
cede0770 376
493ef90d 377//____________________________________________________________________________
378Bool_t track()
379{
380
381 //TrackSegments process
382 AliPHOSGetter * gime = AliPHOSGetter::GetInstance("testPHOS.root") ;
383 const Float_t maxTrackSegments = 1 ;
384
385 TString mess("") ;
386 TString reconame = "test suite" ;
387
388 AliPHOSTrackSegmentMaker * tracks = new AliPHOSTrackSegmentMakerv1("testPHOS.root",reconame.Data()) ;
389
390 tracks->ExecuteTask("deb") ;
cede0770 391
493ef90d 392 Float_t nTrackSegments = static_cast<Float_t> (gime->TrackSegmentMaker(reconame.Data())->GetTrackSegmentsInRun()) /
393 static_cast<Float_t>(gime->MaxEvent()) ;
394
395 if ( nTrackSegments < maxTrackSegments-0.25 ||
396 nTrackSegments > maxTrackSegments+0.25 ) {
25b01074 397 cout << "track() : nTrackSegments = " << nTrackSegments
398 << " maxTrackSegments,widTrackSegments= " << maxTrackSegments << "," << "0.25" << endl ;
493ef90d 399 mess = "Error detected in the TrackSegments process. Sending error file to PHOS director." ;
400 write_info(mess) ;
401 return kFALSE ;
cede0770 402 }
493ef90d 403 else
404 return kTRUE ;
405}
406
407
408//____________________________________________________________________________
409Bool_t particle()
410{
cede0770 411
493ef90d 412 //RecParticles process
413 AliPHOSGetter * gime = AliPHOSGetter::GetInstance("testPHOS.root") ;
414 const Float_t maxRecParticles = 1 ;
cede0770 415
493ef90d 416 TString mess("") ;
417 TString reconame = "test suite" ;
418
419 AliPHOSPID * pid = new AliPHOSPIDv1("testPHOS.root",reconame.Data()) ;
420
421 pid->ExecuteTask("deb") ;
422
423 Float_t nRecParticles = static_cast<Float_t> (gime->PID(reconame.Data())->GetRecParticlesInRun()) /
424 static_cast<Float_t>(gime->MaxEvent()) ;
425
426
427 if ( nRecParticles < maxRecParticles-0.25 ||
428 nRecParticles > maxRecParticles+0.25 ) {
25b01074 429 cout << "particle() : nRecParticles = " << nRecParticles
9e9b9801 430 << " maxRecParticles,widRecParticles= " << maxRecParticles << "," << "0.25" << endl ;
493ef90d 431 mess = "Error detected in the RecParticles process. Sending error file to PHOS director.Stop reconstruction." ;
432 write_info(mess) ;
433 return kFALSE ;
cede0770 434 }
493ef90d 435 else
436 return kTRUE ;
cede0770 437}
493ef90d 438
439
440//____________________________________________________________________________
441void write_info(TString mess)
442{
443 cerr << " _____________________________________________________________ " << endl
444 << " " << endl
445 << "MESS ==> " << mess <<endl
446 << " _____________________________________________________________ " <<endl ;
447}
448