T0/MakeTrendT0.C

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