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       TClonesArray * hita = static_cast<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     write_info(mess) ;
 303     return kFALSE ;
 304   }
 305   else
 306     return kTRUE ;
 307 }
 308
 309
 310 //____________________________________________________________________________
 311 Bool_t digit()
 312 {
 313
 314   //Digits process
 315   AliPHOSGetter * gime = AliPHOSGetter::GetInstance("testPHOS.root") ;
 316   TString reconame = "test suite" ;
 317   const Float_t maxDigits = 2860. ;
 318   const Float_t widDigits = TMath::Sqrt(maxDigits) ;
 319
 320   TString mess("") ;
 321
 322   AliPHOSDigitizer *d = new AliPHOSDigitizer("testPHOS.root",reconame.Data()) ;
 323
 324   d->ExecuteTask("deb") ;
 325
 326   Float_t nDigits = static_cast<Float_t>(gime->Digitizer()->GetDigitsInRun()) / static_cast<Float_t>(gime->MaxEvent()) ;
 327
 328   if ( nDigits < maxDigits-widDigits || nDigits > maxDigits+widDigits ) {
 329     mess = "Error detected in the Digits process. Sending error file to PHOS director." ;
 330     write_info(mess) ;
 331     return kFALSE ;
 332   }
 333   else
 334     return kTRUE ;
 335 }
 336
 337
 338 //____________________________________________________________________________
 339 Bool_t recpoint()
 340 {
 341
 342   //RecPoints process
 343   AliPHOSGetter * gime = AliPHOSGetter::GetInstance("testPHOS.root") ;
 344   TString reconame = "test suite" ;
 345
 346   const Float_t maxRecPoints = 1.0 ;
 347   const Float_t widRecPoints = TMath::Sqrt(maxRecPoints) ;
 348
 349   TString mess("") ;
 350
 351   AliPHOSClusterizer * cluster = new  AliPHOSClusterizerv1("testPHOS.root", reconame.Data()) ;
 352
 353   cluster->ExecuteTask("deb") ;
 354
 355   Float_t nRecPoints =  static_cast<Float_t>(gime->Clusterizer(reconame.Data())->GetRecPointsInRun()) /
 356     static_cast<Float_t>(gime->MaxEvent()) ;
 357
 358   if ( nRecPoints < maxRecPoints-widRecPoints
 359        || nRecPoints > maxRecPoints+widRecPoints ) {
 360     mess = "Error detected in the Clusterizing process. Sending error file to PHOS director." ;
 361     write_info(mess) ;
 362     return kFALSE ;
 363   }
 364   else
 365     return kTRUE ;
 366 }
 367
 368
 369 //____________________________________________________________________________
 370 Bool_t track()
 371 {
 372
 373   //TrackSegments process
 374   AliPHOSGetter * gime = AliPHOSGetter::GetInstance("testPHOS.root") ;
 375   const Float_t maxTrackSegments = 1 ;
 376
 377   TString mess("") ;
 378   TString reconame = "test suite" ;
 379
 380   AliPHOSTrackSegmentMaker * tracks = new AliPHOSTrackSegmentMakerv1("testPHOS.root",reconame.Data()) ;
 381
 382   tracks->ExecuteTask("deb") ;
 383
 384   Float_t nTrackSegments =  static_cast<Float_t> (gime->TrackSegmentMaker(reconame.Data())->GetTrackSegmentsInRun()) /
 385     static_cast<Float_t>(gime->MaxEvent()) ;
 386
 387   if ( nTrackSegments < maxTrackSegments-0.25 ||
 388        nTrackSegments > maxTrackSegments+0.25 ) {
 389     mess = "Error detected in the TrackSegments process. Sending error file to PHOS director." ;
 390     write_info(mess) ;
 391     return kFALSE ;
 392   }
 393   else
 394     return kTRUE ;
 395 }
 396
 397
 398 //____________________________________________________________________________
 399 Bool_t particle()
 400 {
 401
 402   //RecParticles process
 403   AliPHOSGetter * gime = AliPHOSGetter::GetInstance("testPHOS.root") ;
 404   const Float_t maxRecParticles = 1 ;
 405
 406   TString mess("") ;
 407   TString reconame = "test suite" ;
 408
 409   AliPHOSPID * pid = new AliPHOSPIDv1("testPHOS.root",reconame.Data()) ;
 410
 411   pid->ExecuteTask("deb") ;
 412
 413   Float_t nRecParticles =  static_cast<Float_t> (gime->PID(reconame.Data())->GetRecParticlesInRun()) /
 414     static_cast<Float_t>(gime->MaxEvent()) ;
 415
 416
 417   if ( nRecParticles < maxRecParticles-0.25 ||
 418        nRecParticles > maxRecParticles+0.25 ) {
 419     mess = "Error detected in the RecParticles process. Sending error file to PHOS director.Stop reconstruction." ;
 420     write_info(mess) ;
 421     return kFALSE ;
 422   }
 423   else
 424     return kTRUE ;
 425 }
 426
 427
 428 //____________________________________________________________________________
 429 void write_info(TString mess)
 430 {
 431   cerr << " _____________________________________________________________ " << endl
 432        << " " << endl
 433        << "MESS ==> " << mess <<endl
 434        << " _____________________________________________________________ " <<endl ;
 435 }
 436