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
  37
  38 ClassImp(AliT0CalibWalk)
  39
  40 //________________________________________________________________
  41   AliT0CalibWalk::AliT0CalibWalk():   TNamed(),
  42                                       fWalk(0),
  43                                       fAmpLEDRec(0)
  44 {
  45   //
  46 }
  47
  48 //________________________________________________________________
  49 AliT0CalibWalk::AliT0CalibWalk(const char* name):TNamed(),
  50                                       fWalk(0),
  51                                       fAmpLEDRec(0)
  52 {
  53   TString namst = "Calib_";
  54   namst += name;
  55   SetName(namst.Data());
  56   SetTitle(namst.Data());
  57
  58 }
  59
  60 //________________________________________________________________
  61 AliT0CalibWalk::AliT0CalibWalk(const AliT0CalibWalk& calibda) :
  62   TNamed(calibda),
  63   fWalk(0),
  64   fAmpLEDRec(0)
  65
  66 {
  67 // copy constructor
  68   SetName(calibda.GetName());
  69   SetTitle(calibda.GetName());
  70
  71
  72 }
  73
  74 //________________________________________________________________
  75 AliT0CalibWalk &AliT0CalibWalk::operator =(const AliT0CalibWalk& calibda)
  76 {
  77 // assignment operator
  78   SetName(calibda.GetName());
  79   SetTitle(calibda.GetName());
  80
  81   return *this;
  82 }
  83
  84 //________________________________________________________________
  85 AliT0CalibWalk::~AliT0CalibWalk()
  86 {
  87   //
  88 }
  89 //________________________________________________________________
  90 void AliT0CalibWalk::SetWalk(Int_t ipmt)
  91 {
  92   //read QTC walk graph from external file
  93
  94   Int_t mv, ps;
  95   Int_t x[70000], y[70000], index[70000];
  96   Float_t time[10000],amplitude[10000];
  97   string buffer;
  98   Bool_t down=false;
  99
 100   const char * filename = gSystem->ExpandPathName("$ALICE_ROOT/T0/data/CFD-Amp.txt");
 101   ifstream inFile(filename);
 102   if(!inFile) AliError(Form("Cannot open file %s !",filename));
 103
 104   Int_t i=0;
 105   while(getline(inFile,buffer)){
 106     inFile >> ps >> mv;
 107
 108     x[i]=ps; y[i]=mv;
 109     i++;
 110   }
 111   inFile.close();
 112   cout<<" number of data "<<i<<endl;
 113
 114   TMath::Sort(i, y, index,down);
 115   Int_t amp=0, iin=0, isum=0, sum=0;
 116   Int_t ind=0;
 117   for (Int_t ii=0; ii<i; ii++)
 118     {
 119       ind=index[ii];
 120       if(y[ind] == amp)
 121         {
 122           sum +=x[ind];
 123           iin++;
 124         }
 125       else
 126         {
 127           if(iin>0)
 128             time[isum] = Float_t (sum/(iin));
 129           else
 130             time[isum] =Float_t (x[ind]);
 131           amplitude[isum] = Float_t (amp);
 132           amp=y[ind];
 133           iin=0;
 134           isum++;
 135           sum=0;
 136         }
 137     }
 138
 139   inFile.close();
 140
 141   TGraph* gr = new TGraph(isum, amplitude, time);
 142   fWalk.AddAtAndExpand(gr,ipmt);
 143 }
 144
 145 //________________________________________________________________
 146
 147 void AliT0CalibWalk::SetAmpLEDRec(Int_t ipmt)
 148 {
 149   // read LED walk from external file
 150   Float_t mv, ps;
 151   Float_t x[100], y[100];
 152   string buffer;
 153
 154   const char * filename = gSystem->ExpandPathName("$ALICE_ROOT/T0/data/CFD-LED.txt");
 155    ifstream inFile(filename);
 156   if(!inFile) {AliError(Form("Cannot open file %s !",filename));}
 157
 158   inFile >> mv>>ps;
 159   Int_t i=0;
 160
 161   while(getline(inFile,buffer)){
 162     x[i]=mv; y[i]=ps;
 163     inFile >> mv >> ps;
 164     i++;
 165   }
 166   inFile.close();
 167   Float_t y1[100], x1[100];
 168   for (Int_t ir=0; ir<i; ir++){
 169     y1[ir]=y[i-ir]; x1[ir]=x[i-ir];}
 170   TGraph* gr = new TGraph(i,y1,x1);
 171   fAmpLEDRec.AddAtAndExpand(gr,ipmt);
 172
 173 }
 174
 175 //________________________________________________________________
 176
 177 void AliT0CalibWalk::MakeWalkCorrGraph(const char *laserFile)
 178 {
 179   //make walk corerction for preprocessor
 180
 181   TFile *gFile = TFile::Open(laserFile);
 182   //gSystem->Load("libSpectrum");
 183
 184   Int_t npeaks = 20;
 185   Int_t sigma=3.;
 186   Bool_t down = false;
 187
 188   Int_t index[20];
 189   Char_t buf1[10], buf2[10],buf3[20];
 190
 191
 192   for (Int_t i=0; i<24; i++)
 193     {
 194       if(i>11)  {  sprintf(buf1,"T0_A_%i_CFD",i+1-12); }
 195
 196       if(i<12)  { sprintf(buf1,"T0_C_%i_CFD",i+1); }
 197
 198       sprintf(buf2,"CFD_QTC%i",i+1);
 199       sprintf(buf3,"CFD_LED%i",i+1);
 200       cout<<buf1<<" "<<buf2<<" "<<buf3<<endl;
 201       TH2F *qtc_cfd = (TH2F*) gFile->Get(buf2);
 202       TH2F *led_cfd = (TH2F*) gFile->Get(buf3);
 203       TH1F *cfd = (TH1F*) gFile->Get(buf1);
 204       //       cfd->Draw();
 205       TSpectrum *s = new TSpectrum(2*npeaks,1);
 206       Int_t nfound = s->Search(cfd,sigma," ",0.05);
 207       cout<<"Found "<<nfound<<" peaks sigma "<<sigma<<endl;;
 208       if(nfound!=0){
 209         Float_t *xpeak = s->GetPositionX();
 210         TMath::Sort(nfound, xpeak, index,down);
 211         Float_t xp = xpeak[index[0]];
 212         Int_t xbin = cfd->GetXaxis()->FindBin(xp);
 213         Float_t yp = cfd->GetBinContent(xbin);
 214         cout<<"xbin = "<<xbin<<"\txpeak = "<<xpeak[1]<<"\typeak = "<<yp<<endl;
 215         Float_t hmax = xp+10*sigma;
 216         Float_t hmin = xp-10*sigma;
 217         cout<<hmin<< " "<<hmax<<endl;
 218         //          cfd->GetXaxis()->SetRange(hmin,hmax);
 219         //          TF1 *g1 = new TF1("g1", "gaus", hmin, hmax);
 220         // cfd->Fit("g1","R");
 221         //      Int_t hminbin=  qtc_cfd->GetXaxis()->GetFirst();
 222         //      Int_t hmaxbin=  qtc_cfd->GetXaxis()->GetLast();
 223         Int_t nbins= qtc_cfd->GetXaxis()->GetNbins();
 224         cout<<" nbins "<<nbins<<endl;
 225         //      cout<<"  qtc_cfd "<<hminbin<<" "<<hmaxbin<<" "<<nbins<<endl;
 226         //  qtc_cfd->Draw();
 227         TProfile *pr_y = qtc_cfd->ProfileX();
 228         //      Float_t maxHr=pr_y->GetBinCenter(hmaxbin);
 229
 230         pr_y->SetMaximum(hmax);
 231         pr_y->SetMinimum(hmin);
 232         Int_t np=nbins/20;
 233         Double_t *xx = new Double_t[np];
 234         Double_t *yy = new Double_t[np];
 235         Int_t ng=0;
 236         Double_t yg=0;
 237         for (Int_t ip=1; ip<nbins; ip++)
 238           {
 239             if(ip%20 != 0 ) {
 240               if (pr_y->GetBinContent(ip) !=0)
 241                 yg +=pr_y->GetBinContent(ip);
 242               //        cout<<ng<<" "<<pr_y->GetBinContent(ip)<<" "<<yg<<endl;
 243               ng++;}
 244             else {
 245               xx[ip/20] = Float_t (pr_y->GetBinCenter(ip));
 246               yy[ip/20] = yg/ng;
 247               yg=0;
 248               ng=0;
 249               //                      cout<<ip<<" "<<ip/20<<" "<< xx[ip/20]<<" "<< yy[ip/20]<<endl;
 250             }
 251           }
 252         TGraph *gr = new TGraph(np,xx,yy);
 253         gr->SetMinimum(hmin);
 254         gr->SetMaximum(hmax);
 255         fWalk.AddAtAndExpand(gr,i);
 256         //      fWalk.AddAtAndExpand(gr,i+12);
 257
 258         Int_t nbinsled= led_cfd->GetXaxis()->GetNbins();
 259         cout<<" nbins led "<<nbinsled<<endl;
 260         //      cout<<"  qtc_cfd "<<hminbin<<" "<<hmaxbin<<" "<<nbins<<endl;
 261         //  qtc_cfd->Draw();
 262         TProfile *prled_y = led_cfd->ProfileX();
 263         //      Float_t maxHr=pr_y->GetBinCenter(hmaxbin);
 264
 265         prled_y->SetMaximum(hmax);
 266         prled_y->SetMinimum(hmin);
 267         Int_t npled=nbinsled/20;
 268         Double_t *xxled = new Double_t[np];
 269         Double_t *yyled = new Double_t[np];
 270         Int_t ngled=0;
 271         Double_t ygled=0;
 272         for (Int_t ip=1; ip<nbinsled; ip++)
 273           {
 274             if(ip%20 != 0 ) {
 275               if (prled_y->GetBinContent(ip) !=0)
 276                 ygled +=prled_y->GetBinContent(ip);
 277               //        cout<<ng<<" "<<pr_y->GetBinContent(ip)<<" "<<yg<<endl;
 278               ngled++;}
 279             else {
 280               xxled[ip/20] = Float_t (prled_y->GetBinCenter(ip));
 281               yyled[ip/20] = ygled/ngled;
 282               ygled=0;
 283               ngled=0;
 284             }
 285           }
 286         TGraph *grled = new TGraph(npled,xxled,yyled);
 287         grled->SetMinimum(hmin);
 288         grled->SetMaximum(hmax);
 289         fAmpLEDRec.AddAtAndExpand(grled,i);
 290         //      delete [] xx;
 291         //      delete [] yy;
 292
 293       }
 294
 295     }
 296
 297
 298
 299 }
 300