]> git.uio.no Git - u/mrichter/AliRoot.git/blobdiff - TRD/AliTRDCalibraVdriftLinearFit.cxx
Revert "Revert "#103626: Commit DCal geometry to master" since the files are broken."
[u/mrichter/AliRoot.git] / TRD / AliTRDCalibraVdriftLinearFit.cxx
index a2a25815f016c881de5716d6d82000016848e510..855fceb5319c238fef97c25260da9814496f2df8 100644 (file)
 #include <TAxis.h>
 #include <TLinearFitter.h>
 #include <TMath.h>
+#include <TStyle.h>
+#include <TCanvas.h>
 #include <TTreeStream.h>
-
+#include <TGraphErrors.h>
+#include <TDirectory.h>
+#include <TTreeStream.h>
+#include <TF1.h>
 
 //header file
 #include "AliTRDCalibraVdriftLinearFit.h"
@@ -46,9 +51,19 @@ ClassImp(AliTRDCalibraVdriftLinearFit) /*FOLD00*/
 AliTRDCalibraVdriftLinearFit::AliTRDCalibraVdriftLinearFit() : /*FOLD00*/
   TObject(),
   fVersion(0),
+  fNameCalibUsed(""),
   fLinearFitterHistoArray(540),
   fLinearFitterPArray(540),
