warning fixed
authoralla <alla@f7af4fe6-9843-0410-8265-dc069ae4e863>
Wed, 16 Nov 2011 19:11:26 +0000 (19:11 +0000)
committeralla <alla@f7af4fe6-9843-0410-8265-dc069ae4e863>
Wed, 16 Nov 2011 19:11:26 +0000 (19:11 +0000)
T0/AliT0CalibWalk.cxx

index 4f2aa5e..e0fd56b 100644 (file)
@@ -33,7 +33,6 @@
 #include <Riostream.h>
 #include <TSpectrum.h>
 #include <TProfile.h>
-#include <TCanvas.h>
 #include <TF1.h>
 
 ClassImp(AliT0CalibWalk)
@@ -98,33 +97,9 @@ AliT0CalibWalk::~AliT0CalibWalk()
 Bool_t AliT0CalibWalk::MakeWalkCorrGraph(const char *laserFile)
 {
   //make walk corerction for preprocessor
-  Int_t npeaks = 20;
-  Int_t sigma=3;
-  Bool_t down=false;
-  Int_t index[20];
+  Float_t sigma,cfdmean, qtmean, ledmean;
   Bool_t ok=true;
-  Float_t   mips[100];
-  Int_t nfound=0;
-  /*  
-  //061110 
- Float_t timefitth[24] = { 5594, 5618, 5630, 5662,  5623, 5645,
-                           5645, 5602, 5691, 5669,  5636, 5602,
-                           5651, 5654, 5612, 5625,  5669, 5678,
-                           5652, 5673, 5565, 5717,  5692, 5738} ;
-  */
-  //051210
-  Float_t timefitth[24] = { 5587, 5609, 5621, 5653, 5615, 5635, 
-                           5637, 5595, 5683, 5659, 5629, 5601,
-                           5643, 5654, 5601, 5615, 5661, 5667, 
-                           5645, 5665, 5663, 5709, 5685, 5725 };
-  Float_t xdata[4000], ydata[4000];
-  Int_t npoints =0;
-  TH1F* www[24];
-  TFile *ff = new TFile("calibwalk0.root");
-
-  Int_t ipmt, jpoint; 
+  Float_t   mips[50];
   
   gFile = TFile::Open(laserFile);
   if(!gFile) {
@@ -135,7 +110,7 @@ Bool_t AliT0CalibWalk::MakeWalkCorrGraph(const char *laserFile)
       //      gFile->ls();
       TH1F* hAmp = (TH1F*) gFile->Get("hAmpLaser");
       Int_t nmips=0;
-      for (Int_t ibin=0; ibin<1200; ibin++) { 
+      for (Int_t ibin=0; ibin<2000; ibin++) {
        Float_t bincont = hAmp->GetBinContent(ibin);
        if(bincont>0){ 
          mips[nmips] = hAmp->GetXaxis()->GetBinCenter(ibin);
@@ -144,33 +119,25 @@ Bool_t AliT0CalibWalk::MakeWalkCorrGraph(const char *laserFile)
        }       
       }    
       
-      /*      if (nmips<17) {
+      if (nmips<17) {
        ok=false;
        return ok;
-       } */
-      //  TCanvas *c1 = new TCanvas("c1", "LED-CFD C side",0,48,1280,951);
-      // c1->Divide(4,3);
-   //   TCanvas *c2 = new TCanvas("c1", "LED-CFD C side",0,48,1280,951);
-   //   c2->Divide(4,3);
-        
+      } 
+       
       Float_t x1[50], y1[50]; 
       Float_t x2[50], xx2[50],y2[50];
       Float_t xx1[50],yy1[50], xx[50];
       
-      Float_t cfd0[24];
-      Int_t startim = 0;
+      Float_t cfd0 = 0;
       
+      for (Int_t ii=0; ii<nmips; ii++)
+       x1[ii] = y1[ii] = x2[ii] = y2[ii] = 0; 
       
       for (Int_t i=0; i<24; i++)
        {
-         //     if(i==13) continue;
-         //      else c2->cd(i+1-12);
-         for (Int_t im=startim; im<nmips; im++)
+         cfd0 = 0;
+         for (Int_t im=0; im<nmips; im++)
            {         
-
-             x1[im] = y1[im] = x2[im] = y2[im] = 0; 
-             if(i==23 && im<2) continue;
-             xx[im]=mips[im];
              TString cfd = Form("hCFD%i_%i",i+1,im+1);
              TString qtc = Form("hQTC%i_%i",i+1,im+1);
              TString led = Form("hLED%i_%i",i+1,im+1);
@@ -178,110 +145,60 @@ Bool_t AliT0CalibWalk::MakeWalkCorrGraph(const char *laserFile)
              TH1F *hCFD = (TH1F*) gFile->Get(cfd.Data()) ;
              TH1F *hLED = (TH1F*) gFile->Get(led.Data());
              TH1F *hQTC = (TH1F*) gFile->Get(qtc.Data()) ;
+             //              hCFD->SetDirectory(0);
+             //              hQTC->SetDirectory(0);
+             //              hLED->SetDirectory(0);
              if(!hCFD )
                AliWarning(Form(" no CFD data in LASER DA for channel %i for amplitude %f MIPs",i,mips[im]));
              if(!hQTC )
                AliWarning(Form(" no QTC correction data in LASER DA for channel %i for amplitude %f MIPs",i,mips[im]));
              if(!hLED)       
                AliWarning(Form(" no LED correction data in LASER DA for channel %i for amplitude %f MIPs",i,mips[im]));
-             
-             if(hCFD ) {
-               hCFD->SetDirectory(0);
-               hCFD->GetXaxis()->SetRangeUser(timefitth[i]-100,timefitth[i]+100);
-
-               TF1 *fg = new TF1("fg","gaus",timefitth[i]-50,timefitth[i]+50);
-               hCFD->Fit("fg","RQ", " ",timefitth[i]-50,timefitth[i]+50);
-               Double_t par[3];
-               fg->GetParameters(&par[0]);
-
-               TSpectrum *s = new TSpectrum(2*npeaks,1);
-               nfound = s->Search(hCFD,sigma," ",0.1);
-               if(nfound!=0){
-                 Float_t *xpeak = s->GetPositionX();
-                 TMath::Sort(nfound, xpeak, index,down);
-                 Float_t xp = xpeak[index[0]];
-                 Double_t hmax = xp+3*sigma;
-                 Double_t hmin = xp-3*sigma;
-                 hCFD->GetXaxis()->SetRangeUser(hmin-10,hmax+10);
-                 //    if (i==23) 
-                 //              cout<<"PMT "<<i<<" MIPS "<<mips[im]<<" mean "<<hCFD->GetMean()<<
-                 //     " xp "<<xp<<" fit "<<par[1]<<endl;
-                 //              cout<<"PMT "<<i<<" MIPS "<<im<<" xp "<<xp<<endl;
-               }
+             if( hCFD && hCFD->GetEntries()<500 ) {
+               ok=false;
+               printf("no peak in CFD spectrum for PMT %i amplitude %i\n",i,im);
+               return ok;
+             }
+             if(hCFD && hCFD->GetEntries()>500 ) {
+               if( hCFD->GetRMS() >= 1.5) 
+                 GetMeanAndSigma(hCFD, cfdmean, sigma);
                else
-                 {
-                   hCFD->Rebin(2);
-                   TSpectrum *s1 = new TSpectrum(2*npeaks,1);
-                   nfound = s1->Search(hCFD,sigma," ",0.1);
-                   if(nfound!=0){
-                     Float_t *xpeak = s1->GetPositionX();
-                     TMath::Sort(nfound, xpeak, index,down);
-                     Float_t xp = xpeak[index[0]];
-                     Double_t hmax = xp+3*sigma;
-                     Double_t hmin = xp-3*sigma;
-                     hCFD->GetXaxis()->SetRangeUser(hmin-10,hmax+10);
-                   }
-                   else 
-                     {
-                       ok=true;
-                       //                      ok=false;
-                       printf("no peak in CFD spectrum for PMT %i amplitude %i\n",i,im);
-                       return ok;
-                     }
+                 cfdmean = hCFD->GetMean();
                
-                 }
-               // } //if no C12 80MIPs
-               //        if (i<12) cout<<"!!!! PMT "<<i<<" MIPS "<<im<<" fit  "<<par[1]<<
-               //" mean "<<hCFD->GetMean()<<endl;
-                if (im == 0) cfd0[i] = par[1];
-               if(i==23 && im==2) cfd0[i] = par[1]; 
-                  y1[im] =  par[1] - cfd0[i];
-                  if ((i == 5 && im == 15) ||(i == 8 && im == 19) || (i==11 && im==0) || i==23)
-                      y1[im] =  hCFD->GetMean() - cfd0[i];
-                 
-               //              if (im == 0) cfd0[i] = hCFD->GetMean();
-               //              y1[im] =  hCFD->GetMean() - cfd0[i];
-                 
-                  //   if (i==23 && im<2) y1[im]=0;
-               /*              if( i == 11 && im==10 ){        y1[im] = y1[im-1];
-               cout<<im<<" "<<mips[im]<<" mip "<<" pmt "<<i<<" mean "<<hCFD->GetMean()<<" walk "<<     y1[im]<<endl;
-               }*/ 
-
-             }
-             if( hQTC) {
-               x1[im] = hQTC->GetMean();
+               Int_t   maxBin = hCFD->GetMaximumBin(); 
+               Double_t  meanEstimate = hCFD->GetBinCenter( maxBin); 
+               if(TMath::Abs(meanEstimate - cfdmean) > 20 ) cfdmean = meanEstimate; 
+               if (im == 0) cfd0 = cfdmean;
+               y1[im] =  cfdmean - cfd0;
+             } 
+             if(hQTC && hQTC->GetEntries()>500) {
+               GetMeanAndSigma(hQTC, qtmean, sigma);
+               
+               x1[im] = qtmean;
                if( x1[im] == 0) {
                  ok=false;
                  printf("no peak in QTC signal for PMT %i amplitude %i\n",i,im);
                  return ok;
-                 }
-             }
-             
-             if( hLED){
-               hLED->SetDirectory(0);
-               TSpectrum *s = new TSpectrum(2*npeaks,1);
-               nfound = s->Search(hLED,sigma," ",0.1);
-               if(nfound!=0){
-                 Float_t *xpeak = s->GetPositionX();
-                 TMath::Sort(nfound, xpeak, index,down);
-                 Float_t xp = xpeak[index[0]];
-                 Double_t hmax = xp+10*sigma;
-                 Double_t hmin = xp-10*sigma;
-                 hLED->GetXaxis()->SetRangeUser(hmin-10,hmax+10);
                }
-               else
-                 { 
-                   ok=false;
-                   printf("no peak in LED spectrum for PMT %i amplitude %i\n",i,im);
-                   return ok;
-                 }
-               x2[im] = hLED->GetMean();
-               xx2[im] = x2[nmips-im-1];
              }
-             
-             if (hQTC) delete  hQTC;
-             if (hCFD) delete  hCFD;
-             if (hLED) delete  hLED;
+               
+             if( hLED && hLED->GetEntries()>500) {
+               GetMeanAndSigma(hLED, ledmean, sigma);
+             }                            
+             else
+               { 
+                 ok=false;
+                 printf("no peak in LED spectrum for PMT %i amplitude %i\n",i,im);
+                 return ok;
+               }
+             x2[im] = ledmean;
+             xx2[im] = x2[nmips-im-1];
+               
+               xx[im]=mips[im];
+               if (hQTC) delete  hQTC;
+               if (hCFD) delete  hCFD;
+               if (hLED) delete  hLED;
+               
            }
          
          for (Int_t imi=0; imi<nmips; imi++)
@@ -292,49 +209,30 @@ Bool_t AliT0CalibWalk::MakeWalkCorrGraph(const char *laserFile)
            }
          
          if(i==0) cout<<"Making graphs..."<<endl;
-         // TGraph *grwalkqtc = new TGraph (nmips-4,x1,y1);
-         // grwalkqtc->SetTitle(Form("PMT%i",i));
-         TGraph *grwalkled = new TGraph (nmips-4,x2,y1);
-         grwalkled->SetTitle(Form("PMT%i",i));
-         // fWalk.AddAtAndExpand(grwalkqtc,i);
-         fAmpLEDRec.AddAtAndExpand(grwalkled,i);
-         //      cout<<" add walk "<<i<<endl;
-         
-         //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);
-         //      cout<<" add amp "<<i<<endl;
-         
-         if(i==23)
-           cout<<"Graphs created..."<<endl;
-         //       if(i==13) continue;
-         
-         www[i] = (TH1F*)ff->Get(Form("hy%i",i) );
-         www[i]->Print();
-         npoints =0;
-         for(Int_t ibin=1; ibin< (www[i]->GetNbinsX() -2); ibin++) {
-            if( www[i]->GetBinContent(ibin)>0 && 
-               www[i]->GetBinContent(ibin) !=1000 ) {
-             ydata[ibin-1] = www[i]->GetBinContent(ibin) - 1000;
-             xdata[ibin-1] = ibin+799;
-             cout<<ibin<<" "<<xdata[ibin]<<" "<<ydata[ibin]<<endl;
-             npoints++;
-            }
-         } 
-
-
-         TGraph *grwalkqtc = new TGraph (npoints-3,xdata,ydata);
-         grwalkqtc->SetTitle(Form("PMT%i",i));
-         fWalk.AddAtAndExpand(grwalkqtc,i);
-       }
          
+      TGraph *grwalkqtc = new TGraph (nmips,x1,y1);
+      grwalkqtc->SetTitle(Form("PMT%i",i));
+      TGraph *grwalkled = new TGraph (nmips,x2,y1);
+      grwalkled->SetTitle(Form("PMT%i",i));
+      fWalk.AddAtAndExpand(grwalkqtc,i);
+      fAmpLEDRec.AddAtAndExpand(grwalkled,i);
+      //         cout<<" add walk "<<i<<endl;
       
+      //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);
+      //         cout<<" add amp "<<i<<endl;
+      
+      if(i==23)
+       cout<<"Graphs created..."<<endl;   
+       }
     } //if gFile exits
 
   return ok;
 }
+
 void AliT0CalibWalk::GetMeanAndSigma(TH1F* hist, Float_t &mean, Float_t &sigma) 
 {