some coverity defects fixed

[u/mrichter/AliRoot.git] / T0 / AliT0CalibWalk.cxx

diff --git a/T0/AliT0CalibWalk.cxx b/T0/AliT0CalibWalk.cxx
index 1dca2c32e631f33537076d1cf2a4a866d8dc12dd..0d1a6566934f10c07e40a39623a1400a019a005b 100644 (file)
--- a/T0/AliT0CalibWalk.cxx
+++ b/T0/AliT0CalibWalk.cxx
@@ -97,13 +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, xp, cfdmean, qtmean, ledmean;
   Bool_t ok=true;
-  Float_t   mips[20];
-  Int_t nfound=0;
+  Float_t   mips[50];
   
   gFile = TFile::Open(laserFile);
   if(!gFile) {
@@ -149,86 +145,59 @@ 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);
+             //              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 ) {
-               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);
-               }
-               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=false;
-                       printf("no peak in CFD spectrum for PMT %i amplitude %i\n",i,im);
-                       return ok;
-                     }
-
-                 }
-
-               if (im == 0) cfd0[i] = hCFD->GetMean();
-               y1[im] =  hCFD->GetMean() - cfd0[i];
-             }
-             if( hQTC) {
-               x1[im] = hQTC->GetMean();
-               if( x1[im] == 0) {
+             
+             if(hCFD && hCFD->GetEntries()>500) {
+               GetMeanAndSigma(hCFD, xp, sigma);
+               }                  
+             else 
+               {
                  ok=false;
-                 printf("no peak in QTC signal for PMT %i amplitude %i\n",i,im);
+                 printf("no peak in CFD spectrum for PMT %i amplitude %i\n",i,im);
                  return ok;
+               }
+
+             if (hCFD->GetRMS() < 1.5) cfdmean = hCFD->GetMean();
+             else cfdmean=xp;
+             if (im == 0) cfd0[i] = cfdmean;
+               y1[im] =  cfdmean - cfd0[i];
+               
+               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){
-               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);
                }
+               
+               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] = hLED->GetMean();
+               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;
+               
+               xx[im]=mips[im];
+               if (hQTC) delete  hQTC;
+               if (hCFD) delete  hCFD;
+               if (hLED) delete  hLED;
+               
            }
          
          for (Int_t imi=0; imi<nmips; imi++)
@@ -282,6 +251,27 @@ Bool_t AliT0CalibWalk::MakeWalkCorrGraph(const char *laserFile)
   return ok;
 }
 
+void AliT0CalibWalk::GetMeanAndSigma(TH1F* hist, Float_t &mean, Float_t &sigma) 
+{
+
+  const double window = 5.;  //fit window 
+  double norm  = hist->Integral();  // normalize to one count
+  hist->Scale(1./norm); 
+ 
+  double meanEstimate, sigmaEstimate; 
+  int maxBin;
+  maxBin        =  hist->GetMaximumBin(); //position of maximum
+  meanEstimate  =  hist->GetBinCenter( maxBin); // mean of gaussian sitting in maximum
+  sigmaEstimate = hist->GetRMS();
+  TF1* fit= new TF1("fit","gaus", meanEstimate - window*sigmaEstimate, meanEstimate + window*sigmaEstimate);
+  fit->SetParameters(hist->GetBinContent(maxBin), meanEstimate, sigmaEstimate);
+  hist->Fit("fit","RQ","Q");
+
+  mean  = (Float_t) fit->GetParameter(1);
+  sigma = (Float_t) fit->GetParameter(2);
+
+  delete fit;
+}