Adjust needed statistics for Vd calibration and normalization with TPC signal for...
[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),
661c7be6 60 fDebugStreamer(0x0),
61 fDebugLevel(0),
62 fSeeDetector(0)
3a0f6479 63{
ba94744a 64 //
65 // default constructor
66 //
3a0f6479 67}
3a0f6479 68//_____________________________________________________________________
69AliTRDCalibraVdriftLinearFit::AliTRDCalibraVdriftLinearFit(const AliTRDCalibraVdriftLinearFit &ped) : /*FOLD00*/
70 TObject(ped),
71 fVersion(ped.fVersion),
b2277aa2 72 fNameCalibUsed(ped.fNameCalibUsed),
3a0f6479 73 fLinearFitterHistoArray(540),
74 fLinearFitterPArray(540),
2a1a7b36 75 fLinearFitterEArray(540),
661c7be6 76 fRobustFit(kTRUE),
532be2d4 77 fMinNpointsFit(10),
661c7be6 78 fDebugStreamer(0x0),
79 fDebugLevel(0),
80 fSeeDetector(0)
3a0f6479 81{
82 //
83 // copy constructor
84 //
85 for (Int_t idet = 0; idet < 540; idet++){
ba94744a 86
3a0f6479 87 const TVectorD *vectorE = (TVectorD*)ped.fLinearFitterEArray.UncheckedAt(idet);
88 const TVectorD *vectorP = (TVectorD*)ped.fLinearFitterPArray.UncheckedAt(idet);
0ca1720e 89 const TH2S *hped = (TH2S*)ped.fLinearFitterHistoArray.UncheckedAt(idet);
3a0f6479 90
91 if ( vectorE != 0x0 ) fLinearFitterEArray.AddAt(new TVectorD(*vectorE), idet);
92 if ( vectorP != 0x0 ) fLinearFitterPArray.AddAt(new TVectorD(*vectorP), idet);
93 if ( hped != 0x0 ){
0ca1720e 94 TH2S *hNew = (TH2S *)hped->Clone();
ba94744a 95 //hNew->SetDirectory(0);
3a0f6479 96 fLinearFitterHistoArray.AddAt(hNew,idet);
97 }
3a0f6479 98 }
99}
100//_____________________________________________________________________
d0569428 101AliTRDCalibraVdriftLinearFit::AliTRDCalibraVdriftLinearFit(const TObjArray &obja) : /*FOLD00*/
102 TObject(),
103 fVersion(0),
b2277aa2 104 fNameCalibUsed(""),
d0569428 105 fLinearFitterHistoArray(540),
106 fLinearFitterPArray(540),
2a1a7b36 107 fLinearFitterEArray(540),
661c7be6 108 fRobustFit(kTRUE),
532be2d4 109 fMinNpointsFit(10),
661c7be6 110 fDebugStreamer(0x0),
111 fDebugLevel(0),
112 fSeeDetector(0)
d0569428 113{
114 //
115 // constructor from a TObjArray
116 //
117 for (Int_t idet = 0; idet < 540; idet++){
0ca1720e 118 const TH2S *hped = (TH2S*)obja.UncheckedAt(idet);
d0569428 119 if ( hped != 0x0 ){
0ca1720e 120 TH2S *hNew = (TH2S *)hped->Clone();
ba94744a 121 //hNew->SetDirectory(0);
d0569428 122 fLinearFitterHistoArray.AddAt(hNew,idet);
123 }
124 }
125}
126//_____________________________________________________________________
3a0f6479 127AliTRDCalibraVdriftLinearFit& AliTRDCalibraVdriftLinearFit::operator = (const AliTRDCalibraVdriftLinearFit &source)
128{
129 //
130 // assignment operator
131 //
132 if (&source == this) return *this;
133 new (this) AliTRDCalibraVdriftLinearFit(source);
134
135 return *this;
136}
137//_____________________________________________________________________
138AliTRDCalibraVdriftLinearFit::~AliTRDCalibraVdriftLinearFit() /*FOLD00*/
139{
140 //
141 // destructor
142 //
f5cda287 143 fLinearFitterHistoArray.SetOwner();
144 fLinearFitterPArray.SetOwner();
145 fLinearFitterEArray.SetOwner();
146
1ca79a00 147 fLinearFitterHistoArray.Delete();
148 fLinearFitterPArray.Delete();
149 fLinearFitterEArray.Delete();
150
661c7be6 151 if ( fDebugStreamer ) delete fDebugStreamer;
152
3a0f6479 153}
64942b85 154//_____________________________________________________________________________
155void AliTRDCalibraVdriftLinearFit::Copy(TObject &c) const
156{
157 //
158 // Copy function
159 //
160
161 AliTRDCalibraVdriftLinearFit& target = (AliTRDCalibraVdriftLinearFit &) c;
162
163 // Copy only the histos
164 for (Int_t idet = 0; idet < 540; idet++){
165 if(fLinearFitterHistoArray.UncheckedAt(idet)){
0ca1720e 166 TH2S *hped1 = (TH2S *)target.GetLinearFitterHisto(idet,kTRUE);
ba94744a 167 //hped1->SetDirectory(0);
0ca1720e 168 hped1->Add((const TH2S *)fLinearFitterHistoArray.UncheckedAt(idet));
64942b85 169 }
64942b85 170 }
ba94744a 171
64942b85 172 TObject::Copy(c);
173
174}
175//_____________________________________________________________________________
6bbdc11a 176Long64_t AliTRDCalibraVdriftLinearFit::Merge(const TCollection* list)
64942b85 177{
178 // Merge list of objects (needed by PROOF)
179
180 if (!list)
181 return 0;
182
183 if (list->IsEmpty())
184 return 1;
185
186 TIterator* iter = list->MakeIterator();
187 TObject* obj = 0;
188
189 // collection of generated histograms
190 Int_t count=0;
191 while((obj = iter->Next()) != 0)
ba94744a 192 {
193 AliTRDCalibraVdriftLinearFit* entry = dynamic_cast<AliTRDCalibraVdriftLinearFit*>(obj);
194 if (entry == 0) continue;
195
196 // Copy only the histos
197 for (Int_t idet = 0; idet < 540; idet++){
198 if(entry->GetLinearFitterHisto(idet)){
0ca1720e 199 TH2S *hped1 = (TH2S *)GetLinearFitterHisto(idet,kTRUE);
ba94744a 200 Double_t entriesa = hped1->GetEntries();
0ca1720e 201 Double_t entriesb = ((TH2S *)entry->GetLinearFitterHisto(idet))->GetEntries();
ba94744a 202 if((entriesa + entriesb) < 5*32767) hped1->Add(entry->GetLinearFitterHisto(idet));
203 }
64942b85 204 }
205
ba94744a 206 count++;
64942b85 207 }
64942b85 208
209 return count;
210}
3a0f6479 211//_____________________________________________________________________
fe4ee353 212void AliTRDCalibraVdriftLinearFit::Add(const AliTRDCalibraVdriftLinearFit *ped)
3a0f6479 213{
214 //
215 // Add histo
216 //
217
218 fVersion++;
219
220 for (Int_t idet = 0; idet < 540; idet++){
fe4ee353 221 const TH2S *hped = (TH2S*)ped->GetLinearFitterHistoNoForce(idet);
3a0f6479 222 //printf("idet %d\n",idet);
223 if ( hped != 0x0 ){
224 //printf("add\n");
0ca1720e 225 TH2S *hped1 = (TH2S *)GetLinearFitterHisto(idet,kTRUE);
ba94744a 226 Double_t entriesa = hped1->GetEntries();
227 Double_t entriesb = hped->GetEntries();
228 if((entriesa + entriesb) < 5*32767) hped1->Add(hped);
3a0f6479 229 }
230 }
231}
232//______________________________________________________________________________________
0ca1720e 233TH2S* AliTRDCalibraVdriftLinearFit::GetLinearFitterHisto(Int_t detector, Bool_t force)
3a0f6479 234{
235 //
236 // return pointer to TH2F histo
237 // if force is true create a new histo if it doesn't exist allready
238 //
239 if ( !force || fLinearFitterHistoArray.UncheckedAt(detector) )
0ca1720e 240 return (TH2S*)fLinearFitterHistoArray.UncheckedAt(detector);
3a0f6479 241
fe4ee353 242 return GetLinearFitterHistoForce(detector);
243
244}
245//______________________________________________________________________________________
246TH2S* AliTRDCalibraVdriftLinearFit::GetLinearFitterHistoForce(Int_t detector)
247{
248 //
249 // return pointer to TH2F histo
250 // if NULL create a new histo if it doesn't exist allready
251 //
252 if (fLinearFitterHistoArray.UncheckedAt(detector))
253 return (TH2S*)fLinearFitterHistoArray.UncheckedAt(detector);
254
255 // if we are forced and TLinearFitter doesn't yes exist create it
256
257 // new TH2F
258 TString name("LFDV");
259 name += detector;
260 name += "version";
261 name += fVersion;
262
263 TH2S *lfdv = new TH2S((const Char_t *)name,(const Char_t *) name
264 ,36,-0.9,0.9,48
265 ,-1.2,1.2);
266 lfdv->SetXTitle("tan(phi_{track})");
267 lfdv->SetYTitle("dy/dt");
268 lfdv->SetZTitle("Number of clusters");
269 lfdv->SetStats(0);
270 lfdv->SetDirectory(0);
271
272 fLinearFitterHistoArray.AddAt(lfdv,detector);
273 return lfdv;
3a0f6479 274}
275//______________________________________________________________________________________
276Bool_t AliTRDCalibraVdriftLinearFit::GetParam(Int_t detector, TVectorD *param)
277{
278 //
279 // return param for this detector
280 //
281 if ( fLinearFitterPArray.UncheckedAt(detector) ){
282 const TVectorD *vectorP = (TVectorD*)fLinearFitterPArray.UncheckedAt(detector);
283 if(!param) param = new TVectorD(2);
284 for(Int_t k = 0; k < 2; k++){
285 (*param)[k] = (*vectorP)[k];
286 }
287 return kTRUE;
288 }
289 else return kFALSE;
290
291}
292//______________________________________________________________________________________
293Bool_t AliTRDCalibraVdriftLinearFit::GetError(Int_t detector, TVectorD *error)
294{
295 //
296 // return error for this detector
297 //
298 if ( fLinearFitterEArray.UncheckedAt(detector) ){
299 const TVectorD *vectorE = (TVectorD*)fLinearFitterEArray.UncheckedAt(detector);
300 if(!error) error = new TVectorD(3);
301 for(Int_t k = 0; k < 3; k++){
302 (*error)[k] = (*vectorE)[k];
303 }
304 return kTRUE;
305 }
306 else return kFALSE;
307
308}
309//______________________________________________________________________________________
310void AliTRDCalibraVdriftLinearFit::Update(Int_t detector, Float_t tnp, Float_t pars1)
311{
312 //
313 // Fill the 2D histos for debugging
314 //
ba94744a 315
0ca1720e 316 TH2S *h = ((TH2S *) GetLinearFitterHisto(detector,kTRUE));
ba94744a 317 Double_t nbentries = h->GetEntries();
0ca1720e 318 if(nbentries < 5*32767) h->Fill(tnp,pars1);
3a0f6479 319
320}
321//____________Functions fit Online CH2d________________________________________
322void AliTRDCalibraVdriftLinearFit::FillPEArray()
323{
324 //
325 // Fill fLinearFitterPArray and fLinearFitterEArray from inside
326 //
327
328
329 Int_t *arrayI = new Int_t[540];
330 for(Int_t k = 0; k< 540; k++){
331 arrayI[k] = 0;
332 }
333
334 // Loop over histos
335 for(Int_t cb = 0; cb < 540; cb++){
0ca1720e 336 const TH2S *linearfitterhisto = (TH2S*)fLinearFitterHistoArray.UncheckedAt(cb);
3a0f6479 337 //printf("Processing the detector cb %d we find %d\n",cb, (Bool_t) linearfitterhisto);
338
0ca1720e 339 if ( linearfitterhisto != 0 ){
3a0f6479 340
341 // Fill a linearfitter
342 TAxis *xaxis = linearfitterhisto->GetXaxis();
343 TAxis *yaxis = linearfitterhisto->GetYaxis();
344 TLinearFitter linearfitter = TLinearFitter(2,"pol1");
345 //printf("test\n");
a5dcf618 346 Double_t integral = linearfitterhisto->Integral();
347 //printf("Integral is %f\n",integral);
348 Bool_t securitybreaking = kFALSE;
349 if(TMath::Abs(integral-1199) < 0.00001) securitybreaking = kTRUE;
3a0f6479 350 for(Int_t ibinx = 0; ibinx < linearfitterhisto->GetNbinsX(); ibinx++){
351 for(Int_t ibiny = 0; ibiny < linearfitterhisto->GetNbinsY(); ibiny++){
8a606eb5 352 if(linearfitterhisto->GetBinContent(ibinx+1,ibiny+1)>0){
3a0f6479 353 Double_t x = xaxis->GetBinCenter(ibinx+1);
354 Double_t y = yaxis->GetBinCenter(ibiny+1);
a5dcf618 355
3a0f6479 356 for(Int_t k = 0; k < (Int_t)linearfitterhisto->GetBinContent(ibinx+1,ibiny+1); k++){
a5dcf618 357 if(!securitybreaking){
358 linearfitter.AddPoint(&x,y);
359 arrayI[cb]++;
360 }
361 else {
362 if(arrayI[cb]< 1198){
363 linearfitter.AddPoint(&x,y);
364 arrayI[cb]++;
365 }
366 }
3a0f6479 367 }
a5dcf618 368
3a0f6479 369 }
370 }
371 }
372
373 //printf("Find %d entries for the detector %d\n",arrayI[cb],cb);
374
375 // Eval the linear fitter
376 if(arrayI[cb]>10){
377 TVectorD *par = new TVectorD(2);
378 TVectorD pare = TVectorD(2);
379 TVectorD *parE = new TVectorD(3);
a5dcf618 380 //printf("Fit\n");
2a1a7b36 381 if(((fRobustFit) && (linearfitter.EvalRobust(0.8)==0)) || ((!fRobustFit) && (linearfitter.Eval()==0))) {
a5dcf618 382 //if((linearfitter.Eval()==0)) {
383 //printf("Take the param\n");
c1105918 384 linearfitter.GetParameters(*par);
a5dcf618 385 //printf("Done\n");
c1105918 386 //linearfitter.GetErrors(pare);
387 //Float_t ppointError = TMath::Sqrt(TMath::Abs(linearfitter.GetChisquare())/arrayI[cb]);
388 //(*parE)[0] = pare[0]*ppointError;
389 //(*parE)[1] = pare[1]*ppointError;
390
391 (*parE)[0] = 0.0;
392 (*parE)[1] = 0.0;
393 (*parE)[2] = (Double_t) arrayI[cb];
394 fLinearFitterPArray.AddAt(par,cb);
395 fLinearFitterEArray.AddAt(parE,cb);
396
397 //par->Print();
398 //parE->Print();
399 }
a5dcf618 400 //printf("Finish\n");
ba94744a 401 }
c1105918 402
07f9b7f2 403 //delete linearfitterhisto;
ba94744a 404
3a0f6479 405 }// if something
fdc15553 406
3a0f6479 407 }
a987273c 408
409 delete [] arrayI;
3a0f6479 410
d1e1a064 411}
c1105918 412
661c7be6 413//____________Functions fit Online CH2d________________________________________
414void AliTRDCalibraVdriftLinearFit::FillPEArray2()
415{
416 //
417 // Fill fFitterPArray and fFitterEArray from inside
418 //
419
420 // Loop over histos
421 TF1 *f1 = new TF1("f1","[0]+[1]*x",-1,1);
422
423 for(Int_t cb = 0; cb < 540; cb++){
424 //const TH2S *fitterhisto = (TH2S*)fLinearFitterHistoArray.UncheckedAt(cb);
425 TH2S *fitterhisto = (TH2S*)fLinearFitterHistoArray.UncheckedAt(cb);
426 //printf("Processing the detector cb %d we find %d\n",cb, (Bool_t) fitterhisto);
427
428 if ( fitterhisto != 0 ){
429
430 Int_t nEntries=0;
431 TGraphErrors *gg=DrawMS(fitterhisto,nEntries);
ce42378b 432 if (!gg) continue;
661c7be6 433 // Number of points of the TGraphErrors
532be2d4 434 //printf("Number of points %d for detector %d\n",gg->GetN(),cb);
435 if(gg->GetN() < fMinNpointsFit) {
661c7be6 436 if(gg) delete gg;
437 continue;
438 }
439 //printf("det %d, number of points %d, nEntries %d\n",cb,gg->GetN(),nEntries);
440 gg->Fit(f1,"Q0");
441
442 TVectorD *par = new TVectorD(2);
443 TVectorD *parE = new TVectorD(3);
444 (*parE)[0] = 0.0;
445 (*parE)[1] = 0.0;
446 (*parE)[2] = (Double_t) nEntries;
447 (*par)[0] = f1->GetParameter(0);
448 (*par)[1] = f1->GetParameter(1);
449 fLinearFitterPArray.AddAt(par,cb);
450 fLinearFitterEArray.AddAt(parE,cb);
532be2d4 451 //printf("Value %f for detector %d\n",(*par)[1],cb);
661c7be6 452
453 if(fDebugLevel==0) {
454 if(gg) delete gg;
455 }
456 else {
457 if(cb==fSeeDetector) {
458 gStyle->SetPalette(1);
459 gStyle->SetOptStat(1111);
460 gStyle->SetPadBorderMode(0);
461 gStyle->SetCanvasColor(10);
462 gStyle->SetPadLeftMargin(0.13);
463 gStyle->SetPadRightMargin(0.01);
464 TCanvas *stat = new TCanvas("stat","",50,50,600,800);
465 stat->cd(1);
466 fitterhisto->Draw("colz");
467 gg->Draw("P");
468 f1->Draw("same");
469 break;
470 }
471 }
472 }
473 }
474 if(fDebugLevel==0) delete f1;
475}
476
477//_________Helper function__________________________________________________
478TGraphErrors* AliTRDCalibraVdriftLinearFit::DrawMS(const TH2 *const h2, Int_t &nEntries)
479{
92df0132 480 //
481 // Debug function
482 //
483
661c7be6 484 TF1 fg("fg", "gaus", -10., 30.);
485 TGraphErrors *gp = new TGraphErrors();
c1105918 486
661c7be6 487 TAxis *ax(h2->GetXaxis());
488 TAxis *ay(h2->GetYaxis());
489 TH1D *h1(NULL);
490 for(Int_t ipt(0), jpt(1), ig(0); ipt<ax->GetNbins(); ipt++, jpt++){
491 h1 = h2->ProjectionY("py", jpt, jpt);
492 fg.SetParameter(1, h1->GetMean());
493 //Float_t x(ax->GetBinCenter(jpt));
494 Int_t n=Int_t(h1->Integral(1, h1->GetNbinsX()));
495 nEntries+=n;
496 if(n < 15){
497 //Warning("drawMS()", Form("reject x[%d]=%f on n=%d", jpt, x, n));
498 continue;
499 }
500 h1->Fit(&fg, "WWQ0");
501 if(fg.GetNDF()<2){
502 //Warning("drawMS()", Form("reject x[%d]=%f on NDF=%d", jpt, x, fg.GetNDF()));
503 continue;
504 }
64c7fdf8 505 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 506 gp->SetPoint(ig, ax->GetBinCenter(jpt), fg.GetParameter(1));
507 gp->SetPointError(ig, 0, TMath::Sqrt(pow(fg.GetParError(1),2) + (1/pow(fg.GetParameter(0),2))));
508 ig++;
509 }
510 delete h1;
511 return gp;
512}