-  fLinearFitterEArray(540)
+  fLinearFitterEArray(540),
+  fRobustFit(kTRUE),
+  fMinNpointsFit(11),
+  fNbBindx(32),
+  fNbBindy(70),
+  fRangedx(0.8),
+  fRangedy(1.4),
+  fDebugStreamer(0x0),
+  fDebugLevel(0),
+  fSeeDetector(0)
 {
   //
   // default constructor
@@ -58,9 +73,19 @@ AliTRDCalibraVdriftLinearFit::AliTRDCalibraVdriftLinearFit() : /*FOLD00*/
 AliTRDCalibraVdriftLinearFit::AliTRDCalibraVdriftLinearFit(const AliTRDCalibraVdriftLinearFit &ped) : /*FOLD00*/
   TObject(ped),
   fVersion(ped.fVersion),
+  fNameCalibUsed(ped.fNameCalibUsed),
   fLinearFitterHistoArray(540),
   fLinearFitterPArray(540),
-  fLinearFitterEArray(540)
+  fLinearFitterEArray(540),
+  fRobustFit(kTRUE),
+  fMinNpointsFit(10),
+  fNbBindx(32),
+  fNbBindy(70),
+  fRangedx(0.8),
+  fRangedy(1.4),
+  fDebugStreamer(0x0),
+  fDebugLevel(0),
+  fSeeDetector(0)
 {
     //
     // copy constructor
@@ -84,9 +109,19 @@ AliTRDCalibraVdriftLinearFit::AliTRDCalibraVdriftLinearFit(const AliTRDCalibraVd
 AliTRDCalibraVdriftLinearFit::AliTRDCalibraVdriftLinearFit(const TObjArray &obja) : /*FOLD00*/
   TObject(),
   fVersion(0),
+  fNameCalibUsed(""),
   fLinearFitterHistoArray(540),
   fLinearFitterPArray(540),
-  fLinearFitterEArray(540)
+  fLinearFitterEArray(540),
+  fRobustFit(kTRUE),
+  fMinNpointsFit(10),
+  fNbBindx(32),
+  fNbBindy(70),
+  fRangedx(0.8),
+  fRangedy(1.4),
+  fDebugStreamer(0x0),
+  fDebugLevel(0),
+  fSeeDetector(0)
 {
   //
   // constructor from a TObjArray
@@ -125,6 +160,8 @@ AliTRDCalibraVdriftLinearFit::~AliTRDCalibraVdriftLinearFit() /*FOLD00*/
   fLinearFitterPArray.Delete();
   fLinearFitterEArray.Delete();
 
+  if ( fDebugStreamer ) delete fDebugStreamer;
+
 }
 //_____________________________________________________________________________
 void AliTRDCalibraVdriftLinearFit::Copy(TObject &c) const
@@ -205,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)
 {
     //
@@ -236,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");
@@ -353,7 +410,7 @@ void AliTRDCalibraVdriftLinearFit::FillPEArray()
        TVectorD   pare = TVectorD(2);
        TVectorD  *parE = new TVectorD(3);
        //printf("Fit\n");
-       if((linearfitter.EvalRobust(0.8)==0)) {
+       if(((fRobustFit) && (linearfitter.EvalRobust(0.8)==0)) || ((!fRobustFit) && (linearfitter.Eval()==0))) {
          //if((linearfitter.Eval()==0)) {
          //printf("Take the param\n");
          linearfitter.GetParameters(*par);
@@ -385,4 +442,103 @@ void AliTRDCalibraVdriftLinearFit::FillPEArray()
    
 }
 
+//____________Functions fit Online CH2d________________________________________
+void AliTRDCalibraVdriftLinearFit::FillPEArray2()
+{
+  //
+  // Fill fFitterPArray and fFitterEArray from inside
+  //
 
+  // Loop over histos 
+  TF1 *f1 = new TF1("f1","[0]+[1]*x",-1,1);
+      
+  for(Int_t cb = 0; cb < 540; cb++){
+    //const TH2S *fitterhisto = (TH2S*)fLinearFitterHistoArray.UncheckedAt(cb);
+    TH2S *fitterhisto = (TH2S*)fLinearFitterHistoArray.UncheckedAt(cb);
+    //printf("Processing the detector cb %d we find %d\n",cb, (Bool_t) fitterhisto);    
+  
+    if ( fitterhisto != 0 ){
+      
+      Int_t nEntries=0;
+      TGraphErrors *gg=DrawMS(fitterhisto,nEntries);
+      if (!gg) continue;
+      // Number of points of the TGraphErrors
+      //printf("Number of points %d for detector %d\n",gg->GetN(),cb);
+      if(gg->GetN() <  fMinNpointsFit) {
+       if(gg) delete gg;
+       continue;       
+      }
+      //printf("det %d, number of points %d, nEntries %d\n",cb,gg->GetN(),nEntries);
+      gg->Fit(f1,"Q0");
+      
+      TVectorD  *par  = new TVectorD(2);
+      TVectorD  *parE = new TVectorD(3);
+      (*parE)[0] = 0.0;
+      (*parE)[1] = 0.0;
+      (*parE)[2] = (Double_t) nEntries;
+      (*par)[0] = f1->GetParameter(0);
+      (*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;
+      }
+      else {
+       if(cb==fSeeDetector) {
+         gStyle->SetPalette(1);
+         gStyle->SetOptStat(1111);
+         gStyle->SetPadBorderMode(0);
+         gStyle->SetCanvasColor(10);
+         gStyle->SetPadLeftMargin(0.13);
+         gStyle->SetPadRightMargin(0.01);
+         TCanvas *stat = new TCanvas("stat","",50,50,600,800);
+         stat->cd(1);
+         fitterhisto->Draw("colz");
+         gg->Draw("P");
+         f1->Draw("same");
+         break;
+       }
+      } 
+    }
+  }
+  if(fDebugLevel==0) delete f1;
+}
+
+//_________Helper function__________________________________________________
+TGraphErrors* AliTRDCalibraVdriftLinearFit::DrawMS(const TH2 *const h2, Int_t &nEntries)
+{
+  //
+  // Debug function
+  //
+
+  TF1 fg("fg", "gaus", -10., 30.);
+  TGraphErrors *gp = new TGraphErrors();
+
+  TAxis *ax(h2->GetXaxis());
+  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);
+    fg.SetParameter(1, h1->GetMean());
+    //Float_t x(ax->GetBinCenter(jpt)); 
+    Int_t n=Int_t(h1->Integral(1, h1->GetNbinsX()));
+    nEntries+=n;
+    if(n < 15){
+      //Warning("drawMS()", Form("reject x[%d]=%f on n=%d", jpt, x, n));
+      continue;
+    }
+    h1->Fit(&fg, "WWQ0");
+    if(fg.GetNDF()<2){
+      //Warning("drawMS()", Form("reject x[%d]=%f on NDF=%d", jpt, x, fg.GetNDF()));
+      continue;
+    }
+    if(((fg.GetParameter(1)+fg.GetParameter(2)/2)>ay->GetXmax()) || ((fg.GetParameter(1)-fg.GetParameter(2)/2)<ay->GetXmin()) || (TMath::Abs(fg.GetParameter(0))< 0.00001)) continue;
+    gp->SetPoint(ig, ax->GetBinCenter(jpt), fg.GetParameter(1));
+    gp->SetPointError(ig, 0, TMath::Sqrt(pow(fg.GetParError(1),2) + (1/pow(fg.GetParameter(0),2))));
+    ig++;
+  }
+  delete h1;
+  return gp;
+}