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