#include <TAxis.h>
#include <TLinearFitter.h>
#include <TMath.h>
+#include <TStyle.h>
+#include <TCanvas.h>
#include <TTreeStream.h>
-#include <THnSparse.h>
-
+#include <TGraphErrors.h>
+#include <TDirectory.h>
+#include <TTreeStream.h>
+#include <TF1.h>
//header file
#include "AliTRDCalibraVdriftLinearFit.h"
AliTRDCalibraVdriftLinearFit::AliTRDCalibraVdriftLinearFit() : /*FOLD00*/
TObject(),
fVersion(0),
+ fNameCalibUsed(""),
fLinearFitterHistoArray(540),
fLinearFitterPArray(540),
- fLinearFitterEArray(540)
+ fLinearFitterEArray(540),
+ fRobustFit(kTRUE),
+ fMinNpointsFit(11),
+ fDebugStreamer(0x0),
+ fDebugLevel(0),
+ fSeeDetector(0)
{
//
// default constructor
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),
+ fDebugStreamer(0x0),
+ fDebugLevel(0),
+ fSeeDetector(0)
{
//
// copy constructor
const TVectorD *vectorE = (TVectorD*)ped.fLinearFitterEArray.UncheckedAt(idet);
const TVectorD *vectorP = (TVectorD*)ped.fLinearFitterPArray.UncheckedAt(idet);
- const THnSparseS *hped = (THnSparseS*)ped.fLinearFitterHistoArray.UncheckedAt(idet);
+ const TH2S *hped = (TH2S*)ped.fLinearFitterHistoArray.UncheckedAt(idet);
if ( vectorE != 0x0 ) fLinearFitterEArray.AddAt(new TVectorD(*vectorE), idet);
if ( vectorP != 0x0 ) fLinearFitterPArray.AddAt(new TVectorD(*vectorP), idet);
if ( hped != 0x0 ){
- THnSparseS *hNew = (THnSparseS *)hped->Clone();
+ TH2S *hNew = (TH2S *)hped->Clone();
//hNew->SetDirectory(0);
fLinearFitterHistoArray.AddAt(hNew,idet);
}
AliTRDCalibraVdriftLinearFit::AliTRDCalibraVdriftLinearFit(const TObjArray &obja) : /*FOLD00*/
TObject(),
fVersion(0),
+ fNameCalibUsed(""),
fLinearFitterHistoArray(540),
fLinearFitterPArray(540),
- fLinearFitterEArray(540)
+ fLinearFitterEArray(540),
+ fRobustFit(kTRUE),
+ fMinNpointsFit(10),
+ fDebugStreamer(0x0),
+ fDebugLevel(0),
+ fSeeDetector(0)
{
//
// constructor from a TObjArray
//
for (Int_t idet = 0; idet < 540; idet++){
- const THnSparseS *hped = (THnSparseS*)obja.UncheckedAt(idet);
+ const TH2S *hped = (TH2S*)obja.UncheckedAt(idet);
if ( hped != 0x0 ){
- THnSparseS *hNew = (THnSparseS *)hped->Clone();
+ TH2S *hNew = (TH2S *)hped->Clone();
//hNew->SetDirectory(0);
fLinearFitterHistoArray.AddAt(hNew,idet);
}
//
// destructor
//
+ fLinearFitterHistoArray.SetOwner();
+ fLinearFitterPArray.SetOwner();
+ fLinearFitterEArray.SetOwner();
+
fLinearFitterHistoArray.Delete();
fLinearFitterPArray.Delete();
fLinearFitterEArray.Delete();
+ if ( fDebugStreamer ) delete fDebugStreamer;
+
}
//_____________________________________________________________________________
void AliTRDCalibraVdriftLinearFit::Copy(TObject &c) const
// Copy only the histos
for (Int_t idet = 0; idet < 540; idet++){
if(fLinearFitterHistoArray.UncheckedAt(idet)){
- THnSparseS *hped1 = (THnSparseS *)target.GetLinearFitterHisto(idet,kTRUE);
+ TH2S *hped1 = (TH2S *)target.GetLinearFitterHisto(idet,kTRUE);
//hped1->SetDirectory(0);
- hped1->Add((const THnSparseS *)fLinearFitterHistoArray.UncheckedAt(idet));
+ hped1->Add((const TH2S *)fLinearFitterHistoArray.UncheckedAt(idet));
}
}
// Copy only the histos
for (Int_t idet = 0; idet < 540; idet++){
if(entry->GetLinearFitterHisto(idet)){
- THnSparseS *hped1 = (THnSparseS *)GetLinearFitterHisto(idet,kTRUE);
+ TH2S *hped1 = (TH2S *)GetLinearFitterHisto(idet,kTRUE);
Double_t entriesa = hped1->GetEntries();
- Double_t entriesb = ((THnSparseS *)entry->GetLinearFitterHisto(idet))->GetEntries();
+ Double_t entriesb = ((TH2S *)entry->GetLinearFitterHisto(idet))->GetEntries();
if((entriesa + entriesb) < 5*32767) hped1->Add(entry->GetLinearFitterHisto(idet));
}
}
return count;
}
//_____________________________________________________________________
-void AliTRDCalibraVdriftLinearFit::Add(AliTRDCalibraVdriftLinearFit *ped)
+void AliTRDCalibraVdriftLinearFit::Add(const AliTRDCalibraVdriftLinearFit *ped)
{
//
// Add histo
fVersion++;
for (Int_t idet = 0; idet < 540; idet++){
- const THnSparseS *hped = (THnSparseS*)ped->GetLinearFitterHisto(idet);
+ const TH2S *hped = (TH2S*)ped->GetLinearFitterHistoNoForce(idet);
//printf("idet %d\n",idet);
if ( hped != 0x0 ){
//printf("add\n");
- THnSparseS *hped1 = (THnSparseS *)GetLinearFitterHisto(idet,kTRUE);
+ TH2S *hped1 = (TH2S *)GetLinearFitterHisto(idet,kTRUE);
Double_t entriesa = hped1->GetEntries();
Double_t entriesb = hped->GetEntries();
if((entriesa + entriesb) < 5*32767) hped1->Add(hped);
}
}
//______________________________________________________________________________________
-THnSparse* AliTRDCalibraVdriftLinearFit::GetLinearFitterHisto(Int_t detector, Bool_t force)
+TH2S* AliTRDCalibraVdriftLinearFit::GetLinearFitterHisto(Int_t detector, Bool_t force)
{
//
// return pointer to TH2F histo
// if force is true create a new histo if it doesn't exist allready
//
if ( !force || fLinearFitterHistoArray.UncheckedAt(detector) )
- return (THnSparse*)fLinearFitterHistoArray.UncheckedAt(detector);
-
- // if we are forced and TLinearFitter doesn't yes exist create it
-
- // new TH2F
- TString name("LFDV");
- name += detector;
- name += "version";
- name += fVersion;
-
- //create the map
- Int_t thnDim[2];
- thnDim[0] = 36;
- thnDim[1] = 48;
+ return (TH2S*)fLinearFitterHistoArray.UncheckedAt(detector);
- //arrays for lower bounds :
- Double_t* binEdges[2];
- for(Int_t ivar = 0; ivar < 2; ivar++)
- binEdges[ivar] = new Double_t[thnDim[ivar] + 1];
-
- //values for bin lower bounds
- for(Int_t i=0; i<=thnDim[0]; i++) binEdges[0][i]= -0.9 + (2*0.9)/thnDim[0]*(Double_t)i;
- for(Int_t i=0; i<=thnDim[1]; i++) binEdges[1][i]= -1.2 + (2*1.2)/thnDim[1]*(Double_t)i;
-
- THnSparseS *lfdv = new THnSparseS((const Char_t *)name,(const Char_t *) name,2,thnDim);
+ return GetLinearFitterHistoForce(detector);
- for (int k=0; k<2; k++) {
- lfdv->SetBinEdges(k,binEdges[k]);
- }
- lfdv->Sumw2();
-
- //TH2F *lfdv = new TH2F((const Char_t *)name,(const Char_t *) name
- // ,18,-0.9,0.9,24
- // ,-1.2,1.2);
- //lfdv->SetXTitle("tan(phi_{track})");
- //lfdv->SetYTitle("dy/dt");
- //lfdv->SetZTitle("Number of clusters");
- //lfdv->SetStats(0);
- //lfdv->SetDirectory(0);
-
- fLinearFitterHistoArray.AddAt(lfdv,detector);
- return lfdv;
+}
+//______________________________________________________________________________________
+TH2S* AliTRDCalibraVdriftLinearFit::GetLinearFitterHistoForce(Int_t detector)
+{
+ //
+ // return pointer to TH2F histo
+ // if NULL create a new histo if it doesn't exist allready
+ //
+ if (fLinearFitterHistoArray.UncheckedAt(detector))
+ return (TH2S*)fLinearFitterHistoArray.UncheckedAt(detector);
+
+ // if we are forced and TLinearFitter doesn't yes exist create it
+
+ // new TH2F
+ TString name("LFDV");
+ name += detector;
+ name += "version";
+ name += fVersion;
+
+ TH2S *lfdv = new TH2S((const Char_t *)name,(const Char_t *) name
+ ,36,-0.9,0.9,48
+ ,-1.2,1.2);
+ lfdv->SetXTitle("tan(phi_{track})");
+ lfdv->SetYTitle("dy/dt");
+ lfdv->SetZTitle("Number of clusters");
+ lfdv->SetStats(0);
+ lfdv->SetDirectory(0);
+
+ fLinearFitterHistoArray.AddAt(lfdv,detector);
+ return lfdv;
}
//______________________________________________________________________________________
Bool_t AliTRDCalibraVdriftLinearFit::GetParam(Int_t detector, TVectorD *param)
//
// Fill the 2D histos for debugging
//
- Double_t entries[2] = {tnp,pars1};
- THnSparseS *h = ((THnSparseS *) GetLinearFitterHisto(detector,kTRUE));
+ TH2S *h = ((TH2S *) GetLinearFitterHisto(detector,kTRUE));
Double_t nbentries = h->GetEntries();
- if(nbentries < 5*32767) h->Fill(&entries[0]);
+ if(nbentries < 5*32767) h->Fill(tnp,pars1);
}
//____________Functions fit Online CH2d________________________________________
// Loop over histos
for(Int_t cb = 0; cb < 540; cb++){
- const THnSparseS *linearfitterh = (THnSparseS*)fLinearFitterHistoArray.UncheckedAt(cb);
+ const TH2S *linearfitterhisto = (TH2S*)fLinearFitterHistoArray.UncheckedAt(cb);
//printf("Processing the detector cb %d we find %d\n",cb, (Bool_t) linearfitterhisto);
- if ( linearfitterh != 0 ){
-
- TH2D *linearfitterhisto = linearfitterh->Projection(1,0);
-
+ if ( linearfitterhisto != 0 ){
// Fill a linearfitter
TAxis *xaxis = linearfitterhisto->GetXaxis();
TAxis *yaxis = linearfitterhisto->GetYaxis();
TLinearFitter linearfitter = TLinearFitter(2,"pol1");
//printf("test\n");
+ Double_t integral = linearfitterhisto->Integral();
+ //printf("Integral is %f\n",integral);
+ Bool_t securitybreaking = kFALSE;
+ if(TMath::Abs(integral-1199) < 0.00001) securitybreaking = kTRUE;
for(Int_t ibinx = 0; ibinx < linearfitterhisto->GetNbinsX(); ibinx++){
for(Int_t ibiny = 0; ibiny < linearfitterhisto->GetNbinsY(); ibiny++){
if(linearfitterhisto->GetBinContent(ibinx+1,ibiny+1)>0){
Double_t x = xaxis->GetBinCenter(ibinx+1);
Double_t y = yaxis->GetBinCenter(ibiny+1);
+
for(Int_t k = 0; k < (Int_t)linearfitterhisto->GetBinContent(ibinx+1,ibiny+1); k++){
- linearfitter.AddPoint(&x,y);
- arrayI[cb]++;
+ if(!securitybreaking){
+ linearfitter.AddPoint(&x,y);
+ arrayI[cb]++;
+ }
+ else {
+ if(arrayI[cb]< 1198){
+ linearfitter.AddPoint(&x,y);
+ arrayI[cb]++;
+ }
+ }
}
+
}
}
}
TVectorD *par = new TVectorD(2);
TVectorD pare = TVectorD(2);
TVectorD *parE = new TVectorD(3);
- linearfitter.Eval();
- linearfitter.GetParameters(*par);
- linearfitter.GetErrors(pare);
- Float_t ppointError = TMath::Sqrt(TMath::Abs(linearfitter.GetChisquare())/arrayI[cb]);
- (*parE)[0] = pare[0]*ppointError;
- (*parE)[1] = pare[1]*ppointError;
- (*parE)[2] = (Double_t) arrayI[cb];
- fLinearFitterPArray.AddAt(par,cb);
- fLinearFitterEArray.AddAt(parE,cb);
- //par->Print();
- //parE->Print();
+ //printf("Fit\n");
+ if(((fRobustFit) && (linearfitter.EvalRobust(0.8)==0)) || ((!fRobustFit) && (linearfitter.Eval()==0))) {
+ //if((linearfitter.Eval()==0)) {
+ //printf("Take the param\n");
+ linearfitter.GetParameters(*par);
+ //printf("Done\n");
+ //linearfitter.GetErrors(pare);
+ //Float_t ppointError = TMath::Sqrt(TMath::Abs(linearfitter.GetChisquare())/arrayI[cb]);
+ //(*parE)[0] = pare[0]*ppointError;
+ //(*parE)[1] = pare[1]*ppointError;
+
+ (*parE)[0] = 0.0;
+ (*parE)[1] = 0.0;
+ (*parE)[2] = (Double_t) arrayI[cb];
+ fLinearFitterPArray.AddAt(par,cb);
+ fLinearFitterEArray.AddAt(parE,cb);
+
+ //par->Print();
+ //parE->Print();
+ }
+ //printf("Finish\n");
}
-
- delete linearfitterhisto;
+
+ //delete linearfitterhisto;
}// if something
+
}
+
+ delete [] arrayI;
}
+//____________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;
+}
git://git.uio.no / u / mrichter / AliRoot.git / blobdiff