Improved quitting Event Display with close button.
[u/mrichter/AliRoot.git] / TRD / AliTRDCalibraExbAltFit.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: 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 ClassImp(AliTRDCalibraExbAltFit) /*FOLD00*/
46
47 //_____________________________________________________________________
48 AliTRDCalibraExbAltFit::AliTRDCalibraExbAltFit() : /*FOLD00*/
49   TObject(),
50   fVersion(0),
51   fFitterHistoArray(540),
52   fFitterPArray(540),
53   fFitterEArray(540),
54   fRobustFit(kFALSE),
55   fDebugStreamer(0x0),
56   fDebugLevel(0)
57 {
58   //
59   // default constructor
60   //
61 }
62 //_____________________________________________________________________
63 AliTRDCalibraExbAltFit::AliTRDCalibraExbAltFit(const AliTRDCalibraExbAltFit &ped) : /*FOLD00*/
64   TObject(ped),
65   fVersion(ped.fVersion),
66   fFitterHistoArray(540),
67   fFitterPArray(540),
68   fFitterEArray(540),
69   fRobustFit(kFALSE),
70   fDebugStreamer(0x0),
71   fDebugLevel(0)
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 //_____________________________________________________________________
92 AliTRDCalibraExbAltFit::AliTRDCalibraExbAltFit(const TObjArray &obja) : /*FOLD00*/
93   TObject(),
94   fVersion(0),
95   fFitterHistoArray(540),
96   fFitterPArray(540),
97   fFitterEArray(540),
98   fRobustFit(kFALSE),
99   fDebugStreamer(0x0),
100   fDebugLevel(0)
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 //_____________________________________________________________________
115 AliTRDCalibraExbAltFit& 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 //_____________________________________________________________________
126 AliTRDCalibraExbAltFit::~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
139   if ( fDebugStreamer ) delete fDebugStreamer;
140
141 }
142 //_____________________________________________________________________________
143 void 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 //_____________________________________________________________________________
164 Long64_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 //_____________________________________________________________________
200 void 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 //______________________________________________________________________________________
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     if ( !force || fFitterHistoArray.UncheckedAt(detector) )
228         return (TH2S*)fFitterHistoArray.UncheckedAt(detector);
229
230     return GetFitterHistoForce(detector);
231
232 }
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   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 //______________________________________________________________________________________
264 Bool_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 //______________________________________________________________________________________
281 Bool_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 //______________________________________________________________________________________
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   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________________________________________
310 void 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       const TAxis *xaxis = fitterhisto->GetXaxis();
331       const 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");
369         //if((fitter.EvalRobust(0.8)==0)) {
370         //if((fitter.EvalRobust()==0)) {
371         if(((fRobustFit) && (fitter.EvalRobust(0.8)==0)) || ((!fRobustFit) && (fitter.Eval()==0))) {
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
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   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);
420       if (!gg) continue;
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__________________________________________________
473 TGraphErrors* AliTRDCalibraExbAltFit::DrawMS(const TH2 *const h2, Int_t &nEntries)
474 {
475   //
476   // Debug function
477   //
478
479
480   TF1 fg("fg", "gaus", -10., 30.);
481   TGraphErrors *gp = new TGraphErrors();
482
483   const TAxis *ax(h2->GetXaxis());
484   const 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     }
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;
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 }