T0/MakeTrendT0.C

   1 #include "TChain.h"
   2 #include "TTree.h"
   3 #include "TH1F.h"
   4 #include "TF1.h"
   5 #include "TH2F.h"
   6 #include "TCanvas.h"
   7 #include "TObjArray.h"
   8 #include "/home/alla/alice/AliRoot/AliRoot/STEER/STEERBase/TTreeStream.h"
   9
  10 #define NPMTs 24
  11
  12 int MakeTrendT0( char *infile, int run) {
  13
  14   char *outfile = "trending.root";
  15
  16   if(!infile) return -1;
  17   if(!outfile) return -1;
  18
  19   TFile *f = TFile::Open(infile,"read");
  20   if (!f) {
  21     printf("File %s not available\n", infile);
  22     return -1;
  23   }
  24
  25   //            LOAD HISTOGRAMS FROM QAresults.root
  26   TObjArray   *fTzeroObject = (TObjArray*) f->Get("T0_Performance/QAT0chists");
  27
  28   TH1F* fTzeroORAplusORC =(TH1F*) ((TH1F*) fTzeroObject->FindObject("fTzeroORAplusORC"))->Clone("A");
  29   TH1F* fResolution      =(TH1F*) ((TH1F*) fTzeroObject->FindObject("fResolution"))->Clone("B");
  30   TH1F* fTzeroORA        =(TH1F*) ((TH1F*) fTzeroObject->FindObject("fTzeroORA"))->Clone("C");
  31   TH1F* fTzeroORC        =(TH1F*) ((TH1F*) fTzeroObject->FindObject("fTzeroORC"))->Clone("D");
  32
  33   TH2F *fTimeVSAmplitude[NPMTs];//counting PMTs from 0
  34   TH1D *fAmplitude[NPMTs];
  35   TH1D *fTime[NPMTs];
  36   //-----> add new histogram here
  37
  38
  39
  40   // sigma fit of resolution, mean T0A, T0C and T0C&A,  mean amplitude for each PMT
  41   double resolutionSigma = -9999; //dummy init
  42   double tzeroOrAPlusOrC = -9999; //dummy init
  43   double tzeroOrA        = -9999; //dummy init
  44   double tzeroOrC        = -9999; //dummy init
  45   double meanAmplitude[NPMTs]; //dummy init
  46   double meanTime[NPMTs]; //dummy init
  47   double timeDelayOCDB[NPMTs]; //dummy init
  48
  49   //.................................................................................
  50
  51   for(int ipmt=1; ipmt<=NPMTs; ipmt++){ //loop over all PMTs
  52     fTimeVSAmplitude[ipmt-1] =(TH2F*) ((TH2F*) fTzeroObject->FindObject(Form("fTimeVSAmplitude%d",ipmt)))->Clone(Form("E%d",ipmt));
  53     int nbinsX = fTimeVSAmplitude[ipmt-1]->GetNbinsX();
  54     int nbinsY = fTimeVSAmplitude[ipmt-1]->GetNbinsY();
  55     //Amplitude
  56     fAmplitude[ipmt-1] = (TH1D*) fTimeVSAmplitude[ipmt-1]->ProjectionX(Form("fAmplitude%d",ipmt), 1, nbinsY);
  57     meanAmplitude[ipmt-1] = -9999; //dummy init
  58     if(fAmplitude[ipmt-1]->GetEntries()>0){
  59       meanAmplitude[ipmt-1] = fAmplitude[ipmt-1]->GetMean();
  60     }
  61
  62     //Time
  63     fTime[ipmt-1] = (TH1D*) fTimeVSAmplitude[ipmt-1]->ProjectionY(Form("fTime%d",ipmt), 1, nbinsX);
  64     meanTime[ipmt-1] = -9999; //dummy init
  65     if(fTime[ipmt-1]->GetEntries()>0){
  66       if(fTime[ipmt-1]->GetEntries()>20){
  67         meanTime[ipmt-1] = GetParameterGaus((TH1F*) fTime[ipmt-1], 1); // Mean Time
  68       }else if(fTime[ipmt-1]->GetEntries()>0){
  69         meanTime[ipmt-1] = fTime[ipmt-1]->GetMean();
  70       }
  71     }
  72   }
  73
  74
  75
  76   if(fResolution->GetEntries()>20){
  77     resolutionSigma = GetParameterGaus(fResolution, 2); //gaussian sigma
  78   }else if(fResolution->GetEntries()>0){
  79     resolutionSigma = fResolution->GetRMS(); //gaussian sigma
  80   }
  81
  82   if(fTzeroORAplusORC->GetEntries()>20){
  83     tzeroOrAPlusOrC = GetParameterGaus(fTzeroORAplusORC, 1); //gaussian mean
  84   }else if(fTzeroORAplusORC->GetEntries()>0){
  85     tzeroOrAPlusOrC = fTzeroORAplusORC->GetMean();
  86   }
  87
  88   if(fTzeroORA->GetEntries()>20){
  89     tzeroOrA = GetParameterGaus(fTzeroORA, 1); //gaussian mean
  90   }else if(fTzeroORA->GetEntries()>0){
  91     tzeroOrA = fTzeroORA->GetMean();
  92   }
  93
  94   if(fTzeroORC->GetEntries()>20){
  95     tzeroOrC = GetParameterGaus(fTzeroORC, 1); //gaussian mean
  96   }else if(fTzeroORC->GetEntries()>0){
  97     tzeroOrC = fTzeroORC->GetMean(); //gaussian mean
  98   }
  99
 100   //-----> analyze the new histogram here and set mean/sigma
 101
 102
 103   f->Close();
 104
 105   //-------------------- READ OCDB TIME DELAYS ---------------------------
 106   // Arguments:
 107   AliCDBManager* man = AliCDBManager::Instance();
 108   man->SetDefaultStorage("raw://");
 109   man->SetRun(run);
 110   AliCDBEntry *entry = AliCDBManager::Instance()->Get("T0/Calib/TimeDelay");
 111   AliT0CalibTimeEq *clb = (AliT0CalibTimeEq*)entry->GetObject();
 112   for (Int_t i=0; i<NPMTs; i++){
 113     timeDelayOCDB[i] = 0;
 114     if(clb)
 115       timeDelayOCDB[i] = clb->GetCFDvalue(i,0);
 116   }
 117
 118
 119
 120   //--------------- write walues to the output ------------
 121   TTreeSRedirector* pcstream = NULL;
 122   pcstream = new TTreeSRedirector(outfile);
 123   if (!pcstream) return;
 124
 125
 126   TFile *x =  pcstream->GetFile();
 127   x->cd();
 128
 129   TObjString runType;
 130   Int_t startTimeGRP=0;
 131   Int_t stopTimeGRP=0;
 132   Int_t time=0;
 133   Int_t duration=0;
 134
 135   time     = (startTimeGRP+stopTimeGRP)/2;
 136   duration = (stopTimeGRP-startTimeGRP);
 137
 138   (*pcstream)<<"t0QA"<<
 139                   "run="<<run;//<<
 140                  //"time="<<time<<
 141             // "startTimeGRP="<<startTimeGRP<<
 142             //  "stopTimeGRP="<<stopTimeGRP<<
 143               ///    "duration="<<duration<<
 144                //   "runType.="<<&runType;
 145
 146
 147  (*pcstream)<<"t0QA"<<
 148              "resolution="<< resolutionSigma<<
 149              "tzeroOrAPlusOrC="<< tzeroOrAPlusOrC<<
 150              "tzeroOrA="<< tzeroOrA<<
 151              "tzeroOrC="<< tzeroOrC<<
 152              "amplPMT1="<<meanAmplitude[0]<<
 153              "amplPMT2="<<meanAmplitude[1]<<
 154              "amplPMT3="<<meanAmplitude[2]<<
 155              "amplPMT4="<<meanAmplitude[3]<<
 156              "amplPMT5="<<meanAmplitude[4]<<
 157              "amplPMT6="<<meanAmplitude[5]<<
 158              "amplPMT7="<<meanAmplitude[6]<<
 159              "amplPMT8="<<meanAmplitude[7]<<
 160              "amplPMT9="<<meanAmplitude[8]<<
 161              "amplPMT10="<<meanAmplitude[9]<<
 162              "amplPMT11="<<meanAmplitude[10]<<
 163              "amplPMT12="<<meanAmplitude[11]<<
 164              "amplPMT13="<<meanAmplitude[12]<<
 165              "amplPMT14="<<meanAmplitude[13]<<
 166              "amplPMT15="<<meanAmplitude[14]<<
 167              "amplPMT16="<<meanAmplitude[15]<<
 168              "amplPMT17="<<meanAmplitude[16]<<
 169              "amplPMT18="<<meanAmplitude[17]<<
 170              "amplPMT19="<<meanAmplitude[18]<<
 171              "amplPMT20="<<meanAmplitude[19]<<
 172              "amplPMT21="<<meanAmplitude[20]<<
 173              "amplPMT22="<<meanAmplitude[21]<<
 174              "amplPMT23="<<meanAmplitude[22]<<
 175              "amplPMT24="<<meanAmplitude[23]<<
 176              "timePMT1="<<meanTime[0]<<
 177              "timePMT2="<<meanTime[1]<<
 178              "timePMT3="<<meanTime[2]<<
 179              "timePMT4="<<meanTime[3]<<
 180              "timePMT5="<<meanTime[4]<<
 181              "timePMT6="<<meanTime[5]<<
 182              "timePMT7="<<meanTime[6]<<
 183              "timePMT8="<<meanTime[7]<<
 184              "timePMT9="<<meanTime[8]<<
 185              "timePMT10="<<meanTime[9]<<
 186              "timePMT11="<<meanTime[10]<<
 187              "timePMT12="<<meanTime[11]<<
 188              "timePMT13="<<meanTime[12]<<
 189              "timePMT14="<<meanTime[13]<<
 190              "timePMT15="<<meanTime[14]<<
 191              "timePMT16="<<meanTime[15]<<
 192              "timePMT17="<<meanTime[16]<<
 193              "timePMT18="<<meanTime[17]<<
 194              "timePMT19="<<meanTime[18]<<
 195              "timePMT20="<<meanTime[19]<<
 196              "timePMT21="<<meanTime[20]<<
 197              "timePMT22="<<meanTime[21]<<
 198              "timePMT23="<<meanTime[22]<<
 199              "timePMT24="<<meanTime[23];
 200 (*pcstream)<<"t0QA"<<
 201              "timeDelayPMT1="<<timeDelayOCDB[0]<<
 202              "timeDelayPMT2="<<timeDelayOCDB[1]<<
 203              "timeDelayPMT3="<<timeDelayOCDB[2]<<
 204              "timeDelayPMT4="<<timeDelayOCDB[3]<<
 205              "timeDelayPMT5="<<timeDelayOCDB[4]<<
 206              "timeDelayPMT6="<<timeDelayOCDB[5]<<
 207              "timeDelayPMT7="<<timeDelayOCDB[6]<<
 208              "timeDelayPMT8="<<timeDelayOCDB[7]<<
 209              "timeDelayPMT9="<<timeDelayOCDB[8]<<
 210              "timeDelayPMT10="<<timeDelayOCDB[9]<<
 211              "timeDelayPMT11="<<timeDelayOCDB[10]<<
 212              "timeDelayPMT12="<<timeDelayOCDB[11]<<
 213              "timeDelayPMT13="<<timeDelayOCDB[12]<<
 214              "timeDelayPMT14="<<timeDelayOCDB[13]<<
 215              "timeDelayPMT15="<<timeDelayOCDB[14]<<
 216              "timeDelayPMT16="<<timeDelayOCDB[15]<<
 217              "timeDelayPMT17="<<timeDelayOCDB[16]<<
 218              "timeDelayPMT18="<<timeDelayOCDB[17]<<
 219              "timeDelayPMT19="<<timeDelayOCDB[18]<<
 220              "timeDelayPMT20="<<timeDelayOCDB[19]<<
 221              "timeDelayPMT21="<<timeDelayOCDB[20]<<
 222              "timeDelayPMT22="<<timeDelayOCDB[21]<<
 223              "timeDelayPMT23="<<timeDelayOCDB[22]<<
 224              "timeDelayPMT24="<<timeDelayOCDB[23];
 225
 226             //-----> add the mean/sigma of the new histogram here
 227
 228  (*pcstream)<<"t0QA"<<"\n";
 229
 230   pcstream->Close();
 231
 232   delete pcstream;
 233
 234   return;
 235
 236 }
 237
 238 //_____________________________________________________________________________
 239 double GetParameterGaus(TH1F *histo, int whichParameter){
 240
 241    int    maxBin   =  histo->GetMaximumBin();
 242    double max      =  (histo->GetBinContent(maxBin-1) + histo->GetBinContent(maxBin) + histo->GetBinContent(maxBin+1))/3;
 243    double mean     =  histo->GetBinCenter(maxBin); //mean
 244    double lowfwhm  =  histo->GetBinCenter(histo->FindFirstBinAbove(max/2));
 245    double highfwhm =  histo->GetBinCenter(histo->FindLastBinAbove(max/2));
 246    double sigma    = (highfwhm - lowfwhm)/2.35482; //estimate fwhm  FWHM = 2.35482*sigma
 247
 248    TF1 *gaussfit   = new TF1("gaussfit","gaus", mean - 4*sigma, mean + 4*sigma); // fit in +- 4 sigma window
 249    gaussfit->SetParameters(max, mean, sigma);
 250
 251    if(whichParameter==2)
 252      histo->Fit(gaussfit,"RQNI");
 253    else
 254      histo->Fit(gaussfit,"RQN");
 255
 256    double parValue = gaussfit->GetParameter(whichParameter);
 257
 258    delete gaussfit;
 259
 260    return parValue;
 261 }
 262
 263