PHOS/testsimexam.C

   1 #include "AliPHOSGetter.h"
   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 "AliPHOSHit.h"
  11 #include "TF1.h"
  12 #include "TFormula.h"
  13 #include "TFolder.h"
  14 #include "TStopwatch.h"
  15 #include "TObjArray.h"
  16 #include "AliPHOSGeometry.h"
  17
  18
  19
  20
  21 void testsimexam (Int_t nevent=1, const char *config="testconfig.C")
  22 {
  23   //
  24   // Simple macro to run aliroot in a batch mode
  25   //
  26   cerr<<" __________________________________________________________________ "<<endl;
  27   cerr<< " " <<endl;
  28   cerr<<"           MESS ==> Beginning of the simulation examination."<<endl;
  29   cerr<<" __________________________________________________________________ "<<endl;
  30   // Definition of the alarm bounds
  31      const Float_t maxAlaHitsM = 12.79  ;  // total multiplicity
  32      const Float_t maxAlaTotEn = 19.34  ;  // total energy
  33      const Float_t maxAlaTotEnB = 18.35 ;  // per block multiplicity
  34      const Float_t maxAlaHitsMB = 11.1  ;  // per block energy
  35
  36        // boolean which will test if there is an alarm
  37
  38    Int_t boolala;
  39    boolala = 0;
  40    // define the array in which the events that have not reach th EMCA will be put.
  41
  42
  43
  44    TStopwatch timer;
  45    timer.Start();
  46
  47    // Get the number of events generated in the simulation
  48
  49   AliPHOSGetter * gime = AliPHOSGetter::GetInstance("galice.root") ;
  50   Int_t maxevent = gime->MaxEvent() ;
  51
  52   // Examine the alarms
  53
  54    TObjArray * alahm = (TObjArray*)(gime->Alarms()->FindObject("HitsM"));
  55    Float_t ratiohm = 100.0*(Float_t)alahm->GetEntries()/(Float_t)maxevent ;
  56
  57
  58    TObjArray * alaet = (TObjArray*)(gime->Alarms()->FindObject("TotEn"));
  59    Float_t ratioet = 100.0*(Float_t)alaet->GetEntries()/(Float_t)maxevent ;
  60
  61
  62    // Define the alarms per block and examine them
  63    char namemul[80], namen[80];
  64    TObjArray* alahmb[5];
  65    TObjArray* alaenb[5];
  66    Float_t ratiohmb[5], ratioenb[5];
  67    for (Int_t i = 0 ; i < 5 ; i++)
  68      {
  69         sprintf(namemul,"HitsMB%d",i+1);
  70         sprintf(namen,"TotEnB%d",i+1);
  71         alahmb[i] = (TObjArray*) (gime->Alarms()->FindObject(namemul));
  72         ratiohmb[i] = 100.0*(Float_t)alahmb[i]->GetEntries()/(Float_t)maxevent;
  73         alaenb[i] = (TObjArray*)(gime->Alarms()->FindObject(namen));
  74         ratioenb[i] = 100.0*(Float_t)alaenb[i]->GetEntries()/(Float_t)maxevent;
  75         if (ratiohmb[i]>maxAlaHitsMB){
  76           boolala = 1;
  77           cerr<<" _____________________________________________________________ "<<endl;
  78           cerr<< " " <<endl;
  79           cerr << "       MESS ==> Examination detected an error in "<<namemul << endl ;
  80           cerr<<" _____________________________________________________________ "<<endl;
  81         }
  82         if (ratioenb[i]>maxAlaTotEnB) {
  83           boolala = 1;
  84           cerr<<" _____________________________________________________________ "<<endl;
  85           cerr<< " " <<endl;
  86           cerr << "       MESS ==> Examination detected an error in "<<namen << endl ;
  87           cerr<<" _____________________________________________________________ "<<endl;
  88         }
  89
  90      }
  91
  92
  93   timer.Stop();
  94   timer.Print();
  95
  96
  97   if (ratiohm>maxAlaHitsM){
  98      boolala = 1;
  99      cerr<<" _____________________________________________________________ "<<endl;
 100      cerr<< " " <<endl;
 101      cerr << "          MESS ==> Examination detected an error in HitsM." << endl ;
 102      cerr<<" _____________________________________________________________ "<<endl;
 103   }
 104
 105   if (ratioet>maxAlaTotEn){
 106      boolala = 1;
 107      cerr<<" _____________________________________________________________ "<<endl;
 108      cerr<< " " <<endl;
 109      cerr << "         MESS ==> Examination detected an error in TotEn." << endl ;
 110      cerr<<" _____________________________________________________________ "<<endl;
 111   }
 112
 113
 114   // Condition that will launch the general loop that builds the histograms in order to allow a further analysis.
 115
 116     if ( boolala == 1 ) {
 117   cerr<<" _____________________________________________________________ "<<endl;
 118   cerr<< " " <<endl;
 119   cerr << "          MESS ==> Examination sets up the file that will be sent to the responsible for PHOS in ALIROOT development (30s). " << endl ;
 120   cerr<<" _____________________________________________________________ "<<endl;
 121       Int_t index = 0 ;
 122       Int_t nhits = 0 ;
 123
 124       TH1F * his = new TH1F("Total Multiplicity", "Total Multiplicity in PHOS", 200, 0., 200.) ;
 125       TH1F * hisnrg = new TH1F("Total Energy", "Total energy in PHOS", 200, 0., 12.) ;
 126
 127      // name and define the different histograms per block involved in the analysis
 128       TClonesArray hisba("TH1F") ;
 129       TClonesArray hisbanrg("TH1F") ;
 130       for (Int_t i = 0 ; i < 6 ; i++)
 131        {
 132          char name[80], title[80] ;
 133          sprintf(name,"multiplicity for block %d",i) ;
 134          sprintf(title,"multiplicity per blocks, block %d",i) ;
 135          new(hisba[i]) TH1F(name, title, 100, 0., 100.) ;
 136
 137
 138          char namenrg[80], titlenrg[80] ;
 139          sprintf(namenrg,"total energy for block %d",i) ;
 140          sprintf(titlenrg,"total energy per block, block %d",i) ;
 141          new(hisbanrg[i]) TH1F(namenrg, titlenrg, 200, 0., 12.) ;
 142        }
 143       // define the global histograms, the hit pointer and give the means to get the actual block reached by the photon
 144       AliPHOSHit * hit;
 145       TH1F * hist ;
 146       TH1F * histnrg;
 147       TH2F * hbiz = new TH2F ("hbiz","hbiz",200.,0.,200.,200.,0.,12.);
 148       AliPHOSGeometry * geom = AliPHOSGeometry::GetInstance() ;
 149
 150
 151       // the very general loop
 152       for (index = 0 ; index < maxevent ; index++)
 153        {
 154          //get the number of the event
 155          gime->Event(index) ;
 156          // get the number of cells reached during this event and fill the total multiplicity histogram
 157          Int_t n = gime->NHits() ;
 158          nhits += n ;
 159          his->Fill(n) ;
 160          // Get the data per block
 161          TClonesArray * hita = (TClonesArray *) gime -> Hits();
 162          TIter next(hita);
 163          Float_t Et = 0.;
 164          Int_t id = 0, block = 0;
 165          Int_t nhit[6], rid[4];
 166          Float_t etblock[6];
 167          for ( Int_t i = 0; i < 6 ; i++) {
 168            nhit[i] = 0 ;
 169            etblock[i] = 0 ;
 170          }
 171
 172          while ( (hit = (AliPHOSHit *) next()) ) {
 173                  id = hit->GetId();
 174                  if  (geom->IsInEMC(id) ) {
 175                    Et += ( hit -> GetEnergy());
 176                    geom->AbsToRelNumbering(id,rid) ;
 177                    block = rid[0];
 178                    nhit[block]++  ;
 179                    etblock[block] +=  ( hit -> GetEnergy());
 180                  }
 181                 }
 182
 183
 184          //Fill all the histograms but total multiplicity, already done
 185          hist = static_cast<TH1F*>(hisba.At(block)) ;
 186          hist->Fill(nhit[block]) ;
 187          histnrg = static_cast<TH1F*>(hisbanrg.At(block)) ;
 188          histnrg->Fill(etblock[block]);
 189          hisnrg -> Fill(Et);
 190          hbiz->Fill(n,Et);
 191        }
 192
 193       nhits /= maxevent ;
 194       cerr << "av = " << nhits << endl ;
 195       TFile * file = gROOT -> GetFile("galice.root");
 196       file -> Write();
 197       his->Draw() ;
 198       hisnrg->Draw() ;
 199       hbiz->Draw();
 200
 201       //Put the histograms in the root file
 202       for (Int_t i = 0 ; i < 6 ; i++)   {
 203          hist = static_cast<TH1F*>(hisba.At(i)) ;
 204          histnrg = static_cast<TH1F*>(hisbanrg.At(i));
 205
 206
 207          // cout << hist << endl << i << endl ;
 208          hist->Draw();
 209          histnrg->Draw();
 210
 211
 212        }
 213
 214       file->Write() ;
 215       hisba.Delete() ;
 216       hisbanrg.Delete() ;
 217       file->Close();
 218       gSystem->Exec("uuencode $ALICE_ROOT/PHOS/galice.root galice.root | mail -s 'PHOS INSTALLATION ERROR  ' schutz@in2p3.fr");
 219
 220     }
 221   cerr<<" _____________________________________________________________ "<<endl;
 222   cerr<< " " <<endl;
 223   cerr << "           MESS ==> Examination ended successfully. " << endl ;
 224   cerr<<" _____________________________________________________________ "<<endl;
 225 }