1 /**************************************************************************
2 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
4 * Author: The ALICE Off-line Project. *
5 * Contributors are mentioned in the code where appropriate. *
7 * Permission to use, copy, modify and distribute this software and its *
8 * documentation strictly for non-commercial purposes is hereby granted *
9 * without fee, provided that the above copyright notice appears in all *
10 * copies and that both the copyright notice and this permission notice *
11 * appear in the supporting documentation. The authors make no claims *
12 * about the suitability of this software for any purpose. It is *
13 * provided "as is" without express or implied warranty. *
14 **************************************************************************/
18 ////////////////////////////////////////////////////////////////////////////
20 // AliTRDCalibraVdriftLinearFit //
22 // Does the Vdrift an ExB calibration by applying a linear fit //
25 // R. Bailhache (R.Bailhache@gsi.de) //
27 ////////////////////////////////////////////////////////////////////////////
30 #include <TObjArray.h>
35 #include <TLinearFitter.h>
39 #include <TTreeStream.h>
40 #include <TGraphErrors.h>
41 #include <TDirectory.h>
42 #include <TTreeStream.h>
46 #include "AliTRDCalibraVdriftLinearFit.h"
48 ClassImp(AliTRDCalibraVdriftLinearFit) /*FOLD00*/
50 //_____________________________________________________________________
51 AliTRDCalibraVdriftLinearFit::AliTRDCalibraVdriftLinearFit() : /*FOLD00*/
55 fLinearFitterHistoArray(540),
56 fLinearFitterPArray(540),
57 fLinearFitterEArray(540),
69 // default constructor
72 //_____________________________________________________________________
73 AliTRDCalibraVdriftLinearFit::AliTRDCalibraVdriftLinearFit(const AliTRDCalibraVdriftLinearFit &ped) : /*FOLD00*/
75 fVersion(ped.fVersion),
76 fNameCalibUsed(ped.fNameCalibUsed),
77 fLinearFitterHistoArray(540),
78 fLinearFitterPArray(540),
79 fLinearFitterEArray(540),
93 for (Int_t idet = 0; idet < 540; idet++){
95 const TVectorD *vectorE = (TVectorD*)ped.fLinearFitterEArray.UncheckedAt(idet);
96 const TVectorD *vectorP = (TVectorD*)ped.fLinearFitterPArray.UncheckedAt(idet);
97 const TH2S *hped = (TH2S*)ped.fLinearFitterHistoArray.UncheckedAt(idet);
99 if ( vectorE != 0x0 ) fLinearFitterEArray.AddAt(new TVectorD(*vectorE), idet);
100 if ( vectorP != 0x0 ) fLinearFitterPArray.AddAt(new TVectorD(*vectorP), idet);
102 TH2S *hNew = (TH2S *)hped->Clone();
103 //hNew->SetDirectory(0);
104 fLinearFitterHistoArray.AddAt(hNew,idet);
108 //_____________________________________________________________________
109 AliTRDCalibraVdriftLinearFit::AliTRDCalibraVdriftLinearFit(const TObjArray &obja) : /*FOLD00*/
113 fLinearFitterHistoArray(540),
114 fLinearFitterPArray(540),
115 fLinearFitterEArray(540),
127 // constructor from a TObjArray
129 for (Int_t idet = 0; idet < 540; idet++){
130 const TH2S *hped = (TH2S*)obja.UncheckedAt(idet);
132 TH2S *hNew = (TH2S *)hped->Clone();
133 //hNew->SetDirectory(0);
134 fLinearFitterHistoArray.AddAt(hNew,idet);
138 //_____________________________________________________________________
139 AliTRDCalibraVdriftLinearFit& AliTRDCalibraVdriftLinearFit::operator = (const AliTRDCalibraVdriftLinearFit &source)
142 // assignment operator
144 if (&source == this) return *this;
145 new (this) AliTRDCalibraVdriftLinearFit(source);
149 //_____________________________________________________________________
150 AliTRDCalibraVdriftLinearFit::~AliTRDCalibraVdriftLinearFit() /*FOLD00*/
155 fLinearFitterHistoArray.SetOwner();
156 fLinearFitterPArray.SetOwner();
157 fLinearFitterEArray.SetOwner();
159 fLinearFitterHistoArray.Delete();
160 fLinearFitterPArray.Delete();
161 fLinearFitterEArray.Delete();
163 if ( fDebugStreamer ) delete fDebugStreamer;
166 //_____________________________________________________________________________
167 void AliTRDCalibraVdriftLinearFit::Copy(TObject &c) const
173 AliTRDCalibraVdriftLinearFit& target = (AliTRDCalibraVdriftLinearFit &) c;
175 // Copy only the histos
176 for (Int_t idet = 0; idet < 540; idet++){
177 if(fLinearFitterHistoArray.UncheckedAt(idet)){
178 TH2S *hped1 = (TH2S *)target.GetLinearFitterHisto(idet,kTRUE);
179 //hped1->SetDirectory(0);
180 hped1->Add((const TH2S *)fLinearFitterHistoArray.UncheckedAt(idet));
187 //_____________________________________________________________________________
188 Long64_t AliTRDCalibraVdriftLinearFit::Merge(const TCollection* list)
190 // Merge list of objects (needed by PROOF)
198 TIterator* iter = list->MakeIterator();
201 // collection of generated histograms
203 while((obj = iter->Next()) != 0)
205 AliTRDCalibraVdriftLinearFit* entry = dynamic_cast<AliTRDCalibraVdriftLinearFit*>(obj);
206 if (entry == 0) continue;
208 // Copy only the histos
209 for (Int_t idet = 0; idet < 540; idet++){
210 if(entry->GetLinearFitterHisto(idet)){
211 TH2S *hped1 = (TH2S *)GetLinearFitterHisto(idet,kTRUE);
212 Double_t entriesa = hped1->GetEntries();
213 Double_t entriesb = ((TH2S *)entry->GetLinearFitterHisto(idet))->GetEntries();
214 if((entriesa + entriesb) < 5*32767) hped1->Add(entry->GetLinearFitterHisto(idet));
223 //_____________________________________________________________________
224 void AliTRDCalibraVdriftLinearFit::Add(const AliTRDCalibraVdriftLinearFit *ped)
232 for (Int_t idet = 0; idet < 540; idet++){
233 const TH2S *hped = (TH2S*)ped->GetLinearFitterHistoNoForce(idet);
234 //printf("idet %d\n",idet);
237 TH2S *hped1 = (TH2S *)GetLinearFitterHisto(idet,kTRUE);
238 Double_t entriesa = hped1->GetEntries();
239 Double_t entriesb = hped->GetEntries();
240 if((entriesa + entriesb) < 5*32767) hped1->Add(hped);
244 //______________________________________________________________________________________
245 TH2S* AliTRDCalibraVdriftLinearFit::AddAll()
248 // return pointer to TH2S of all added histos
252 for(Int_t k=0; k < 540; k++) {
253 TH2S * u = GetLinearFitterHistoForce(k);
255 histo = new TH2S(*u);
256 histo->SetName("sum");
264 //______________________________________________________________________________________
265 TH2S* AliTRDCalibraVdriftLinearFit::GetLinearFitterHisto(Int_t detector, Bool_t force)
268 // return pointer to TH2F histo
269 // if force is true create a new histo if it doesn't exist allready
271 if ( !force || fLinearFitterHistoArray.UncheckedAt(detector) )
272 return (TH2S*)fLinearFitterHistoArray.UncheckedAt(detector);
274 return GetLinearFitterHistoForce(detector);
277 //______________________________________________________________________________________
278 TH2S* AliTRDCalibraVdriftLinearFit::GetLinearFitterHistoForce(Int_t detector)
281 // return pointer to TH2F histo
282 // if NULL create a new histo if it doesn't exist allready
284 if (fLinearFitterHistoArray.UncheckedAt(detector))
285 return (TH2S*)fLinearFitterHistoArray.UncheckedAt(detector);
287 // if we are forced and TLinearFitter doesn't yes exist create it
290 TString name("LFDV");
295 TH2S *lfdv = new TH2S((const Char_t *)name,(const Char_t *) name
296 ,(Int_t)fNbBindx,-fRangedx,fRangedx,(Int_t)fNbBindy
297 ,-fRangedy,fRangedy);
298 lfdv->SetXTitle("tan(phi_{track})");
299 lfdv->SetYTitle("dy/dt");
300 lfdv->SetZTitle("Number of clusters");
302 lfdv->SetDirectory(0);
304 fLinearFitterHistoArray.AddAt(lfdv,detector);
307 //______________________________________________________________________________________
308 Bool_t AliTRDCalibraVdriftLinearFit::GetParam(Int_t detector, TVectorD *param)
311 // return param for this detector
313 if ( fLinearFitterPArray.UncheckedAt(detector) ){
314 const TVectorD *vectorP = (TVectorD*)fLinearFitterPArray.UncheckedAt(detector);
315 if(!param) param = new TVectorD(2);
316 for(Int_t k = 0; k < 2; k++){
317 (*param)[k] = (*vectorP)[k];
324 //______________________________________________________________________________________
325 Bool_t AliTRDCalibraVdriftLinearFit::GetError(Int_t detector, TVectorD *error)
328 // return error for this detector
330 if ( fLinearFitterEArray.UncheckedAt(detector) ){
331 const TVectorD *vectorE = (TVectorD*)fLinearFitterEArray.UncheckedAt(detector);
332 if(!error) error = new TVectorD(3);
333 for(Int_t k = 0; k < 3; k++){
334 (*error)[k] = (*vectorE)[k];
341 //______________________________________________________________________________________
342 void AliTRDCalibraVdriftLinearFit::Update(Int_t detector, Float_t tnp, Float_t pars1)
345 // Fill the 2D histos for debugging
348 TH2S *h = ((TH2S *) GetLinearFitterHisto(detector,kTRUE));
349 Double_t nbentries = h->GetEntries();
350 if(nbentries < 5*32767) h->Fill(tnp,pars1);
353 //____________Functions fit Online CH2d________________________________________
354 void AliTRDCalibraVdriftLinearFit::FillPEArray()
357 // Fill fLinearFitterPArray and fLinearFitterEArray from inside
361 Int_t *arrayI = new Int_t[540];
362 for(Int_t k = 0; k< 540; k++){
367 for(Int_t cb = 0; cb < 540; cb++){
368 const TH2S *linearfitterhisto = (TH2S*)fLinearFitterHistoArray.UncheckedAt(cb);
369 //printf("Processing the detector cb %d we find %d\n",cb, (Bool_t) linearfitterhisto);
371 if ( linearfitterhisto != 0 ){
373 // Fill a linearfitter
374 TAxis *xaxis = linearfitterhisto->GetXaxis();
375 TAxis *yaxis = linearfitterhisto->GetYaxis();
376 TLinearFitter linearfitter = TLinearFitter(2,"pol1");
378 Double_t integral = linearfitterhisto->Integral();
379 //printf("Integral is %f\n",integral);
380 Bool_t securitybreaking = kFALSE;
381 if(TMath::Abs(integral-1199) < 0.00001) securitybreaking = kTRUE;
382 for(Int_t ibinx = 0; ibinx < linearfitterhisto->GetNbinsX(); ibinx++){
383 for(Int_t ibiny = 0; ibiny < linearfitterhisto->GetNbinsY(); ibiny++){
384 if(linearfitterhisto->GetBinContent(ibinx+1,ibiny+1)>0){
385 Double_t x = xaxis->GetBinCenter(ibinx+1);
386 Double_t y = yaxis->GetBinCenter(ibiny+1);
388 for(Int_t k = 0; k < (Int_t)linearfitterhisto->GetBinContent(ibinx+1,ibiny+1); k++){
389 if(!securitybreaking){
390 linearfitter.AddPoint(&x,y);
394 if(arrayI[cb]< 1198){
395 linearfitter.AddPoint(&x,y);
405 //printf("Find %d entries for the detector %d\n",arrayI[cb],cb);
407 // Eval the linear fitter
409 TVectorD *par = new TVectorD(2);
410 TVectorD pare = TVectorD(2);
411 TVectorD *parE = new TVectorD(3);
413 if(((fRobustFit) && (linearfitter.EvalRobust(0.8)==0)) || ((!fRobustFit) && (linearfitter.Eval()==0))) {
414 //if((linearfitter.Eval()==0)) {
415 //printf("Take the param\n");
416 linearfitter.GetParameters(*par);
418 //linearfitter.GetErrors(pare);
419 //Float_t ppointError = TMath::Sqrt(TMath::Abs(linearfitter.GetChisquare())/arrayI[cb]);
420 //(*parE)[0] = pare[0]*ppointError;
421 //(*parE)[1] = pare[1]*ppointError;
425 (*parE)[2] = (Double_t) arrayI[cb];
426 fLinearFitterPArray.AddAt(par,cb);
427 fLinearFitterEArray.AddAt(parE,cb);
432 //printf("Finish\n");
435 //delete linearfitterhisto;
445 //____________Functions fit Online CH2d________________________________________
446 void AliTRDCalibraVdriftLinearFit::FillPEArray2()
449 // Fill fFitterPArray and fFitterEArray from inside
453 TF1 *f1 = new TF1("f1","[0]+[1]*x",-1,1);
455 for(Int_t cb = 0; cb < 540; cb++){
456 //const TH2S *fitterhisto = (TH2S*)fLinearFitterHistoArray.UncheckedAt(cb);
457 TH2S *fitterhisto = (TH2S*)fLinearFitterHistoArray.UncheckedAt(cb);
458 //printf("Processing the detector cb %d we find %d\n",cb, (Bool_t) fitterhisto);
460 if ( fitterhisto != 0 ){
463 TGraphErrors *gg=DrawMS(fitterhisto,nEntries);
465 // Number of points of the TGraphErrors
466 //printf("Number of points %d for detector %d\n",gg->GetN(),cb);
467 if(gg->GetN() < fMinNpointsFit) {
471 //printf("det %d, number of points %d, nEntries %d\n",cb,gg->GetN(),nEntries);
474 TVectorD *par = new TVectorD(2);
475 TVectorD *parE = new TVectorD(3);
478 (*parE)[2] = (Double_t) nEntries;
479 (*par)[0] = f1->GetParameter(0);
480 (*par)[1] = f1->GetParameter(1);
481 fLinearFitterPArray.AddAt(par,cb);
482 fLinearFitterEArray.AddAt(parE,cb);
483 //printf("Value %f for detector %d\n",(*par)[1],cb);
489 if(cb==fSeeDetector) {
490 gStyle->SetPalette(1);
491 gStyle->SetOptStat(1111);
492 gStyle->SetPadBorderMode(0);
493 gStyle->SetCanvasColor(10);
494 gStyle->SetPadLeftMargin(0.13);
495 gStyle->SetPadRightMargin(0.01);
496 TCanvas *stat = new TCanvas("stat","",50,50,600,800);
498 fitterhisto->Draw("colz");
506 if(fDebugLevel==0) delete f1;
509 //_________Helper function__________________________________________________
510 TGraphErrors* AliTRDCalibraVdriftLinearFit::DrawMS(const TH2 *const h2, Int_t &nEntries)
516 TF1 fg("fg", "gaus", -10., 30.);
517 TGraphErrors *gp = new TGraphErrors();
519 TAxis *ax(h2->GetXaxis());
520 TAxis *ay(h2->GetYaxis());
522 for(Int_t ipt(0), jpt(1), ig(0); ipt<ax->GetNbins(); ipt++, jpt++){
523 h1 = h2->ProjectionY("py", jpt, jpt);
524 fg.SetParameter(1, h1->GetMean());
525 //Float_t x(ax->GetBinCenter(jpt));
526 Int_t n=Int_t(h1->Integral(1, h1->GetNbinsX()));
529 //Warning("drawMS()", Form("reject x[%d]=%f on n=%d", jpt, x, n));
532 h1->Fit(&fg, "WWQ0");
534 //Warning("drawMS()", Form("reject x[%d]=%f on NDF=%d", jpt, x, fg.GetNDF()));
537 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;
538 gp->SetPoint(ig, ax->GetBinCenter(jpt), fg.GetParameter(1));
539 gp->SetPointError(ig, 0, TMath::Sqrt(pow(fg.GetParError(1),2) + (1/pow(fg.GetParameter(0),2))));