Inheritance from TObject. Automatic streamers.
[u/mrichter/AliRoot.git] / PHOS / testsimglobal.C
1 // Root
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"
10 #include "TString.h"
11 #include "TF1.h"
12 #include "TFormula.h"
13 #include "TFolder.h"
14 #include "TStopwatch.h"
15 #include "TObjArray.h"
16
17 //AliRoot
18 #include "AliPHOSHit.h"
19 #include "AliPHOSGetter.h"
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
31 Bool_t sim_exam() ;
32 Bool_t sdigit() ;
33 Bool_t digit() ;
34 Bool_t recpoint() ;
35 Bool_t track() ;
36 Bool_t particle() ;
37 void write_info(TString) ;
38
39
40 //____________________________________________________________________________ 
41 void testsimglobal() 
42 {
43   
44   gSystem->SetIncludePath("-I$ROOTSYS/include -I$ALICE_ROOT/include -I.");
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 {
107     gSystem->Exec("uuencode $ALICE_ROOT/PHOS/testPHOS.root testPHOS.root | mail -s 'PHOS INSTALLATION ERROR' schutz@in2p3.fr") ;
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 }
112
113
114 //____________________________________________________________________________ 
115 Bool_t sim_exam()
116 {
117   
118   // Definition of the alarm bounds for the examination
119   const Float_t maxAlaHitsM = 12.79  ;  // total multiplicity
120   const Float_t maxAlaTotEn = 19.34  ;    // total energy
121   const Float_t maxAlaTotEnB = 18.35 ;     // per block multiplicity
122   const Float_t maxAlaHitsMB = 11.1  ;  // per block energy
123   
124   TString mess("") ;
125
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   
140   // Examine the alarms
141   TObjArray * alahm = dynamic_cast<TObjArray*>(gime->Alarms()->FindObject("HitsM")) ;
142   Float_t ratiohm = 100.0*static_cast<Float_t>(alahm->GetEntries())/static_cast<Float_t>(maxevent) ;
143   
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) ;
146   
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 ;
164       mess = "Examination detected an error in " + TString(namemul) ;
165       write_info(mess) ;
166     }
167     
168     if (ratioenb[i] > maxAlaTotEnB) {
169       error = kTRUE ;
170       mess = "Examination detected an error in " + TString(namen)  ;
171       write_info(mess) ;
172     }
173   }
174   
175   
176   timer.Stop() ;
177   timer.Print() ;
178   
179   
180   if (ratiohm > maxAlaHitsM){
181     error = kTRUE ;
182     mess = "Examination detected an error in HitsM." ;
183     write_info(mess) ;
184   }
185   
186   if (ratioet>maxAlaTotEn){
187     error = kTRUE ;
188     mess = "Examination detected an error in TotEn." ;
189     write_info(mess) ;
190   }
191   
192   // Condition that will launch the general loop that builds the histograms in order to allow a further analysis.
193   
194   if (!error)
195     return kTRUE ;
196   else {
197     mess = "Examination sets up the file that will be sent to PHOS director (30s)." ;
198     write_info(mess) ;
199     
200     Int_t index = 0 ;
201     Int_t nhits = 0 ;
202     
203     TH1F * his    = new TH1F("Total Multiplicity", "Total Multiplicity in PHOS", 200, 0., 200.) ;
204     TH1F * hisnrg = new TH1F("Total Energy", "Total energy in PHOS", 200, 0., 12.) ;
205     
206     // name and define the different histograms per block involved in the analysis
207     
208     TClonesArray hisba("TH1F") ;
209     TClonesArray hisbanrg("TH1F") ;
210     Int_t i = 0 ;
211     
212     char name[80], title[80] ;
213     for (i = 0 ; i < 6 ; i++) {
214       sprintf(name,"multiplicity for block %d",i) ;
215       sprintf(title,"multiplicity per blocks, block %d",i) ;
216       new(hisba[i]) TH1F(name, title, 100, 0., 100.) ;
217       
218       sprintf(name,"total energy for block %d",i) ;
219       sprintf(title,"total energy per block, block %d",i) ;
220       new(hisbanrg[i]) TH1F(name, title, 200, 0., 12.) ;
221     }
222     
223     // define the global histograms, the hit pointer and give the means to get the actual block reached by the photon
224     
225     AliPHOSHit * hit ;
226     TH1F * hist ;
227     TH1F * histnrg ;
228     TH2F * hbiz = new TH2F ("hbiz","hbiz", 200, 0., 200., 200, 0., 12.) ;
229     const AliPHOSGeometry * geom = gime->PHOSGeometry() ;
230     
231     // the very general loop
232     
233     for (index = 0 ; index < maxevent ; index++) {
234       //get the number of the event
235       gime->Event(index) ;
236       // get the number of cells reached during this event and fill the total multiplicity histogram
237       Int_t n = gime->Hits()->GetEntries() ;
238       nhits += n ;
239       his->Fill(n) ;
240       // Get the data per block
241       const TClonesArray * hita = static_cast<const TClonesArray *>(gime -> Hits()) ;
242       TIter next(hita) ;
243       Float_t Et = 0. ;
244       Int_t id = 0, block = 0 ;
245       Int_t nhit[6], rid[4] ;
246       Float_t etblock[6] ;
247       Int_t i = 0 ;
248       
249       for ( i = 0; i < 6 ; i++) {
250         nhit[i] = 0 ;
251         etblock[i] = 0 ;
252       }
253       
254       while ( (hit = static_cast<AliPHOSHit *>(next())) ) {
255         id = hit->GetId() ;
256         if  (geom->IsInEMC(id) ) {
257           Et += ( hit -> GetEnergy()) ;
258           geom->AbsToRelNumbering(id,rid) ;
259           block = rid[0] ;
260           nhit[block]++ ;
261           etblock[block] +=  ( hit -> GetEnergy()) ;
262         }
263       }
264       
265       //Fill all the histograms but total multiplicity, already done
266       hist = static_cast<TH1F*>(hisba.At(block)) ;
267       hist->Fill(nhit[block]) ;
268       histnrg = static_cast<TH1F*>(hisbanrg.At(block)) ;
269       histnrg->Fill(etblock[block]) ;
270       hisnrg -> Fill(Et) ;
271       hbiz->Fill(n,Et) ;
272     }
273          
274     TFile * file = gROOT -> GetFile("testPHOS.root") ;
275     file -> Write() ;
276     file->Close() ;
277     return kFALSE ;
278   }
279 }
280
281
282 //____________________________________________________________________________ 
283 Bool_t sdigit()
284 {
285   //SDigits process
286   
287   const Float_t maxSDigits = 62.89 ;
288   const Float_t widSDigits = TMath::Sqrt(maxSDigits) ;
289
290   TString mess("") ;
291   TString reconame = "test suite" ;
292   
293   AliPHOSSDigitizer *sd = new AliPHOSSDigitizer("testPHOS.root",reconame.Data()) ;
294   AliPHOSGetter * gime = AliPHOSGetter::GetInstance() ;
295   
296   sd->ExecuteTask("deb") ;
297   
298   Float_t nSDigits =  static_cast<Float_t>(gime->SDigitizer()->GetSDigitsInRun()) / static_cast<Float_t>(gime->MaxEvent()) ;
299   if ( nSDigits < maxSDigits-widSDigits ||
300        nSDigits > maxSDigits+widSDigits ) {
301     mess = "Error detected in the SDigits process. Sending error file to PHOS director." ;
302     cout <<  "sdigit() : nsDigits = " << nSDigits 
303          << " maxSDigits,widSDigits= " << maxSDigits << "," << widSDigits << endl ;    
304     write_info(mess) ;
305     return kFALSE ;
306   }
307   else
308     return kTRUE ;
309 }
310
311
312 //____________________________________________________________________________ 
313 Bool_t digit()
314 {
315   
316   //Digits process
317   AliPHOSGetter * gime = AliPHOSGetter::GetInstance("testPHOS.root") ;
318   TString reconame = "test suite" ;
319   const Float_t maxDigits = 2860. ;
320   const Float_t widDigits = TMath::Sqrt(maxDigits) ;
321   
322   TString mess("") ;
323
324   AliPHOSDigitizer *d = new AliPHOSDigitizer("testPHOS.root",reconame.Data()) ;
325   
326   d->ExecuteTask("deb") ;
327   
328   Float_t nDigits = static_cast<Float_t>(gime->Digitizer()->GetDigitsInRun()) / static_cast<Float_t>(gime->MaxEvent()) ;
329   
330   if ( nDigits < maxDigits-widDigits || nDigits > maxDigits+widDigits ) {
331     cout <<  "digit() : nDigits = " << nDigits 
332          << " maxDigits,widDigits= " << maxDigits << "," << widDigits << endl ;    
333     mess = "Error detected in the Digits process. Sending error file to PHOS director." ;
334     write_info(mess) ;
335     return kFALSE ;
336   }
337   else
338     return kTRUE ;
339 }
340
341
342 //____________________________________________________________________________ 
343 Bool_t recpoint()
344 {
345   
346   //RecPoints process
347   AliPHOSGetter * gime = AliPHOSGetter::GetInstance("testPHOS.root") ;
348   TString reconame = "test suite" ;
349   
350   const Float_t maxRecPoints = 1.0 ;
351   const Float_t widRecPoints = TMath::Sqrt(maxRecPoints) ;
352   
353   TString mess("") ;
354
355   AliPHOSClusterizer * cluster = new  AliPHOSClusterizerv1("testPHOS.root", reconame.Data()) ;
356   
357   cluster->ExecuteTask("deb") ;
358   
359   Float_t nRecPoints =  static_cast<Float_t>(gime->Clusterizer(reconame.Data())->GetRecPointsInRun()) /
360     static_cast<Float_t>(gime->MaxEvent()) ;
361   
362   if ( nRecPoints < maxRecPoints-widRecPoints
363        || nRecPoints > maxRecPoints+widRecPoints ) {
364     cout <<  "recpoint() : nRecPoints = " << nRecPoints 
365          << " maxRecPoints,widRecPoints= " << maxRecPoints << "," << widRecPoints << endl ;    
366     mess = "Error detected in the Clusterizing process. Sending error file to PHOS director." ;
367     write_info(mess) ;
368     return kFALSE ;
369   }
370   else
371     return kTRUE ;
372 }
373
374
375 //____________________________________________________________________________ 
376 Bool_t track()
377 {
378   
379   //TrackSegments process
380   AliPHOSGetter * gime = AliPHOSGetter::GetInstance("testPHOS.root") ;
381   const Float_t maxTrackSegments = 1 ;
382
383   TString mess("") ;
384   TString reconame = "test suite" ;
385   
386   AliPHOSTrackSegmentMaker * tracks = new AliPHOSTrackSegmentMakerv1("testPHOS.root",reconame.Data()) ;
387   
388   tracks->ExecuteTask("deb") ;
389   
390   Float_t nTrackSegments =  static_cast<Float_t> (gime->TrackSegmentMaker(reconame.Data())->GetTrackSegmentsInRun()) /
391     static_cast<Float_t>(gime->MaxEvent()) ;
392   
393   if ( nTrackSegments < maxTrackSegments-0.25 ||
394        nTrackSegments > maxTrackSegments+0.25 ) {
395     cout <<  "track() : nTrackSegments = " << nTrackSegments
396          << " maxTrackSegments,widTrackSegments= " << maxTrackSegments << "," << "0.25" << endl ;    
397     mess = "Error detected in the TrackSegments process. Sending error file to PHOS director." ;
398     write_info(mess) ;
399     return kFALSE ;
400   }
401   else
402     return kTRUE ;
403 }
404
405
406 //____________________________________________________________________________ 
407 Bool_t particle()
408 {
409   
410   //RecParticles process
411   AliPHOSGetter * gime = AliPHOSGetter::GetInstance("testPHOS.root") ;
412   const Float_t maxRecParticles = 1 ;
413
414   TString mess("") ;
415   TString reconame = "test suite" ;
416   
417   AliPHOSPID * pid = new AliPHOSPIDv1("testPHOS.root",reconame.Data()) ;
418   
419   pid->ExecuteTask("deb") ;
420   
421   Float_t nRecParticles =  static_cast<Float_t> (gime->PID(reconame.Data())->GetRecParticlesInRun()) /
422     static_cast<Float_t>(gime->MaxEvent()) ;
423   
424   
425   if ( nRecParticles < maxRecParticles-0.25 ||
426        nRecParticles > maxRecParticles+0.25 ) {
427     cout <<  "particle() : nRecParticles = " << nRecParticles 
428          << " maxRecParticles,widRecParticles= " << maxRecParticles << "," << "0.25" << endl ;    
429     mess = "Error detected in the RecParticles process. Sending error file to PHOS director.Stop reconstruction." ;
430     write_info(mess) ;
431     return kFALSE ;
432   }
433   else
434     return kTRUE ;
435 }
436
437
438 //____________________________________________________________________________ 
439 void write_info(TString mess)
440 {
441   cerr << " _____________________________________________________________ " << endl
442        << " " << endl
443        << "MESS ==> " << mess <<endl
444        << " _____________________________________________________________ " <<endl ;
445 }
446