Removing useless flag.
[u/mrichter/AliRoot.git] / TRD / AliTRDCalibraVdriftLinearFit.cxx
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 ClassImp(AliTRDCalibraVdriftLinearFit) /*FOLD00*/
49
50 //_____________________________________________________________________
51 AliTRDCalibraVdriftLinearFit::AliTRDCalibraVdriftLinearFit() : /*FOLD00*/
52   TObject(),
53   fVersion(0),
54   fNameCalibUsed(""),
55   fLinearFitterHistoArray(540),
56   fLinearFitterPArray(540),
57   fLinearFitterEArray(540),
58   fRobustFit(kTRUE),
59   fMinNpointsFit(11),
60   fNbBindx(32),
61   fNbBindy(70),
62   fRangedx(0.8),
63   fRangedy(1.4),
64   fDebugStreamer(0x0),
65   fDebugLevel(0),
66   fSeeDetector(0)
67 {
68   //
69   // default constructor
70   //
71 }
72 //_____________________________________________________________________
73 AliTRDCalibraVdriftLinearFit::AliTRDCalibraVdriftLinearFit(const AliTRDCalibraVdriftLinearFit &ped) : /*FOLD00*/
74   TObject(ped),
75   fVersion(ped.fVersion),
76   fNameCalibUsed(ped.fNameCalibUsed),
77   fLinearFitterHistoArray(540),
78   fLinearFitterPArray(540),
79   fLinearFitterEArray(540),
80   fRobustFit(kTRUE),
81   fMinNpointsFit(10),
82   fNbBindx(32),
83   fNbBindy(70),
84   fRangedx(0.8),
85   fRangedy(1.4),
86   fDebugStreamer(0x0),
87   fDebugLevel(0),
88   fSeeDetector(0)
89 {
90     //
91     // copy constructor
92     //
93   for (Int_t idet = 0; idet < 540; idet++){
94    
95     const TVectorD     *vectorE     = (TVectorD*)ped.fLinearFitterEArray.UncheckedAt(idet);
96     const TVectorD     *vectorP     = (TVectorD*)ped.fLinearFitterPArray.UncheckedAt(idet);
97     const TH2S         *hped        = (TH2S*)ped.fLinearFitterHistoArray.UncheckedAt(idet);
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 ){
102       TH2S *hNew = (TH2S *)hped->Clone();
103       //hNew->SetDirectory(0);
104       fLinearFitterHistoArray.AddAt(hNew,idet);
105     }
106   }
107 }
108 //_____________________________________________________________________
109 AliTRDCalibraVdriftLinearFit::AliTRDCalibraVdriftLinearFit(const TObjArray &obja) : /*FOLD00*/
110   TObject(),
111   fVersion(0),
112   fNameCalibUsed(""),
113   fLinearFitterHistoArray(540),
114   fLinearFitterPArray(540),
115   fLinearFitterEArray(540),
116   fRobustFit(kTRUE),
117   fMinNpointsFit(10),
118   fNbBindx(32),
119   fNbBindy(70),
120   fRangedx(0.8),
121   fRangedy(1.4),
122   fDebugStreamer(0x0),
123   fDebugLevel(0),
124   fSeeDetector(0)
125 {
126   //
127   // constructor from a TObjArray
128   //
129   for (Int_t idet = 0; idet < 540; idet++){
130     const TH2S         *hped        = (TH2S*)obja.UncheckedAt(idet);
131     if ( hped != 0x0 ){
132       TH2S *hNew = (TH2S *)hped->Clone();
133       //hNew->SetDirectory(0);
134       fLinearFitterHistoArray.AddAt(hNew,idet);
135     }
136   }
137 }
138 //_____________________________________________________________________
139 AliTRDCalibraVdriftLinearFit& 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 //_____________________________________________________________________
150 AliTRDCalibraVdriftLinearFit::~AliTRDCalibraVdriftLinearFit() /*FOLD00*/
151 {
152   //
153   // destructor
154   //
155   fLinearFitterHistoArray.SetOwner();
156   fLinearFitterPArray.SetOwner();
157   fLinearFitterEArray.SetOwner();
158
159   fLinearFitterHistoArray.Delete();
160   fLinearFitterPArray.Delete();
161   fLinearFitterEArray.Delete();
162
163   if ( fDebugStreamer ) delete fDebugStreamer;
164
165 }
166 //_____________________________________________________________________________
167 void 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)){
178       TH2S *hped1 = (TH2S *)target.GetLinearFitterHisto(idet,kTRUE);
179       //hped1->SetDirectory(0);
180       hped1->Add((const TH2S *)fLinearFitterHistoArray.UncheckedAt(idet));
181     }
182   }
183   
184   TObject::Copy(c);
185
186 }
187 //_____________________________________________________________________________
188 Long64_t AliTRDCalibraVdriftLinearFit::Merge(const TCollection* list) 
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) 
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)){
211           TH2S *hped1 = (TH2S *)GetLinearFitterHisto(idet,kTRUE);
212           Double_t entriesa = hped1->GetEntries();
213           Double_t entriesb = ((TH2S *)entry->GetLinearFitterHisto(idet))->GetEntries();
214           if((entriesa + entriesb) < 5*32767) hped1->Add(entry->GetLinearFitterHisto(idet));
215         }
216       }
217       
218       count++;
219     }
220   
221   return count;
222 }
223 //_____________________________________________________________________
224 void AliTRDCalibraVdriftLinearFit::Add(const AliTRDCalibraVdriftLinearFit *ped)
225 {
226   //
227   // Add histo
228   //
229
230   fVersion++;
231
232   for (Int_t idet = 0; idet < 540; idet++){
233     const TH2S         *hped        = (TH2S*)ped->GetLinearFitterHistoNoForce(idet);
234     //printf("idet %d\n",idet);
235     if ( hped != 0x0 ){
236       //printf("add\n");
237       TH2S *hped1 = (TH2S *)GetLinearFitterHisto(idet,kTRUE);
238       Double_t entriesa = hped1->GetEntries();
239       Double_t entriesb = hped->GetEntries();
240       if((entriesa + entriesb) < 5*32767) hped1->Add(hped);
241     }
242   }
243 }
244 //______________________________________________________________________________________
245 TH2S* 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 //______________________________________________________________________________________
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     if ( !force || fLinearFitterHistoArray.UncheckedAt(detector) )
272         return (TH2S*)fLinearFitterHistoArray.UncheckedAt(detector);
273
274     return GetLinearFitterHistoForce(detector);
275
276 }
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   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
296                         ,(Int_t)fNbBindx,-fRangedx,fRangedx,(Int_t)fNbBindy
297                         ,-fRangedy,fRangedy);
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;
306 }
307 //______________________________________________________________________________________
308 Bool_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 //______________________________________________________________________________________
325 Bool_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 //______________________________________________________________________________________
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   TH2S *h = ((TH2S *) GetLinearFitterHisto(detector,kTRUE));
349   Double_t nbentries = h->GetEntries();
350   if(nbentries < 5*32767) h->Fill(tnp,pars1);
351
352 }
353 //____________Functions fit Online CH2d________________________________________
354 void 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++){
368     const TH2S *linearfitterhisto = (TH2S*)fLinearFitterHistoArray.UncheckedAt(cb);
369     //printf("Processing the detector cb %d we find %d\n",cb, (Bool_t) linearfitterhisto);    
370
371     if ( linearfitterhisto != 0 ){
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");
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;
382       for(Int_t ibinx = 0; ibinx < linearfitterhisto->GetNbinsX(); ibinx++){
383         for(Int_t ibiny = 0; ibiny < linearfitterhisto->GetNbinsY(); ibiny++){
384           if(linearfitterhisto->GetBinContent(ibinx+1,ibiny+1)>0){
385             Double_t x = xaxis->GetBinCenter(ibinx+1);
386             Double_t y = yaxis->GetBinCenter(ibiny+1);
387             
388             for(Int_t k = 0; k < (Int_t)linearfitterhisto->GetBinContent(ibinx+1,ibiny+1); k++){
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               }
399             }
400             
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);
412         //printf("Fit\n");
413         if(((fRobustFit) && (linearfitter.EvalRobust(0.8)==0)) || ((!fRobustFit) && (linearfitter.Eval()==0))) {
414           //if((linearfitter.Eval()==0)) {
415           //printf("Take the param\n");
416           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           (*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         }
432         //printf("Finish\n");
433       }
434       
435       //delete linearfitterhisto;
436       
437     }// if something
438
439   }
440
441   delete [] arrayI;
442    
443 }
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   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);
464       if (!gg) continue;
465       // Number of points of the TGraphErrors
466       //printf("Number of points %d for detector %d\n",gg->GetN(),cb);
467       if(gg->GetN() <  fMinNpointsFit) {
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);
483       //printf("Value %f for detector %d\n",(*par)[1],cb);
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__________________________________________________
510 TGraphErrors* AliTRDCalibraVdriftLinearFit::DrawMS(const TH2 *const h2, Int_t &nEntries)
511 {
512   //
513   // Debug function
514   //
515
516   TF1 fg("fg", "gaus", -10., 30.);
517   TGraphErrors *gp = new TGraphErrors();
518
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     }
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;
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 }