]> git.uio.no Git - u/mrichter/AliRoot.git/blobdiff - T0/AliT0CalibWalk.cxx
Fix in the file open option
[u/mrichter/AliRoot.git] / T0 / AliT0CalibWalk.cxx
index 130b3ca35f2a5d3271049db7d8b7871086240276..43f2910bb88393ee96cad7b7a31a1d6626512023 100644 (file)
@@ -12,7 +12,7 @@
  * about the suitability of this software for any purpose. It is          *
  * provided "as is" without express or implied warranty.                  *
  **************************************************************************/
-
 /* $Id$ */
 
 ///////////////////////////////////////////////////////////////////////////////
 #include <Riostream.h>
 #include <TSpectrum.h>
 #include <TProfile.h>
-
+#include <TF1.h>
 
 ClassImp(AliT0CalibWalk)
 
 //________________________________________________________________
   AliT0CalibWalk::AliT0CalibWalk():   TNamed(),
                                      fWalk(0),
-                                     fAmpLEDRec(0)
+                                     fAmpLEDRec(0),
+                                     fQTC(0),
+                                     fAmpLED(0)
 {
   //
 }
@@ -48,7 +50,9 @@ ClassImp(AliT0CalibWalk)
 //________________________________________________________________
 AliT0CalibWalk::AliT0CalibWalk(const char* name):TNamed(),
                                      fWalk(0),
-                                     fAmpLEDRec(0)                                   
+                                     fAmpLEDRec(0),                                                                  fQTC(0),
+                                     fAmpLED(0)
+    
 {
   TString namst = "Calib_";
   namst += name;
@@ -61,8 +65,10 @@ AliT0CalibWalk::AliT0CalibWalk(const char* name):TNamed(),
 AliT0CalibWalk::AliT0CalibWalk(const AliT0CalibWalk& calibda) :
   TNamed(calibda),             
   fWalk(0),
-  fAmpLEDRec(0)
-
+  fAmpLEDRec(0),
+  fQTC(0),
+  fAmpLED(0)
 {
 // copy constructor
   SetName(calibda.GetName());
@@ -109,7 +115,7 @@ void AliT0CalibWalk::SetWalk(Int_t ipmt)
     i++;
   }
   inFile.close();
-  cout<<" number of data "<<i<<endl;
+  //  cout<<" number of data "<<i<<endl;
  
   TMath::Sort(i, y, index,down);
   Int_t amp=0, iin=0, isum=0, sum=0;
@@ -140,6 +146,12 @@ void AliT0CalibWalk::SetWalk(Int_t ipmt)
 
   TGraph* gr = new TGraph(isum, amplitude, time);
   fWalk.AddAtAndExpand(gr,ipmt);
+
+  //should be change to real
+  Double_t xq[10] = { 1220, 1370, 1542, 1697, 1860, 2023,2171,2331,2495,2684};
+  Double_t yq[10] = {1,2,3,4,5,6,7,8,9,10};
+  TGraph* gr1 = new TGraph(10, xq, yq);
+  fQTC.AddAtAndExpand(gr1,ipmt);
 }
 
 //________________________________________________________________
@@ -169,6 +181,11 @@ void AliT0CalibWalk::SetAmpLEDRec(Int_t ipmt)
     y1[ir]=y[i-ir]; x1[ir]=x[i-ir];}
   TGraph* gr = new TGraph(i,y1,x1);
   fAmpLEDRec.AddAtAndExpand(gr,ipmt);
+  //should be change to real
+  Double_t xq[10] = { 411, 412,413,415,417,419,422,428,437,452};
+  Double_t yq[10] = {1,2,3,4,5,6,7,8,9,10};
+  TGraph* gr1 = new TGraph(10, xq, yq);
+  fAmpLED.AddAtAndExpand(gr1,ipmt);
   
 }
 
