X-Git-Url: http://git.uio.no/git/?a=blobdiff_plain;f=TRD%2FAliTRDCalibraVdriftLinearFit.cxx;h=855fceb5319c238fef97c25260da9814496f2df8;hb=57fa828cce8c40e08fa2d69490209a0ffa386a7b;hp=3ac80ee67fe86c1c597bb16a588cc672e19a5142;hpb=a4d55cb0539cb977bf7c143cac461ea57aea2079;p=u%2Fmrichter%2FAliRoot.git diff --git a/TRD/AliTRDCalibraVdriftLinearFit.cxx b/TRD/AliTRDCalibraVdriftLinearFit.cxx index 3ac80ee67fe..855fceb5319 100644 --- a/TRD/AliTRDCalibraVdriftLinearFit.cxx +++ b/TRD/AliTRDCalibraVdriftLinearFit.cxx @@ -34,7 +34,13 @@ #include #include #include +#include +#include #include +#include +#include +#include +#include //header file #include "AliTRDCalibraVdriftLinearFit.h" @@ -45,56 +51,86 @@ 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 - // + // + // 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), + fNbBindx(32), + fNbBindy(70), + fRangedx(0.8), + fRangedy(1.4), + fDebugStreamer(0x0), + fDebugLevel(0), + fSeeDetector(0) { // // copy constructor // for (Int_t idet = 0; idet < 540; idet++){ + const TVectorD *vectorE = (TVectorD*)ped.fLinearFitterEArray.UncheckedAt(idet); const TVectorD *vectorP = (TVectorD*)ped.fLinearFitterPArray.UncheckedAt(idet); - const TH2F *hped = (TH2F*)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 ){ - TH2F *hNew = new TH2F(*hped); - hNew->SetDirectory(0); + 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), + fNbBindx(32), + fNbBindy(70), + fRangedx(0.8), + fRangedy(1.4), + fDebugStreamer(0x0), + fDebugLevel(0), + fSeeDetector(0) { // // constructor from a TObjArray // for (Int_t idet = 0; idet < 540; idet++){ - const TH2F *hped = (TH2F*)obja.UncheckedAt(idet); + const TH2S *hped = (TH2S*)obja.UncheckedAt(idet); if ( hped != 0x0 ){ - TH2F *hNew = new TH2F(*hped); - hNew->SetDirectory(0); + TH2S *hNew = (TH2S *)hped->Clone(); + //hNew->SetDirectory(0); fLinearFitterHistoArray.AddAt(hNew,idet); } } @@ -116,10 +152,16 @@ AliTRDCalibraVdriftLinearFit::~AliTRDCalibraVdriftLinearFit() /*FOLD00*/ // // destructor // + fLinearFitterHistoArray.SetOwner(); + fLinearFitterPArray.SetOwner(); + fLinearFitterEArray.SetOwner(); + fLinearFitterHistoArray.Delete(); fLinearFitterPArray.Delete(); fLinearFitterEArray.Delete(); + if ( fDebugStreamer ) delete fDebugStreamer; + } //_____________________________________________________________________________ void AliTRDCalibraVdriftLinearFit::Copy(TObject &c) const @@ -133,18 +175,17 @@ void AliTRDCalibraVdriftLinearFit::Copy(TObject &c) const // Copy only the histos for (Int_t idet = 0; idet < 540; idet++){ if(fLinearFitterHistoArray.UncheckedAt(idet)){ - TH2F *hped1 = target.GetLinearFitterHisto(idet,kTRUE); - hped1->SetDirectory(0); - hped1->Add((const TH1F *)fLinearFitterHistoArray.UncheckedAt(idet)); + TH2S *hped1 = (TH2S *)target.GetLinearFitterHisto(idet,kTRUE); + //hped1->SetDirectory(0); + hped1->Add((const TH2S *)fLinearFitterHistoArray.UncheckedAt(idet)); } - } - + TObject::Copy(c); } //_____________________________________________________________________________ -Long64_t AliTRDCalibraVdriftLinearFit::Merge(TCollection* list) +Long64_t AliTRDCalibraVdriftLinearFit::Merge(const TCollection* list) { // Merge list of objects (needed by PROOF) @@ -160,27 +201,27 @@ Long64_t AliTRDCalibraVdriftLinearFit::Merge(TCollection* list) // collection of generated histograms Int_t count=0; while((obj = iter->Next()) != 0) - { - AliTRDCalibraVdriftLinearFit* entry = dynamic_cast(obj); - if (entry == 0) continue; - - // Copy only the histos - for (Int_t idet = 0; idet < 540; idet++){ - if(entry->GetLinearFitterHisto(idet)){ - TH2F *hped1 = GetLinearFitterHisto(idet,kTRUE); - hped1->SetDirectory(0); - hped1->Add(entry->GetLinearFitterHisto(idet)); + { + AliTRDCalibraVdriftLinearFit* entry = dynamic_cast(obj); + if (entry == 0) continue; + + // Copy only the histos + for (Int_t idet = 0; idet < 540; idet++){ + if(entry->GetLinearFitterHisto(idet)){ + TH2S *hped1 = (TH2S *)GetLinearFitterHisto(idet,kTRUE); + Double_t entriesa = hped1->GetEntries(); + Double_t entriesb = ((TH2S *)entry->GetLinearFitterHisto(idet))->GetEntries(); + if((entriesa + entriesb) < 5*32767) hped1->Add(entry->GetLinearFitterHisto(idet)); + } } + count++; } - - count++; - } return count; } //_____________________________________________________________________ -void AliTRDCalibraVdriftLinearFit::Add(AliTRDCalibraVdriftLinearFit *ped) +void AliTRDCalibraVdriftLinearFit::Add(const AliTRDCalibraVdriftLinearFit *ped) { // // Add histo @@ -189,48 +230,79 @@ void AliTRDCalibraVdriftLinearFit::Add(AliTRDCalibraVdriftLinearFit *ped) fVersion++; for (Int_t idet = 0; idet < 540; idet++){ - const TH2F *hped = (TH2F*)ped->GetLinearFitterHisto(idet); + const TH2S *hped = (TH2S*)ped->GetLinearFitterHistoNoForce(idet); //printf("idet %d\n",idet); if ( hped != 0x0 ){ //printf("add\n"); - TH2F *hped1 = GetLinearFitterHisto(idet,kTRUE); - //printf("test2\n"); - hped1->SetDirectory(0); - //printf("test4\n"); - hped1->Add(hped); - //printf("test3\n"); + TH2S *hped1 = (TH2S *)GetLinearFitterHisto(idet,kTRUE); + Double_t entriesa = hped1->GetEntries(); + Double_t entriesb = hped->GetEntries(); + if((entriesa + entriesb) < 5*32767) hped1->Add(hped); + } + } +} +//______________________________________________________________________________________ +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; + } //______________________________________________________________________________________ -TH2F* 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 (TH2F*)fLinearFitterHistoArray.UncheckedAt(detector); + return (TH2S*)fLinearFitterHistoArray.UncheckedAt(detector); - // if we are forced and TLinearFitter doesn't yes exist create it + return GetLinearFitterHistoForce(detector); - // new TH2F - TString name("LFDV"); - name += detector; - name += "version"; - name += fVersion; - - TH2F *lfdv = new TH2F((const Char_t *)name,(const Char_t *) name - ,40,-1.0,1.0,50 - ,-2.0,2.0); - 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 + ,(Int_t)fNbBindx,-fRangedx,fRangedx,(Int_t)fNbBindy + ,-fRangedy,fRangedy); + 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) @@ -272,8 +344,10 @@ void AliTRDCalibraVdriftLinearFit::Update(Int_t detector, Float_t tnp, Float_t p // // Fill the 2D histos for debugging // - - ((TH2F *) GetLinearFitterHisto(detector,kTRUE))->Fill(tnp,pars1); + + TH2S *h = ((TH2S *) GetLinearFitterHisto(detector,kTRUE)); + Double_t nbentries = h->GetEntries(); + if(nbentries < 5*32767) h->Fill(tnp,pars1); } //____________Functions fit Online CH2d________________________________________ @@ -291,7 +365,7 @@ void AliTRDCalibraVdriftLinearFit::FillPEArray() // Loop over histos for(Int_t cb = 0; cb < 540; cb++){ - const TH2F *linearfitterhisto = (TH2F*)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 ( linearfitterhisto != 0 ){ @@ -301,15 +375,29 @@ void AliTRDCalibraVdriftLinearFit::FillPEArray() 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]++; + } + } } + } } } @@ -321,20 +409,136 @@ void AliTRDCalibraVdriftLinearFit::FillPEArray() 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; + }// 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); iptGetNbins(); 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)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; +}