T0/AliT0CalibWalk.cxx

   1 /**************************************************************************
   2  * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
   3  *                                                                        *
   4  * Author: The ALICE Off-line Project.                                    *
   5  * Contributors are mentioned in the code where appropriate.              *
   6  *                                                                        *
   7  * Permission to use, copy, modify and distribute this software and its   *
   8  * documentation strictly for non-commercial purposes is hereby granted   *
   9  * without fee, provided that the above copyright notice appears in all   *
  10  * copies and that both the copyright notice and this permission notice   *
  11  * appear in the supporting documentation. The authors make no claims     *
  12  * about the suitability of this software for any purpose. It is          *
  13  * provided "as is" without express or implied warranty.                  *
  14  **************************************************************************/
  15
  16 /* $Id$ */
  17
  18 ///////////////////////////////////////////////////////////////////////////////
  19 //                                                                           //
  20 // class T0 walk correction      Alla Maevskaya alla@inr.ru 20.11.2007        //
  21 //                                                                           //
  22 ///////////////////////////////////////////////////////////////////////////////
  23
  24 #include "AliT0CalibWalk.h"
  25 #include "AliLog.h"
  26
  27 #include <TObjArray.h>
  28 #include <TGraph.h>
  29 #include <TFile.h>
  30 #include <TH2F.h>
  31 #include <TMath.h>
  32 #include <TSystem.h>
  33 #include <Riostream.h>
  34 #include <TSpectrum.h>
  35 #include <TProfile.h>
  36 #include <TCanvas.h>
  37 #include <TF1.h>
  38
  39 ClassImp(AliT0CalibWalk)
  40
  41 //________________________________________________________________
  42   AliT0CalibWalk::AliT0CalibWalk():   TNamed(),
  43                                       fWalk(0),
  44                                       fAmpLEDRec(0),
  45                                       fQTC(0),
  46                                       fAmpLED(0)
  47 {
  48   //
  49 }
  50
  51 //________________________________________________________________
  52 AliT0CalibWalk::AliT0CalibWalk(const char* name):TNamed(),
  53                                       fWalk(0),
  54                                       fAmpLEDRec(0),                                                                  fQTC(0),
  55                                       fAmpLED(0)
  56
  57 {
  58   TString namst = "Calib_";
  59   namst += name;
  60   SetName(namst.Data());
  61   SetTitle(namst.Data());
  62
  63 }
  64
  65 //________________________________________________________________
  66 AliT0CalibWalk::AliT0CalibWalk(const AliT0CalibWalk& calibda) :
  67   TNamed(calibda),
  68   fWalk(0),
  69   fAmpLEDRec(0),
  70   fQTC(0),
  71   fAmpLED(0)
  72
  73 {
  74 // copy constructor
  75   SetName(calibda.GetName());
  76   SetTitle(calibda.GetName());
  77
  78
  79 }
  80
  81 //________________________________________________________________
  82 AliT0CalibWalk &AliT0CalibWalk::operator =(const AliT0CalibWalk& calibda)
  83 {
  84 // assignment operator
  85   SetName(calibda.GetName());
  86   SetTitle(calibda.GetName());
  87
  88   return *this;
  89 }
  90
  91 //________________________________________________________________
  92 AliT0CalibWalk::~AliT0CalibWalk()
  93 {
  94   //
  95 }
  96 //________________________________________________________________
  97
  98 Bool_t AliT0CalibWalk::MakeWalkCorrGraph(const char *laserFile)
  99 {
 100   //make walk corerction for preprocessor
 101   Int_t npeaks = 20;
 102   Int_t sigma=3;
 103   Bool_t down=false;
 104   Int_t index[20];
 105   Bool_t ok=true;
 106   Float_t   mips[100];
 107   Int_t nfound=0;
 108   /*
 109   //061110
 110  Float_t timefitth[24] = { 5594, 5618, 5630, 5662,  5623, 5645,
 111                             5645, 5602, 5691, 5669,  5636, 5602,
 112                             5651, 5654, 5612, 5625,  5669, 5678,
 113                             5652, 5673, 5565, 5717,  5692, 5738} ;
 114   */
 115
 116   //051210
 117   Float_t timefitth[24] = { 5587, 5609, 5621, 5653, 5615, 5635,
 118                             5637, 5595, 5683, 5659, 5629, 5601,
 119                             5643, 5654, 5601, 5615, 5661, 5667,
 120                             5645, 5665, 5663, 5709, 5685, 5725 };
 121
 122   Float_t xdata[4000], ydata[4000];
 123   Int_t npoints =0;
 124   TH1F* www[24];
 125   TFile *ff = new TFile("calibwalk0.root");
 126
 127   Int_t ipmt, jpoint;
 128
 129   gFile = TFile::Open(laserFile);
 130   if(!gFile) {
 131     AliError("No input laser data found ");
 132   }
 133   else
 134     {
 135       //      gFile->ls();
 136       TH1F* hAmp = (TH1F*) gFile->Get("hAmpLaser");
 137       Int_t nmips=0;
 138       for (Int_t ibin=0; ibin<1200; ibin++) {
 139         Float_t bincont = hAmp->GetBinContent(ibin);
 140         if(bincont>0){
 141           mips[nmips] = hAmp->GetXaxis()->GetBinCenter(ibin);
 142           cout<<ibin<<" bincont "<<bincont<<" amp "<< mips[nmips]<<endl;
 143           nmips++;
 144         }
 145       }
 146
 147       /*      if (nmips<17) {
 148         ok=false;
 149         return ok;
 150         } */
 151       //  TCanvas *c1 = new TCanvas("c1", "LED-CFD C side",0,48,1280,951);
 152       // c1->Divide(4,3);
 153    //   TCanvas *c2 = new TCanvas("c1", "LED-CFD C side",0,48,1280,951);
 154    //   c2->Divide(4,3);
 155
 156       Float_t x1[50], y1[50];
 157       Float_t x2[50], xx2[50],y2[50];
 158       Float_t xx1[50],yy1[50], xx[50];
 159
 160       Float_t cfd0[24];
 161       Int_t startim = 0;
 162
 163
 164       for (Int_t i=0; i<24; i++)
 165         {
 166           //     if(i==13) continue;
 167           //      else c2->cd(i+1-12);
 168           for (Int_t im=startim; im<nmips; im++)
 169             {
 170
 171               x1[im] = y1[im] = x2[im] = y2[im] = 0;
 172               if(i==23 && im<2) continue;
 173               xx[im]=mips[im];
 174               TString cfd = Form("hCFD%i_%i",i+1,im+1);
 175               TString qtc = Form("hQTC%i_%i",i+1,im+1);
 176               TString led = Form("hLED%i_%i",i+1,im+1);
 177
 178               TH1F *hCFD = (TH1F*) gFile->Get(cfd.Data()) ;
 179               TH1F *hLED = (TH1F*) gFile->Get(led.Data());
 180               TH1F *hQTC = (TH1F*) gFile->Get(qtc.Data()) ;
 181               if(!hCFD )
 182                 AliWarning(Form(" no CFD data in LASER DA for channel %i for amplitude %f MIPs",i,mips[im]));
 183               if(!hQTC )
 184                 AliWarning(Form(" no QTC correction data in LASER DA for channel %i for amplitude %f MIPs",i,mips[im]));
 185               if(!hLED)
 186                 AliWarning(Form(" no LED correction data in LASER DA for channel %i for amplitude %f MIPs",i,mips[im]));
 187
 188               if(hCFD ) {
 189                 hCFD->SetDirectory(0);
 190                 hCFD->GetXaxis()->SetRangeUser(timefitth[i]-100,timefitth[i]+100);
 191
 192                 TF1 *fg = new TF1("fg","gaus",timefitth[i]-50,timefitth[i]+50);
 193                 hCFD->Fit("fg","RQ", " ",timefitth[i]-50,timefitth[i]+50);
 194                 Double_t par[3];
 195                 fg->GetParameters(&par[0]);
 196
 197                 TSpectrum *s = new TSpectrum(2*npeaks,1);
 198                 nfound = s->Search(hCFD,sigma," ",0.1);
 199                 if(nfound!=0){
 200                   Float_t *xpeak = s->GetPositionX();
 201                   TMath::Sort(nfound, xpeak, index,down);
 202                   Float_t xp = xpeak[index[0]];
 203                   Double_t hmax = xp+3*sigma;
 204                   Double_t hmin = xp-3*sigma;
 205                   hCFD->GetXaxis()->SetRangeUser(hmin-10,hmax+10);
 206                   //    if (i==23)
 207                   //              cout<<"PMT "<<i<<" MIPS "<<mips[im]<<" mean "<<hCFD->GetMean()<<
 208                   //     " xp "<<xp<<" fit "<<par[1]<<endl;
 209                   //              cout<<"PMT "<<i<<" MIPS "<<im<<" xp "<<xp<<endl;
 210                 }
 211                 else
 212                   {
 213                     hCFD->Rebin(2);
 214                     TSpectrum *s1 = new TSpectrum(2*npeaks,1);
 215                     nfound = s1->Search(hCFD,sigma," ",0.1);
 216                     if(nfound!=0){
 217                       Float_t *xpeak = s1->GetPositionX();
 218                       TMath::Sort(nfound, xpeak, index,down);
 219                       Float_t xp = xpeak[index[0]];
 220                       Double_t hmax = xp+3*sigma;
 221                       Double_t hmin = xp-3*sigma;
 222                       hCFD->GetXaxis()->SetRangeUser(hmin-10,hmax+10);
 223                     }
 224                     else
 225                       {
 226                         ok=true;
 227                         //                      ok=false;
 228                         printf("no peak in CFD spectrum for PMT %i amplitude %i\n",i,im);
 229                         return ok;
 230                       }
 231
 232                   }
 233                 // } //if no C12 80MIPs
 234                 //        if (i<12) cout<<"!!!! PMT "<<i<<" MIPS "<<im<<" fit  "<<par[1]<<
 235                 //" mean "<<hCFD->GetMean()<<endl;
 236                  if (im == 0) cfd0[i] = par[1];
 237                 if(i==23 && im==2) cfd0[i] = par[1];
 238                    y1[im] =  par[1] - cfd0[i];
 239                    if ((i == 5 && im == 15) ||(i == 8 && im == 19) || (i==11 && im==0) || i==23)
 240                        y1[im] =  hCFD->GetMean() - cfd0[i];
 241
 242                 //              if (im == 0) cfd0[i] = hCFD->GetMean();
 243                 //              y1[im] =  hCFD->GetMean() - cfd0[i];
 244
 245                    //   if (i==23 && im<2) y1[im]=0;
 246                 /*              if( i == 11 && im==10 ){        y1[im] = y1[im-1];
 247                 cout<<im<<" "<<mips[im]<<" mip "<<" pmt "<<i<<" mean "<<hCFD->GetMean()<<" walk "<<     y1[im]<<endl;
 248                 }*/
 249
 250               }
 251               if( hQTC) {
 252                 x1[im] = hQTC->GetMean();
 253                 if( x1[im] == 0) {
 254                   ok=false;
 255                   printf("no peak in QTC signal for PMT %i amplitude %i\n",i,im);
 256                   return ok;
 257                   }
 258               }
 259
 260               if( hLED){
 261                 hLED->SetDirectory(0);
 262                 TSpectrum *s = new TSpectrum(2*npeaks,1);
 263                 nfound = s->Search(hLED,sigma," ",0.1);
 264                 if(nfound!=0){
 265                   Float_t *xpeak = s->GetPositionX();
 266                   TMath::Sort(nfound, xpeak, index,down);
 267                   Float_t xp = xpeak[index[0]];
 268                   Double_t hmax = xp+10*sigma;
 269                   Double_t hmin = xp-10*sigma;
 270                   hLED->GetXaxis()->SetRangeUser(hmin-10,hmax+10);
 271                 }
 272                 else
 273                   {
 274                     ok=false;
 275                     printf("no peak in LED spectrum for PMT %i amplitude %i\n",i,im);
 276                     return ok;
 277                   }
 278                 x2[im] = hLED->GetMean();
 279                 xx2[im] = x2[nmips-im-1];
 280               }
 281
 282               if (hQTC) delete  hQTC;
 283               if (hCFD) delete  hCFD;
 284               if (hLED) delete  hLED;
 285             }
 286
 287           for (Int_t imi=0; imi<nmips; imi++)
 288             {
 289               yy1[imi] = Float_t (mips[nmips-imi-1]);
 290               xx1[imi] = x2[nmips-imi-1];
 291               //              cout<<" LED rev "<<i<<" "<<imi<<" "<<xx1[imi]<<" "<< yy1[imi]<<"nmips-imi-1 = "<< nmips-imi-1<<endl;
 292             }
 293
 294           if(i==0) cout<<"Making graphs..."<<endl;
 295           // TGraph *grwalkqtc = new TGraph (nmips-4,x1,y1);
 296           // grwalkqtc->SetTitle(Form("PMT%i",i));
 297           TGraph *grwalkled = new TGraph (nmips-4,x2,y1);
 298           grwalkled->SetTitle(Form("PMT%i",i));
 299           // fWalk.AddAtAndExpand(grwalkqtc,i);
 300           fAmpLEDRec.AddAtAndExpand(grwalkled,i);
 301           //      cout<<" add walk "<<i<<endl;
 302
 303           //fit amplitude graphs to make comparison wth new one
 304           TGraph *grampled = new TGraph (nmips,xx1,yy1);
 305           TGraph *grqtc = new TGraph (nmips,x1,xx);
 306           fQTC.AddAtAndExpand(grqtc,i);
 307           fAmpLED.AddAtAndExpand(grampled,i);
 308           //      cout<<" add amp "<<i<<endl;
 309
 310           if(i==23)
 311             cout<<"Graphs created..."<<endl;
 312           //       if(i==13) continue;
 313
 314           www[i] = (TH1F*)ff->Get(Form("hy%i",i) );
 315           www[i]->Print();
 316           npoints =0;
 317           for(Int_t ibin=1; ibin< (www[i]->GetNbinsX() -2); ibin++) {
 318              if( www[i]->GetBinContent(ibin)>0 &&
 319                 www[i]->GetBinContent(ibin) !=1000 ) {
 320               ydata[ibin-1] = www[i]->GetBinContent(ibin) - 1000;
 321               xdata[ibin-1] = ibin+799;
 322               cout<<ibin<<" "<<xdata[ibin]<<" "<<ydata[ibin]<<endl;
 323               npoints++;
 324              }
 325           }
 326
 327
 328           TGraph *grwalkqtc = new TGraph (npoints-3,xdata,ydata);
 329           grwalkqtc->SetTitle(Form("PMT%i",i));
 330           fWalk.AddAtAndExpand(grwalkqtc,i);
 331         }
 332
 333
 334     } //if gFile exits
 335
 336   return ok;
 337 }
 338 void AliT0CalibWalk::GetMeanAndSigma(TH1F* hist, Float_t &mean, Float_t &sigma)
 339 {
 340
 341   const double window = 2.;  //fit window
 342
 343   double meanEstimate, sigmaEstimate;
 344   int maxBin;
 345   maxBin        =  hist->GetMaximumBin(); //position of maximum
 346   meanEstimate  =  hist->GetBinCenter( maxBin); // mean of gaussian sitting in maximum
 347   sigmaEstimate = hist->GetRMS();
 348   TF1* fit= new TF1("fit","gaus", meanEstimate - window*sigmaEstimate, meanEstimate + window*sigmaEstimate);
 349   fit->SetParameters(hist->GetBinContent(maxBin), meanEstimate, sigmaEstimate);
 350   hist->Fit("fit","RQ","Q");
 351
 352   mean  = (Float_t) fit->GetParameter(1);
 353   sigma = (Float_t) fit->GetParameter(2);
 354   printf(" mean %f max %f sigma %f \n",mean, meanEstimate , sigma);
 355
 356   delete fit;
 357 }
 358
 359
 360