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