Print centrality estimator when providing centrality info (Diego)
[u/mrichter/AliRoot.git] / TRD / AliTRDCalibraVdriftLinearFit.cxx
CommitLineData
3a0f6479 1/**************************************************************************
2 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
3 * *
4 * Author: The ALICE Off-line Project. *
5 * Contributors are mentioned in the code where appropriate. *
6 * *
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 **************************************************************************/
15
16/* $Id$ */
17
0bc7827a 18////////////////////////////////////////////////////////////////////////////
19// //
20// AliTRDCalibraVdriftLinearFit //
21// //
22// Does the Vdrift an ExB calibration by applying a linear fit //
23// //
24// Author: //
25// R. Bailhache (R.Bailhache@gsi.de) //
26// //
27////////////////////////////////////////////////////////////////////////////
28
3a0f6479 29//Root includes
30#include <TObjArray.h>
31#include <TH2F.h>
32#include <TString.h>
33#include <TVectorD.h>
34#include <TAxis.h>
35#include <TLinearFitter.h>
36#include <TMath.h>
661c7be6 37#include <TStyle.h>
38#include <TCanvas.h>
3a0f6479 39#include <TTreeStream.h>
661c7be6 40#include <TGraphErrors.h>
41#include <TDirectory.h>
42#include <TTreeStream.h>
43#include <TF1.h>
3a0f6479 44
45//header file
46#include "AliTRDCalibraVdriftLinearFit.h"
47
48ClassImp(AliTRDCalibraVdriftLinearFit) /*FOLD00*/
49
50//_____________________________________________________________________
51AliTRDCalibraVdriftLinearFit::AliTRDCalibraVdriftLinearFit() : /*FOLD00*/
52 TObject(),
53 fVersion(0),
b2277aa2 54 fNameCalibUsed(""),
3a0f6479 55 fLinearFitterHistoArray(540),
56 fLinearFitterPArray(540),
2a1a7b36 57 fLinearFitterEArray(540),
661c7be6 58 fRobustFit(kTRUE),
532be2d4 59 fMinNpointsFit(11),
4f3bd513 60 fNbBindx(32),
61 fNbBindy(70),
62 fRangedx(0.8),
63 fRangedy(1.4),
661c7be6 64 fDebugStreamer(0x0),
65 fDebugLevel(0),
66 fSeeDetector(0)
3a0f6479 67{
ba94744a 68 //
69 // default constructor
70 //
3a0f6479 71}
3a0f6479 72//_____________________________________________________________________
73AliTRDCalibraVdriftLinearFit::AliTRDCalibraVdriftLinearFit(const AliTRDCalibraVdriftLinearFit &ped) : /*FOLD00*/
74 TObject(ped),
75 fVersion(ped.fVersion),
b2277aa2 76 fNameCalibUsed(ped.fNameCalibUsed),
3a0f6479 77 fLinearFitterHistoArray(540),
78 fLinearFitterPArray(540),
2a1a7b36 79 fLinearFitterEArray(540),
661c7be6 80 fRobustFit(kTRUE),
532be2d4 81 fMinNpointsFit(10),
4f3bd513 82 fNbBindx(32),
83 fNbBindy(70),
84 fRangedx(0.8),
85 fRangedy(1.4),
661c7be6 86 fDebugStreamer(0x0),
87 fDebugLevel(0),
88 fSeeDetector(0)
3a0f6479 89{
90 //
91 // copy constructor
92 //
93 for (Int_t idet = 0; idet < 540; idet++){
ba94744a 94
3a0f6479 95 const TVectorD *vectorE = (TVectorD*)ped.fLinearFitterEArray.UncheckedAt(idet);
96 const TVectorD *vectorP = (TVectorD*)ped.fLinearFitterPArray.UncheckedAt(idet);
0ca1720e 97 const TH2S *hped = (TH2S*)ped.fLinearFitterHistoArray.UncheckedAt(idet);
3a0f6479 98
99 if ( vectorE != 0x0 ) fLinearFitterEArray.AddAt(new TVectorD(*vectorE), idet);
100 if ( vectorP != 0x0 ) fLinearFitterPArray.AddAt(new TVectorD(*vectorP), idet);
101 if ( hped != 0x0 ){
0ca1720e 102 TH2S *hNew = (TH2S *)hped->Clone();
ba94744a 103 //hNew->SetDirectory(0);
3a0f6479 104 fLinearFitterHistoArray.AddAt(hNew,idet);
105 }
3a0f6479 106 }
107}
108//_____________________________________________________________________
d0569428 109AliTRDCalibraVdriftLinearFit::AliTRDCalibraVdriftLinearFit(const TObjArray &obja) : /*FOLD00*/
110 TObject(),
111 fVersion(0),
b2277aa2 112 fNameCalibUsed(""),
d0569428 113 fLinearFitterHistoArray(540),
114 fLinearFitterPArray(540),
2a1a7b36 115 fLinearFitterEArray(540),
661c7be6 116 fRobustFit(kTRUE),
532be2d4 117 fMinNpointsFit(10),
4f3bd513 118 fNbBindx(32),
119 fNbBindy(70),
120 fRangedx(0.8),
121 fRangedy(1.4),
661c7be6 122 fDebugStreamer(0x0),
123 fDebugLevel(0),
124 fSeeDetector(0)
d0569428 125{
126 //
127 // constructor from a TObjArray
128 //
129 for (Int_t idet = 0; idet < 540; idet++){
0ca1720e 130 const TH2S *hped = (TH2S*)obja.UncheckedAt(idet);
d0569428 131 if ( hped != 0x0 ){
0ca1720e 132 TH2S *hNew = (TH2S *)hped->Clone();
ba94744a 133 //hNew->SetDirectory(0);
d0569428 134 fLinearFitterHistoArray.AddAt(hNew,idet);
135 }
136 }
137}
138//_____________________________________________________________________
3a0f6479 139AliTRDCalibraVdriftLinearFit& AliTRDCalibraVdriftLinearFit::operator = (const AliTRDCalibraVdriftLinearFit &source)
140{
141 //
142 // assignment operator
143 //
144 if (&source == this) return *this;
145 new (this) AliTRDCalibraVdriftLinearFit(source);
146
147 return *this;
148}
149//_____________________________________________________________________
150AliTRDCalibraVdriftLinearFit::~AliTRDCalibraVdriftLinearFit() /*FOLD00*/
151{
152 //
153 // destructor
154 //
f5cda287 155 fLinearFitterHistoArray.SetOwner();
156 fLinearFitterPArray.SetOwner();
157 fLinearFitterEArray.SetOwner();
158
1ca79a00 159 fLinearFitterHistoArray.Delete();
160 fLinearFitterPArray.Delete();
161 fLinearFitterEArray.Delete();
162
661c7be6 163 if ( fDebugStreamer ) delete fDebugStreamer;
164
3a0f6479 165}
64942b85 166//_____________________________________________________________________________
167void AliTRDCalibraVdriftLinearFit::Copy(TObject &c) const
168{
169 //
170 // Copy function
171 //
172
173 AliTRDCalibraVdriftLinearFit& target = (AliTRDCalibraVdriftLinearFit &) c;
174
175 // Copy only the histos
176 for (Int_t idet = 0; idet < 540; idet++){
177 if(fLinearFitterHistoArray.UncheckedAt(idet)){
0ca1720e 178 TH2S *hped1 = (TH2S *)target.GetLinearFitterHisto(idet,kTRUE);
ba94744a 179 //hped1->SetDirectory(0);
0ca1720e 180 hped1->Add((const TH2S *)fLinearFitterHistoArray.UncheckedAt(idet));
64942b85 181 }
64942b85 182 }
ba94744a 183
64942b85 184 TObject::Copy(c);
185
186}
187//_____________________________________________________________________________
6bbdc11a 188Long64_t AliTRDCalibraVdriftLinearFit::Merge(const TCollection* list)
64942b85 189{
190 // Merge list of objects (needed by PROOF)
191
192 if (!list)
193 return 0;
194
195 if (list->IsEmpty())
196 return 1;
197
198 TIterator* iter = list->MakeIterator();
199 TObject* obj = 0;
200
201 // collection of generated histograms
202 Int_t count=0;
203 while((obj = iter->Next()) != 0)
ba94744a 204 {
205 AliTRDCalibraVdriftLinearFit* entry = dynamic_cast<AliTRDCalibraVdriftLinearFit*>(obj);
206 if (entry == 0) continue;
207
208 // Copy only the histos
209 for (Int_t idet = 0; idet < 540; idet++){
210 if(entry->GetLinearFitterHisto(idet)){
0ca1720e 211 TH2S *hped1 = (TH2S *)GetLinearFitterHisto(idet,kTRUE);
ba94744a 212 Double_t entriesa = hped1->GetEntries();
0ca1720e 213 Double_t entriesb = ((TH2S *)entry->GetLinearFitterHisto(idet))->GetEntries();
ba94744a 214 if((entriesa + entriesb) < 5*32767) hped1->Add(entry->GetLinearFitterHisto(idet));
215 }
64942b85 216 }
217
ba94744a 218 count++;
64942b85 219 }
64942b85 220
221 return count;
222}
3a0f6479 223//_____________________________________________________________________
fe4ee353 224void AliTRDCalibraVdriftLinearFit::Add(const AliTRDCalibraVdriftLinearFit *ped)
3a0f6479 225{
226 //
227 // Add histo
228 //
229
230 fVersion++;
231
232 for (Int_t idet = 0; idet < 540; idet++){
fe4ee353 233 const TH2S *hped = (TH2S*)ped->GetLinearFitterHistoNoForce(idet);
3a0f6479 234 //printf("idet %d\n",idet);
235 if ( hped != 0x0 ){
236 //printf("add\n");
0ca1720e 237 TH2S *hped1 = (TH2S *)GetLinearFitterHisto(idet,kTRUE);
ba94744a 238 Double_t entriesa = hped1->GetEntries();
239 Double_t entriesb = hped->GetEntries();
240 if((entriesa + entriesb) < 5*32767) hped1->Add(hped);
3a0f6479 241 }
242 }
243}
244//______________________________________________________________________________________
4f3bd513 245TH2S* AliTRDCalibraVdriftLinearFit::AddAll()
246{
247 //
248 // return pointer to TH2S of all added histos
249 //
250
251 TH2S *histo = 0x0;
252 for(Int_t k=0; k < 540; k++) {
253 TH2S * u = GetLinearFitterHistoForce(k);
254 if(k == 0) {
255 histo = new TH2S(*u);
256 histo->SetName("sum");
257 }
258 else histo->Add(u);
259 }
260
261 return histo;
262
263}
264//______________________________________________________________________________________
0ca1720e 265TH2S* AliTRDCalibraVdriftLinearFit::GetLinearFitterHisto(Int_t detector, Bool_t force)
3a0f6479 266{
267 //
268 // return pointer to TH2F histo
269 // if force is true create a new histo if it doesn't exist allready
270 //
271 if ( !force || fLinearFitterHistoArray.UncheckedAt(detector) )
0ca1720e 272 return (TH2S*)fLinearFitterHistoArray.UncheckedAt(detector);
3a0f6479 273
fe4ee353 274 return GetLinearFitterHistoForce(detector);
275
276}
277//______________________________________________________________________________________
278TH2S* AliTRDCalibraVdriftLinearFit::GetLinearFitterHistoForce(Int_t detector)
279{
280 //
281 // return pointer to TH2F histo
282 // if NULL create a new histo if it doesn't exist allready
283 //
284 if (fLinearFitterHistoArray.UncheckedAt(detector))
285 return (TH2S*)fLinearFitterHistoArray.UncheckedAt(detector);
286
287 // if we are forced and TLinearFitter doesn't yes exist create it
288
289 // new TH2F
290 TString name("LFDV");
291 name += detector;
292 name += "version";
293 name += fVersion;
294
295 TH2S *lfdv = new TH2S((const Char_t *)name,(const Char_t *) name
4f3bd513 296 ,(Int_t)fNbBindx,-fRangedx,fRangedx,(Int_t)fNbBindy
297 ,-fRangedy,fRangedy);
fe4ee353 298 lfdv->SetXTitle("tan(phi_{track})");
299 lfdv->SetYTitle("dy/dt");
300 lfdv->SetZTitle("Number of clusters");
301 lfdv->SetStats(0);
302 lfdv->SetDirectory(0);
303
304 fLinearFitterHistoArray.AddAt(lfdv,detector);
305 return lfdv;
3a0f6479 306}
307//______________________________________________________________________________________
308Bool_t AliTRDCalibraVdriftLinearFit::GetParam(Int_t detector, TVectorD *param)
309{
310 //
311 // return param for this detector
312 //
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];
318 }
319 return kTRUE;
320 }
321 else return kFALSE;
322
323}
324//______________________________________________________________________________________
325Bool_t AliTRDCalibraVdriftLinearFit::GetError(Int_t detector, TVectorD *error)
326{
327 //
328 // return error for this detector
329 //
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];
335 }
336 return kTRUE;
337 }
338 else return kFALSE;
339
340}
341//______________________________________________________________________________________
342void AliTRDCalibraVdriftLinearFit::Update(Int_t detector, Float_t tnp, Float_t pars1)
343{
344 //
345 // Fill the 2D histos for debugging
346 //
ba94744a 347
0ca1720e 348 TH2S *h = ((TH2S *) GetLinearFitterHisto(detector,kTRUE));
ba94744a 349 Double_t nbentries = h->GetEntries();
0ca1720e 350 if(nbentries < 5*32767) h->Fill(tnp,pars1);
3a0f6479 351
352}
353//____________Functions fit Online CH2d________________________________________
354void AliTRDCalibraVdriftLinearFit::FillPEArray()
355{
356 //
357 // Fill fLinearFitterPArray and fLinearFitterEArray from inside
358 //
359
360
361 Int_t *arrayI = new Int_t[540];
362 for(Int_t k = 0; k< 540; k++){
363 arrayI[k] = 0;
364 }
365
366 // Loop over histos
367 for(Int_t cb = 0; cb < 540; cb++){
0ca1720e 368 const TH2S *linearfitterhisto = (TH2S*)fLinearFitterHistoArray.UncheckedAt(cb);
3a0f6479 369 //printf("Processing the detector cb %d we find %d\n",cb, (Bool_t) linearfitterhisto);
370
0ca1720e 371 if ( linearfitterhisto != 0 ){
3a0f6479 372
373 // Fill a linearfitter
374 TAxis *xaxis = linearfitterhisto->GetXaxis();
375 TAxis *yaxis = linearfitterhisto->GetYaxis();
376 TLinearFitter linearfitter = TLinearFitter(2,"pol1");
377 //printf("test\n");
a5dcf618 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;
3a0f6479 382 for(Int_t ibinx = 0; ibinx < linearfitterhisto->GetNbinsX(); ibinx++){
383 for(Int_t ibiny = 0; ibiny < linearfitterhisto->GetNbinsY(); ibiny++){
8a606eb5 384 if(linearfitterhisto->GetBinContent(ibinx+1,ibiny+1)>0){
3a0f6479 385 Double_t x = xaxis->GetBinCenter(ibinx+1);
386 Double_t y = yaxis->GetBinCenter(ibiny+1);
a5dcf618 387
3a0f6479 388 for(Int_t k = 0; k < (Int_t)linearfitterhisto->GetBinContent(ibinx+1,ibiny+1); k++){
a5dcf618 389 if(!securitybreaking){
390 linearfitter.AddPoint(&x,y);
391 arrayI[cb]++;
392 }
393 else {
394 if(arrayI[cb]< 1198){
395 linearfitter.AddPoint(&x,y);
396 arrayI[cb]++;
397 }
398 }
3a0f6479 399 }
a5dcf618 400
3a0f6479 401 }
402 }
403 }
404
405 //printf("Find %d entries for the detector %d\n",arrayI[cb],cb);
406
407 // Eval the linear fitter
408 if(arrayI[cb]>10){
409 TVectorD *par = new TVectorD(2);
410 TVectorD pare = TVectorD(2);
411 TVectorD *parE = new TVectorD(3);
a5dcf618 412 //printf("Fit\n");
2a1a7b36 413 if(((fRobustFit) && (linearfitter.EvalRobust(0.8)==0)) || ((!fRobustFit) && (linearfitter.Eval()==0))) {
a5dcf618 414 //if((linearfitter.Eval()==0)) {
415 //printf("Take the param\n");
c1105918 416 linearfitter.GetParameters(*par);
a5dcf618 417 //printf("Done\n");
c1105918 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;
422
423 (*parE)[0] = 0.0;
424 (*parE)[1] = 0.0;
425 (*parE)[2] = (Double_t) arrayI[cb];
426 fLinearFitterPArray.AddAt(par,cb);
427 fLinearFitterEArray.AddAt(parE,cb);
428
429 //par->Print();
430 //parE->Print();
431 }
a5dcf618 432 //printf("Finish\n");
ba94744a 433 }
c1105918 434
07f9b7f2 435 //delete linearfitterhisto;
ba94744a 436
3a0f6479 437 }// if something
fdc15553 438
3a0f6479 439 }
a987273c 440
441 delete [] arrayI;
3a0f6479 442
d1e1a064 443}
c1105918 444
661c7be6 445//____________Functions fit Online CH2d________________________________________
446void AliTRDCalibraVdriftLinearFit::FillPEArray2()
447{
448 //
449 // Fill fFitterPArray and fFitterEArray from inside
450 //
451
452 // Loop over histos
453 TF1 *f1 = new TF1("f1","[0]+[1]*x",-1,1);
454
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);
459
460 if ( fitterhisto != 0 ){
461
462 Int_t nEntries=0;
463 TGraphErrors *gg=DrawMS(fitterhisto,nEntries);
ce42378b 464 if (!gg) continue;
661c7be6 465 // Number of points of the TGraphErrors
532be2d4 466 //printf("Number of points %d for detector %d\n",gg->GetN(),cb);
467 if(gg->GetN() < fMinNpointsFit) {
661c7be6 468 if(gg) delete gg;
469 continue;
470 }
471 //printf("det %d, number of points %d, nEntries %d\n",cb,gg->GetN(),nEntries);
472 gg->Fit(f1,"Q0");
473
474 TVectorD *par = new TVectorD(2);
475 TVectorD *parE = new TVectorD(3);
476 (*parE)[0] = 0.0;
477 (*parE)[1] = 0.0;
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);
532be2d4 483 //printf("Value %f for detector %d\n",(*par)[1],cb);
661c7be6 484
485 if(fDebugLevel==0) {
486 if(gg) delete gg;
487 }
488 else {
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);
497 stat->cd(1);
498 fitterhisto->Draw("colz");
499 gg->Draw("P");
500 f1->Draw("same");
501 break;
502 }
503 }
504 }
505 }
506 if(fDebugLevel==0) delete f1;
507}
508
509//_________Helper function__________________________________________________
510TGraphErrors* AliTRDCalibraVdriftLinearFit::DrawMS(const TH2 *const h2, Int_t &nEntries)
511{
92df0132 512 //
513 // Debug function
514 //
515
661c7be6 516 TF1 fg("fg", "gaus", -10., 30.);
517 TGraphErrors *gp = new TGraphErrors();
c1105918 518
661c7be6 519 TAxis *ax(h2->GetXaxis());
520 TAxis *ay(h2->GetYaxis());
521 TH1D *h1(NULL);
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()));
527 nEntries+=n;
528 if(n < 15){
529 //Warning("drawMS()", Form("reject x[%d]=%f on n=%d", jpt, x, n));
530 continue;
531 }
532 h1->Fit(&fg, "WWQ0");
533 if(fg.GetNDF()<2){
534 //Warning("drawMS()", Form("reject x[%d]=%f on NDF=%d", jpt, x, fg.GetNDF()));
535 continue;
536 }
64c7fdf8 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;
661c7be6 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))));
540 ig++;
541 }
542 delete h1;
543 return gp;
544}