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