@@ -177,124 +194,93 @@ void AliT0CalibWalk::SetAmpLEDRec(Int_t ipmt)
 void AliT0CalibWalk::MakeWalkCorrGraph(const char *laserFile)
 {
   //make walk corerction for preprocessor
+  Int_t nmips=14;
+  gFile = TFile::Open(laserFile);
+  if(!gFile) {
+    AliError("No input laser data found ");
+  }
+  else
+    {
+      //      gFile->ls();
 
-  TFile *gFile = TFile::Open(laserFile);
-  //gSystem->Load("libSpectrum");
-  Int_t npeaks = 20;
-  Int_t sigma=3.;
-  Bool_t down = false;
-
-  Int_t index[20];
-  Char_t buf1[10], buf2[10],buf3[20];
-  
+      Float_t x1[14], y1[14]; 
+      Float_t x2[14], y2[14];
   
-  for (Int_t i=0; i<24; i++)
-    {
-      if(i>11) {  sprintf(buf1,"T0_A_%i_CFD",i+1-12); }
+      Float_t mips[14] = {1.,2,3,4,5,6,7,8,9,10,20,30,40,50};
+     
+      Float_t xx1[14],yy1[14], xx[14];
+      for (Int_t ii=0; ii<14; ii++)
+       {      x1[ii]=0; y1[ii]=0; x2[ii]=0; y2[ii]=0; }
+      
+      
+      TH2F*  hCFDvsQTC[24][14]; TH2F*  hCFDvsLED[24][14];
+      
+      for (Int_t i=0; i<24; i++)
+       {
+         for (Int_t im=0; im<nmips; im++)
+           {         
+             TString qtc = Form("CFDvsQTC%i_%i",i+1,im+1);
+             TString led = Form("CFDvsLED%i_%i",i+1,im+1);
+             hCFDvsQTC[i][im] = (TH2F*) gFile->Get(qtc.Data()) ;
+             hCFDvsLED[i][im] = (TH2F*) gFile->Get(led.Data());
+             
+             if(!hCFDvsQTC[i][im] || !hCFDvsLED[i][im]) 
+               AliWarning(Form(" no walk correction data in LASER DA for channel %i for amplitude %f MIPs",i,mips[im]));
+             
+             Float_t qtc0 = hCFDvsQTC[i][0]->GetMean(2);             
+             if(hCFDvsQTC[i][im])
+               {         
+                 x1[im] = hCFDvsQTC[i][im]->GetMean(1);
+                 y1[im] =  hCFDvsQTC[i][im]->GetMean(2) - qtc0;
+               }
+             Float_t led0 = hCFDvsLED[i][0]->GetMean(2);
+             if( hCFDvsLED[i][im]){
+                 x2[im] = hCFDvsLED[i][im]->GetMean(1);
+                 y2[im] =  hCFDvsLED[i][im]->GetMean(2) - led0;
+               }
+             xx[im]=mips[im];
+           }
+         for (Int_t imi=0; imi<nmips; imi++)
+           {
+             yy1[imi] = Float_t (mips[nmips-imi-1]);
+             xx1[imi] = x2[nmips-imi-1]; 
+           }
+         if(i==0) cout<<"Making graphs..."<<endl;
+         
+         TGraph *grwalkqtc = new TGraph (nmips,x1,y1);
+         TGraph *grwalkled = new TGraph (nmips,x2,y2);
+         fWalk.AddAtAndExpand(grwalkqtc,i);
+         fAmpLEDRec.AddAtAndExpand(grwalkled,i);
+        
+         //fit amplitude graphs to make comparison wth new one   
+         TGraph *grampled = new TGraph (nmips,xx1,yy1);
+         TGraph *grqtc = new TGraph (nmips,x1,xx);
+         fQTC.AddAtAndExpand(grqtc,i);  
+         fAmpLED.AddAtAndExpand(grampled,i);
 
-      if(i<12) { sprintf(buf1,"T0_C_%i_CFD",i+1); }
+         //fit amplitude graphs
+         Double_t parled [2], parqtc[2];
 
-      sprintf(buf2,"CFD_QTC%i",i+1);
-      sprintf(buf3,"CFD_LED%i",i+1);
-      cout<<buf1<<" "<<buf2<<" "<<buf3<<endl;
-      TH2F *qtc_cfd = (TH2F*) gFile->Get(buf2);
-      TH2F *led_cfd = (TH2F*) gFile->Get(buf3);
-      TH1F *cfd = (TH1F*) gFile->Get(buf1);
-      //       cfd->Draw();
-      TSpectrum *s = new TSpectrum(2*npeaks,1);
-      Int_t nfound = s->Search(cfd,sigma," ",0.05);
-      cout<<"Found "<<nfound<<" peaks sigma "<<sigma<<endl;;
-      if(nfound!=0){
-       Float_t *xpeak = s->GetPositionX();
-       TMath::Sort(nfound, xpeak, index,down);
-       Float_t xp = xpeak[index[0]];
-       Int_t xbin = cfd->GetXaxis()->FindBin(xp);
-       Float_t yp = cfd->GetBinContent(xbin);
-       cout<<"xbin = "<<xbin<<"\txpeak = "<<xpeak[1]<<"\typeak = "<<yp<<endl;
-       Float_t hmax = xp+10*sigma;
-       Float_t hmin = xp-10*sigma;
-       cout<<hmin<< " "<<hmax<<endl;
-       //          cfd->GetXaxis()->SetRange(hmin,hmax);
-       //          TF1 *g1 = new TF1("g1", "gaus", hmin, hmax);
-       // cfd->Fit("g1","R");
-       //      Int_t hminbin=  qtc_cfd->GetXaxis()->GetFirst();
-       //      Int_t hmaxbin=  qtc_cfd->GetXaxis()->GetLast();
-       Int_t nbins= qtc_cfd->GetXaxis()->GetNbins();
-       cout<<" nbins "<<nbins<<endl;
-       //      cout<<"  qtc_cfd "<<hminbin<<" "<<hmaxbin<<" "<<nbins<<endl;
-       //  qtc_cfd->Draw();
-       TProfile *pr_y = qtc_cfd->ProfileX();
-       //      Float_t maxHr=pr_y->GetBinCenter(hmaxbin);
-       
-       pr_y->SetMaximum(hmax);
-       pr_y->SetMinimum(hmin);
-       Int_t np=nbins/20;
-       Double_t *xx = new Double_t[np];
-       Double_t *yy = new Double_t[np];
-       Int_t ng=0;
-       Double_t yg=0;
-       for (Int_t ip=1; ip<nbins; ip++)
-         {
-           if(ip%20 != 0 ) {
-             if (pr_y->GetBinContent(ip) !=0)
-               yg +=pr_y->GetBinContent(ip);
-             //        cout<<ng<<" "<<pr_y->GetBinContent(ip)<<" "<<yg<<endl;
-             ng++;}
-           else {
-             xx[ip/20] = Float_t (pr_y->GetBinCenter(ip));
-             yy[ip/20] = yg/ng;
-             yg=0;
-             ng=0;
-             //                      cout<<ip<<" "<<ip/20<<" "<< xx[ip/20]<<" "<< yy[ip/20]<<endl;
-           }
+         TF1 *fuled = new TF1("fuled","expo",x1[0],x1[nmips]);
+         TF1 *fuqtc = new TF1 ("fuqtc","pol1",x1[0],x1[nmips]);
+         grqtc->Fit("fuqtc","Q"," ",x1[0],x1[nmips]);
+         grampled->Fit("fuled","Q"," ",xx1[0],xx1[nmips]);
+         fuled->GetParameters(&parled[0]);
+         fuqtc->GetParameters(&parqtc[0]);
+
+         for(Int_t ipar=0; ipar<2; ipar++) {
+           SetQTCpar(i,ipar,parqtc[ipar]);
+           SetAmpLEDpar(i,ipar,parled[ipar]);
+           //      cout<<" pars :::: "<<i<<" "<<ipar<<" qtc "<<parqtc[ipar]<<" led "<<parled[ipar]<<endl;
          }
-       TGraph *gr = new TGraph(np,xx,yy);
-       gr->SetMinimum(hmin);
-       gr->SetMaximum(hmax);
-       fWalk.AddAtAndExpand(gr,i);       
-       //      fWalk.AddAtAndExpand(gr,i+12);
 
-       Int_t nbinsled= led_cfd->GetXaxis()->GetNbins();
-       cout<<" nbins led "<<nbinsled<<endl;
-       //      cout<<"  qtc_cfd "<<hminbin<<" "<<hmaxbin<<" "<<nbins<<endl;
-       //  qtc_cfd->Draw();
-       TProfile *prled_y = led_cfd->ProfileX();
-       //      Float_t maxHr=pr_y->GetBinCenter(hmaxbin);
-       
-       prled_y->SetMaximum(hmax);
-       prled_y->SetMinimum(hmin);
-       Int_t npled=nbinsled/20;
-       Double_t *xxled = new Double_t[np];
-       Double_t *yyled = new Double_t[np];
-       Int_t ngled=0;
-       Double_t ygled=0;
-       for (Int_t ip=1; ip<nbinsled; ip++)
-         {
-           if(ip%20 != 0 ) {
-             if (prled_y->GetBinContent(ip) !=0)
-               ygled +=prled_y->GetBinContent(ip);
-             //        cout<<ng<<" "<<pr_y->GetBinContent(ip)<<" "<<yg<<endl;
-             ngled++;}
-           else {
-             xxled[ip/20] = Float_t (prled_y->GetBinCenter(ip));
-             yyled[ip/20] = ygled/ngled;
-             ygled=0;
-             ngled=0;
-           }
+         if(i==23){
+           cout<<"Graphs created..."<<endl;
          }
-       TGraph *grled = new TGraph(npled,xxled,yyled);
-       grled->SetMinimum(hmin);
-       grled->SetMaximum(hmax);
-       fAmpLEDRec.AddAtAndExpand(grled,i);       
-       //      delete [] xx;
-       //      delete [] yy;
-       
-      }
-      
-    }
-  
+       }
+    } //if gFile exits
+}
+
 
 
-}