Option to fit with chi2 or likelihood

[u/mrichter/AliRoot.git] / TRD / AliTRDCalibraVdriftLinearFit.cxx

diff --git a/TRD/AliTRDCalibraVdriftLinearFit.cxx b/TRD/AliTRDCalibraVdriftLinearFit.cxx
index c3d5c5dd02eddcd57854a5472caf15266d88b72f..a57524d55f739fd3c9a3afe850c33d48a8483c43 100644 (file)
--- a/TRD/AliTRDCalibraVdriftLinearFit.cxx
+++ b/TRD/AliTRDCalibraVdriftLinearFit.cxx
@@ -56,6 +56,11 @@ AliTRDCalibraVdriftLinearFit::AliTRDCalibraVdriftLinearFit() : /*FOLD00*/
   fLinearFitterPArray(540),
   fLinearFitterEArray(540),
   fRobustFit(kTRUE),
+  fMinNpointsFit(11),
+  fNbBindx(32),
+  fNbBindy(70),
+  fRangedx(0.8),
+  fRangedy(1.4),
   fDebugStreamer(0x0),
   fDebugLevel(0),
   fSeeDetector(0)
@@ -73,6 +78,11 @@ AliTRDCalibraVdriftLinearFit::AliTRDCalibraVdriftLinearFit(const AliTRDCalibraVd
   fLinearFitterPArray(540),
   fLinearFitterEArray(540),
   fRobustFit(kTRUE),
+  fMinNpointsFit(10),
+  fNbBindx(32),
+  fNbBindy(70),
+  fRangedx(0.8),
+  fRangedy(1.4),
   fDebugStreamer(0x0),
   fDebugLevel(0),
   fSeeDetector(0)
@@ -104,6 +114,11 @@ AliTRDCalibraVdriftLinearFit::AliTRDCalibraVdriftLinearFit(const TObjArray &obja
   fLinearFitterPArray(540),
   fLinearFitterEArray(540),
   fRobustFit(kTRUE),
+  fMinNpointsFit(10),
+  fNbBindx(32),
+  fNbBindy(70),
+  fRangedx(0.8),
+  fRangedy(1.4),
   fDebugStreamer(0x0),
   fDebugLevel(0),
   fSeeDetector(0)
@@ -227,6 +242,26 @@ void AliTRDCalibraVdriftLinearFit::Add(const AliTRDCalibraVdriftLinearFit *ped)
   }
 }
 //______________________________________________________________________________________
+TH2S* AliTRDCalibraVdriftLinearFit::AddAll() 
+{
+    //
+    // return pointer to TH2S of all added histos 
+  //
+
+  TH2S *histo = 0x0;
+  for(Int_t k=0; k < 540; k++) {
+    TH2S * u = GetLinearFitterHistoForce(k);
+    if(k == 0) {
+      histo = new TH2S(*u);
+      histo->SetName("sum");
+    }
+    else histo->Add(u);
+  }
+
+  return histo;
+
+}
+//______________________________________________________________________________________
 TH2S* AliTRDCalibraVdriftLinearFit::GetLinearFitterHisto(Int_t detector, Bool_t force)
 {
     //
@@ -258,8 +293,8 @@ TH2S* AliTRDCalibraVdriftLinearFit::GetLinearFitterHistoForce(Int_t detector)
   name +=  fVersion;
   
   TH2S *lfdv = new TH2S((const Char_t *)name,(const Char_t *) name
-                       ,36,-0.9,0.9,48
-                       ,-1.2,1.2);
+                       ,(Int_t)fNbBindx,-fRangedx,fRangedx,(Int_t)fNbBindy
+                       ,-fRangedy,fRangedy);
   lfdv->SetXTitle("tan(phi_{track})");
   lfdv->SetYTitle("dy/dt");
   lfdv->SetZTitle("Number of clusters");
@@ -336,8 +371,8 @@ void AliTRDCalibraVdriftLinearFit::FillPEArray()
     if ( linearfitterhisto != 0 ){
       
       // Fill a linearfitter
-      TAxis *xaxis = linearfitterhisto->GetXaxis();
-      TAxis *yaxis = linearfitterhisto->GetYaxis();
+      const TAxis *xaxis = linearfitterhisto->GetXaxis();
+      const TAxis *yaxis = linearfitterhisto->GetYaxis();
       TLinearFitter linearfitter = TLinearFitter(2,"pol1");
       //printf("test\n");
       Double_t integral = linearfitterhisto->Integral();
@@ -428,7 +463,8 @@ void AliTRDCalibraVdriftLinearFit::FillPEArray2()
       TGraphErrors *gg=DrawMS(fitterhisto,nEntries);
       if (!gg) continue;
       // Number of points of the TGraphErrors
-      if(gg->GetN() < 20) {
+      //printf("Number of points %d for detector %d\n",gg->GetN(),cb);
+      if(gg->GetN() <  fMinNpointsFit) {
        if(gg) delete gg;
        continue;       
       }
@@ -444,6 +480,7 @@ void AliTRDCalibraVdriftLinearFit::FillPEArray2()
       (*par)[1] = f1->GetParameter(1);
       fLinearFitterPArray.AddAt(par,cb);
       fLinearFitterEArray.AddAt(parE,cb);
+      //printf("Value %f for detector %d\n",(*par)[1],cb);
             
       if(fDebugLevel==0) {
        if(gg) delete gg;
@@ -479,8 +516,8 @@ TGraphErrors* AliTRDCalibraVdriftLinearFit::DrawMS(const TH2 *const h2, Int_t &n
   TF1 fg("fg", "gaus", -10., 30.);
   TGraphErrors *gp = new TGraphErrors();
 
-  TAxis *ax(h2->GetXaxis());
-  TAxis *ay(h2->GetYaxis());
+  const TAxis *ax(h2->GetXaxis());
+  const TAxis *ay(h2->GetYaxis());
   TH1D *h1(NULL);
   for(Int_t ipt(0), jpt(1), ig(0); ipt<ax->GetNbins(); ipt++, jpt++){
     h1 = h2->ProjectionY("py", jpt, jpt);