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