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