Be sure to load mapping when needed
[u/mrichter/AliRoot.git] / TRD / AliTRDCalibraVdriftLinearFit.cxx
CommitLineData
3a0f6479 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
0bc7827a 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
3a0f6479 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 <TTreeStream.h>
ba94744a 38
3a0f6479 39
40//header file
41#include "AliTRDCalibraVdriftLinearFit.h"
42
43ClassImp(AliTRDCalibraVdriftLinearFit) /*FOLD00*/
44
45//_____________________________________________________________________
46AliTRDCalibraVdriftLinearFit::AliTRDCalibraVdriftLinearFit() : /*FOLD00*/
47 TObject(),
48 fVersion(0),
49 fLinearFitterHistoArray(540),
50 fLinearFitterPArray(540),
51 fLinearFitterEArray(540)
52{
ba94744a 53 //
54 // default constructor
55 //
3a0f6479 56}
3a0f6479 57//_____________________________________________________________________
58AliTRDCalibraVdriftLinearFit::AliTRDCalibraVdriftLinearFit(const AliTRDCalibraVdriftLinearFit &ped) : /*FOLD00*/
59 TObject(ped),
60 fVersion(ped.fVersion),
61 fLinearFitterHistoArray(540),
62 fLinearFitterPArray(540),
63 fLinearFitterEArray(540)
64{
65 //
66 // copy constructor
67 //
68 for (Int_t idet = 0; idet < 540; idet++){
ba94744a 69
3a0f6479 70 const TVectorD *vectorE = (TVectorD*)ped.fLinearFitterEArray.UncheckedAt(idet);
71 const TVectorD *vectorP = (TVectorD*)ped.fLinearFitterPArray.UncheckedAt(idet);
0ca1720e 72 const TH2S *hped = (TH2S*)ped.fLinearFitterHistoArray.UncheckedAt(idet);
3a0f6479 73
74 if ( vectorE != 0x0 ) fLinearFitterEArray.AddAt(new TVectorD(*vectorE), idet);
75 if ( vectorP != 0x0 ) fLinearFitterPArray.AddAt(new TVectorD(*vectorP), idet);
76 if ( hped != 0x0 ){
0ca1720e 77 TH2S *hNew = (TH2S *)hped->Clone();
ba94744a 78 //hNew->SetDirectory(0);
3a0f6479 79 fLinearFitterHistoArray.AddAt(hNew,idet);
80 }
3a0f6479 81 }
82}
83//_____________________________________________________________________
d0569428 84AliTRDCalibraVdriftLinearFit::AliTRDCalibraVdriftLinearFit(const TObjArray &obja) : /*FOLD00*/
85 TObject(),
86 fVersion(0),
87 fLinearFitterHistoArray(540),
88 fLinearFitterPArray(540),
89 fLinearFitterEArray(540)
90{
91 //
92 // constructor from a TObjArray
93 //
94 for (Int_t idet = 0; idet < 540; idet++){
0ca1720e 95 const TH2S *hped = (TH2S*)obja.UncheckedAt(idet);
d0569428 96 if ( hped != 0x0 ){
0ca1720e 97 TH2S *hNew = (TH2S *)hped->Clone();
ba94744a 98 //hNew->SetDirectory(0);
d0569428 99 fLinearFitterHistoArray.AddAt(hNew,idet);
100 }
101 }
102}
103//_____________________________________________________________________
3a0f6479 104AliTRDCalibraVdriftLinearFit& AliTRDCalibraVdriftLinearFit::operator = (const AliTRDCalibraVdriftLinearFit &source)
105{
106 //
107 // assignment operator
108 //
109 if (&source == this) return *this;
110 new (this) AliTRDCalibraVdriftLinearFit(source);
111
112 return *this;
113}
114//_____________________________________________________________________
115AliTRDCalibraVdriftLinearFit::~AliTRDCalibraVdriftLinearFit() /*FOLD00*/
116{
117 //
118 // destructor
119 //
f5cda287 120 fLinearFitterHistoArray.SetOwner();
121 fLinearFitterPArray.SetOwner();
122 fLinearFitterEArray.SetOwner();
123
1ca79a00 124 fLinearFitterHistoArray.Delete();
125 fLinearFitterPArray.Delete();
126 fLinearFitterEArray.Delete();
127
3a0f6479 128}
64942b85 129//_____________________________________________________________________________
130void AliTRDCalibraVdriftLinearFit::Copy(TObject &c) const
131{
132 //
133 // Copy function
134 //
135
136 AliTRDCalibraVdriftLinearFit& target = (AliTRDCalibraVdriftLinearFit &) c;
137
138 // Copy only the histos
139 for (Int_t idet = 0; idet < 540; idet++){
140 if(fLinearFitterHistoArray.UncheckedAt(idet)){
0ca1720e 141 TH2S *hped1 = (TH2S *)target.GetLinearFitterHisto(idet,kTRUE);
ba94744a 142 //hped1->SetDirectory(0);
0ca1720e 143 hped1->Add((const TH2S *)fLinearFitterHistoArray.UncheckedAt(idet));
64942b85 144 }
64942b85 145 }
ba94744a 146
64942b85 147 TObject::Copy(c);
148
149}
150//_____________________________________________________________________________
6bbdc11a 151Long64_t AliTRDCalibraVdriftLinearFit::Merge(const TCollection* list)
64942b85 152{
153 // Merge list of objects (needed by PROOF)
154
155 if (!list)
156 return 0;
157
158 if (list->IsEmpty())
159 return 1;
160
161 TIterator* iter = list->MakeIterator();
162 TObject* obj = 0;
163
164 // collection of generated histograms
165 Int_t count=0;
166 while((obj = iter->Next()) != 0)
ba94744a 167 {
168 AliTRDCalibraVdriftLinearFit* entry = dynamic_cast<AliTRDCalibraVdriftLinearFit*>(obj);
169 if (entry == 0) continue;
170
171 // Copy only the histos
172 for (Int_t idet = 0; idet < 540; idet++){
173 if(entry->GetLinearFitterHisto(idet)){
0ca1720e 174 TH2S *hped1 = (TH2S *)GetLinearFitterHisto(idet,kTRUE);
ba94744a 175 Double_t entriesa = hped1->GetEntries();
0ca1720e 176 Double_t entriesb = ((TH2S *)entry->GetLinearFitterHisto(idet))->GetEntries();
ba94744a 177 if((entriesa + entriesb) < 5*32767) hped1->Add(entry->GetLinearFitterHisto(idet));
178 }
64942b85 179 }
180
ba94744a 181 count++;
64942b85 182 }
64942b85 183
184 return count;
185}
3a0f6479 186//_____________________________________________________________________
fe4ee353 187void AliTRDCalibraVdriftLinearFit::Add(const AliTRDCalibraVdriftLinearFit *ped)
3a0f6479 188{
189 //
190 // Add histo
191 //
192
193 fVersion++;
194
195 for (Int_t idet = 0; idet < 540; idet++){
fe4ee353 196 const TH2S *hped = (TH2S*)ped->GetLinearFitterHistoNoForce(idet);
3a0f6479 197 //printf("idet %d\n",idet);
198 if ( hped != 0x0 ){
199 //printf("add\n");
0ca1720e 200 TH2S *hped1 = (TH2S *)GetLinearFitterHisto(idet,kTRUE);
ba94744a 201 Double_t entriesa = hped1->GetEntries();
202 Double_t entriesb = hped->GetEntries();
203 if((entriesa + entriesb) < 5*32767) hped1->Add(hped);
3a0f6479 204 }
205 }
206}
207//______________________________________________________________________________________
0ca1720e 208TH2S* AliTRDCalibraVdriftLinearFit::GetLinearFitterHisto(Int_t detector, Bool_t force)
3a0f6479 209{
210 //
211 // return pointer to TH2F histo
212 // if force is true create a new histo if it doesn't exist allready
213 //
214 if ( !force || fLinearFitterHistoArray.UncheckedAt(detector) )
0ca1720e 215 return (TH2S*)fLinearFitterHistoArray.UncheckedAt(detector);
3a0f6479 216
fe4ee353 217 return GetLinearFitterHistoForce(detector);
218
219}
220//______________________________________________________________________________________
221TH2S* AliTRDCalibraVdriftLinearFit::GetLinearFitterHistoForce(Int_t detector)
222{
223 //
224 // return pointer to TH2F histo
225 // if NULL create a new histo if it doesn't exist allready
226 //
227 if (fLinearFitterHistoArray.UncheckedAt(detector))
228 return (TH2S*)fLinearFitterHistoArray.UncheckedAt(detector);
229
230 // if we are forced and TLinearFitter doesn't yes exist create it
231
232 // new TH2F
233 TString name("LFDV");
234 name += detector;
235 name += "version";
236 name += fVersion;
237
238 TH2S *lfdv = new TH2S((const Char_t *)name,(const Char_t *) name
239 ,36,-0.9,0.9,48
240 ,-1.2,1.2);
241 lfdv->SetXTitle("tan(phi_{track})");
242 lfdv->SetYTitle("dy/dt");
243 lfdv->SetZTitle("Number of clusters");
244 lfdv->SetStats(0);
245 lfdv->SetDirectory(0);
246
247 fLinearFitterHistoArray.AddAt(lfdv,detector);
248 return lfdv;
3a0f6479 249}
250//______________________________________________________________________________________
251Bool_t AliTRDCalibraVdriftLinearFit::GetParam(Int_t detector, TVectorD *param)
252{
253 //
254 // return param for this detector
255 //
256 if ( fLinearFitterPArray.UncheckedAt(detector) ){
257 const TVectorD *vectorP = (TVectorD*)fLinearFitterPArray.UncheckedAt(detector);
258 if(!param) param = new TVectorD(2);
259 for(Int_t k = 0; k < 2; k++){
260 (*param)[k] = (*vectorP)[k];
261 }
262 return kTRUE;
263 }
264 else return kFALSE;
265
266}
267//______________________________________________________________________________________
268Bool_t AliTRDCalibraVdriftLinearFit::GetError(Int_t detector, TVectorD *error)
269{
270 //
271 // return error for this detector
272 //
273 if ( fLinearFitterEArray.UncheckedAt(detector) ){
274 const TVectorD *vectorE = (TVectorD*)fLinearFitterEArray.UncheckedAt(detector);
275 if(!error) error = new TVectorD(3);
276 for(Int_t k = 0; k < 3; k++){
277 (*error)[k] = (*vectorE)[k];
278 }
279 return kTRUE;
280 }
281 else return kFALSE;
282
283}
284//______________________________________________________________________________________
285void AliTRDCalibraVdriftLinearFit::Update(Int_t detector, Float_t tnp, Float_t pars1)
286{
287 //
288 // Fill the 2D histos for debugging
289 //
ba94744a 290
0ca1720e 291 TH2S *h = ((TH2S *) GetLinearFitterHisto(detector,kTRUE));
ba94744a 292 Double_t nbentries = h->GetEntries();
0ca1720e 293 if(nbentries < 5*32767) h->Fill(tnp,pars1);
3a0f6479 294
295}
296//____________Functions fit Online CH2d________________________________________
297void AliTRDCalibraVdriftLinearFit::FillPEArray()
298{
299 //
300 // Fill fLinearFitterPArray and fLinearFitterEArray from inside
301 //
302
303
304 Int_t *arrayI = new Int_t[540];
305 for(Int_t k = 0; k< 540; k++){
306 arrayI[k] = 0;
307 }
308
309 // Loop over histos
310 for(Int_t cb = 0; cb < 540; cb++){
0ca1720e 311 const TH2S *linearfitterhisto = (TH2S*)fLinearFitterHistoArray.UncheckedAt(cb);
3a0f6479 312 //printf("Processing the detector cb %d we find %d\n",cb, (Bool_t) linearfitterhisto);
313
0ca1720e 314 if ( linearfitterhisto != 0 ){
3a0f6479 315
316 // Fill a linearfitter
317 TAxis *xaxis = linearfitterhisto->GetXaxis();
318 TAxis *yaxis = linearfitterhisto->GetYaxis();
319 TLinearFitter linearfitter = TLinearFitter(2,"pol1");
320 //printf("test\n");
a5dcf618 321 Double_t integral = linearfitterhisto->Integral();
322 //printf("Integral is %f\n",integral);
323 Bool_t securitybreaking = kFALSE;
324 if(TMath::Abs(integral-1199) < 0.00001) securitybreaking = kTRUE;
3a0f6479 325 for(Int_t ibinx = 0; ibinx < linearfitterhisto->GetNbinsX(); ibinx++){
326 for(Int_t ibiny = 0; ibiny < linearfitterhisto->GetNbinsY(); ibiny++){
8a606eb5 327 if(linearfitterhisto->GetBinContent(ibinx+1,ibiny+1)>0){
3a0f6479 328 Double_t x = xaxis->GetBinCenter(ibinx+1);
329 Double_t y = yaxis->GetBinCenter(ibiny+1);
a5dcf618 330
3a0f6479 331 for(Int_t k = 0; k < (Int_t)linearfitterhisto->GetBinContent(ibinx+1,ibiny+1); k++){
a5dcf618 332 if(!securitybreaking){
333 linearfitter.AddPoint(&x,y);
334 arrayI[cb]++;
335 }
336 else {
337 if(arrayI[cb]< 1198){
338 linearfitter.AddPoint(&x,y);
339 arrayI[cb]++;
340 }
341 }
3a0f6479 342 }
a5dcf618 343
3a0f6479 344 }
345 }
346 }
347
348 //printf("Find %d entries for the detector %d\n",arrayI[cb],cb);
349
350 // Eval the linear fitter
351 if(arrayI[cb]>10){
352 TVectorD *par = new TVectorD(2);
353 TVectorD pare = TVectorD(2);
354 TVectorD *parE = new TVectorD(3);
a5dcf618 355 //printf("Fit\n");
c1105918 356 if((linearfitter.EvalRobust(0.8)==0)) {
a5dcf618 357 //if((linearfitter.Eval()==0)) {
358 //printf("Take the param\n");
c1105918 359 linearfitter.GetParameters(*par);
a5dcf618 360 //printf("Done\n");
c1105918 361 //linearfitter.GetErrors(pare);
362 //Float_t ppointError = TMath::Sqrt(TMath::Abs(linearfitter.GetChisquare())/arrayI[cb]);
363 //(*parE)[0] = pare[0]*ppointError;
364 //(*parE)[1] = pare[1]*ppointError;
365
366 (*parE)[0] = 0.0;
367 (*parE)[1] = 0.0;
368 (*parE)[2] = (Double_t) arrayI[cb];
369 fLinearFitterPArray.AddAt(par,cb);
370 fLinearFitterEArray.AddAt(parE,cb);
371
372 //par->Print();
373 //parE->Print();
374 }
a5dcf618 375 //printf("Finish\n");
ba94744a 376 }
c1105918 377
07f9b7f2 378 //delete linearfitterhisto;
ba94744a 379
3a0f6479 380 }// if something
fdc15553 381
3a0f6479 382 }
a987273c 383
384 delete [] arrayI;
3a0f6479 385
d1e1a064 386}
c1105918 387
388