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