]> git.uio.no Git - u/mrichter/AliRoot.git/blame_incremental - TRD/AliTRDCalibraFit.cxx
git://git.uio.no / u / mrichter / AliRoot.git / blame_incremental
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (non-incremental) | history | HEAD
Fixing some documentation and comments.
[u/mrichter/AliRoot.git] / TRD / AliTRDCalibraFit.cxx
... / ...
CommitLineData
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// AliTRDCalibraFit
21//
22// This class is for the TRD calibration of the relative gain factor, the drift velocity,
23// the time 0 and the pad response function. It fits the histos.
24// The 2D histograms or vectors (first converted in 1D histos) will be fitted
25// if they have enough entries, otherwise the (default) value of the choosen database
26// will be put. For the relative gain calibration the resulted factors will be globally
27// normalized to the gain factors of the choosen database. It unables to precise
28// previous calibration procedure.
29// The function SetDebug enables the user to see:
30// _fDebug = 0: nothing, only the values are written in the tree if wanted
31// _fDebug = 1: a comparaison of the coefficients found and the default values
32// in the choosen database.
33// fCoef , histogram of the coefs as function of the calibration group number
34// fDelta , histogram of the relative difference of the coef with the default
35// value in the database as function of the calibration group number
36// fError , dirstribution of this relative difference
37// _fDebug = 2: only the fit of the choosen calibration group fFitVoir (SetFitVoir)
38// _fDebug = 3: The coefficients in the choosen detector fDet (SetDet) as function of the
39// pad row and col number
40// _fDebug = 4; The coeffcicients in the choosen detector fDet (SetDet) like in the 3 but with
41// also the comparaison histograms of the 1 for this detector
42//
43//
44// Author:
45// R. Bailhache (R.Bailhache@gsi.de)
46//
47//////////////////////////////////////////////////////////////////////////////////////
48
49#include <TLine.h>
50#include <TH1I.h>
51#include <TStyle.h>
52#include <TProfile2D.h>
53#include <TCanvas.h>
54#include <TGraphErrors.h>
55#include <TObjArray.h>
56#include <TH1.h>
57#include <TH1F.h>
58#include <TF1.h>
59#include <TAxis.h>
60#include <TMath.h>
61#include <TDirectory.h>
62#include <TROOT.h>
63#include <TTreeStream.h>
64#include <TLinearFitter.h>
65#include <TVectorD.h>
66#include <TArrayF.h>
67
68#include "AliLog.h"
69#include "AliMathBase.h"
70
71#include "AliTRDCalibraFit.h"
72#include "AliTRDCalibraMode.h"
73#include "AliTRDCalibraVector.h"
74#include "AliTRDCalibraVdriftLinearFit.h"
75#include "AliTRDcalibDB.h"
76#include "AliTRDgeometry.h"
77#include "AliTRDpadPlane.h"
78#include "AliTRDgeometry.h"
79#include "./Cal/AliTRDCalROC.h"
80#include "./Cal/AliTRDCalPad.h"
81#include "./Cal/AliTRDCalDet.h"
82
83
84ClassImp(AliTRDCalibraFit)
85
86AliTRDCalibraFit* AliTRDCalibraFit::fgInstance = 0;
87Bool_t AliTRDCalibraFit::fgTerminated = kFALSE;
88
89//_____________singleton implementation_________________________________________________
90AliTRDCalibraFit *AliTRDCalibraFit::Instance()
91{
92 //
93 // Singleton implementation
94 //
95
96 if (fgTerminated != kFALSE) {
97 return 0;
98 }
99
100 if (fgInstance == 0) {
101 fgInstance = new AliTRDCalibraFit();
102 }
103
104 return fgInstance;
105
106}
107
108//______________________________________________________________________________________
109void AliTRDCalibraFit::Terminate()
110{
111 //
112 // Singleton implementation
113 // Deletes the instance of this class
114 //
115
116 fgTerminated = kTRUE;
117
118 if (fgInstance != 0) {
119 delete fgInstance;
120 fgInstance = 0;
121 }
122
123}
124
125//______________________________________________________________________________________
126AliTRDCalibraFit::AliTRDCalibraFit()
127 :TObject()
128 ,fGeo(0)
129 ,fNumberOfBinsExpected(0)
130 ,fMethod(0)
131 ,fBeginFitCharge(3.5)
132 ,fFitPHPeriode(1)
133 ,fTakeTheMaxPH(kTRUE)
134 ,fT0Shift0(0.124797)
135 ,fT0Shift1(0.267451)
136 ,fRangeFitPRF(1.0)
137 ,fAccCDB(kFALSE)
138 ,fMinEntries(800)
139 ,fRebin(1)
140 ,fNumberFit(0)
141 ,fNumberFitSuccess(0)
142 ,fNumberEnt(0)
143 ,fStatisticMean(0.0)
144 ,fDebugStreamer(0x0)
145 ,fDebugLevel(0)
146 ,fFitVoir(0)
147 ,fMagneticField(0.5)
148 ,fCalibraMode(new AliTRDCalibraMode())
149 ,fCurrentCoefE(0.0)
150 ,fCurrentCoefE2(0.0)
151 ,fDect1(0)
152 ,fDect2(0)
153 ,fScaleFitFactor(0.0)
154 ,fEntriesCurrent(0)
155 ,fCountDet(0)
156 ,fCount(0)
157 ,fCalDet(0x0)
158 ,fCalROC(0x0)
159 ,fCalDet2(0x0)
160 ,fCalROC2(0x0)
161 ,fCurrentCoefDetector(0x0)
162 ,fCurrentCoefDetector2(0x0)
163 ,fVectorFit(0)
164 ,fVectorFit2(0)
165{
166 //
167 // Default constructor
168 //
169
170 fGeo = new AliTRDgeometry();
171
172 // Current variables initialised
173 for (Int_t k = 0; k < 2; k++) {
174 fCurrentCoef[k] = 0.0;
175 fCurrentCoef2[k] = 0.0;
176 }
177 for (Int_t i = 0; i < 3; i++) {
178 fPhd[i] = 0.0;
179 fDet[i] = 0;
180 }
181
182}
183//______________________________________________________________________________________
184AliTRDCalibraFit::AliTRDCalibraFit(const AliTRDCalibraFit &c)
185:TObject(c)
186,fGeo(0)
187,fNumberOfBinsExpected(c.fNumberOfBinsExpected)
188,fMethod(c.fMethod)
189,fBeginFitCharge(c.fBeginFitCharge)
190,fFitPHPeriode(c.fFitPHPeriode)
191,fTakeTheMaxPH(c.fTakeTheMaxPH)
192,fT0Shift0(c.fT0Shift0)
193,fT0Shift1(c.fT0Shift1)
194,fRangeFitPRF(c.fRangeFitPRF)
195,fAccCDB(c.fAccCDB)
196,fMinEntries(c.fMinEntries)
197,fRebin(c.fRebin)
198,fNumberFit(c.fNumberFit)
199,fNumberFitSuccess(c.fNumberFitSuccess)
200,fNumberEnt(c.fNumberEnt)
201,fStatisticMean(c.fStatisticMean)
202,fDebugStreamer(0x0)
203,fDebugLevel(c.fDebugLevel)
204,fFitVoir(c.fFitVoir)
205,fMagneticField(c.fMagneticField)
206,fCalibraMode(0x0)
207,fCurrentCoefE(c.fCurrentCoefE)
208,fCurrentCoefE2(c.fCurrentCoefE2)
209,fDect1(c.fDect1)
210,fDect2(c.fDect2)
211,fScaleFitFactor(c.fScaleFitFactor)
212,fEntriesCurrent(c.fEntriesCurrent)
213,fCountDet(c.fCountDet)
214,fCount(c.fCount)
215,fCalDet(0x0)
216,fCalROC(0x0)
217,fCalDet2(0x0)
218,fCalROC2(0x0)
219,fCurrentCoefDetector(0x0)
220,fCurrentCoefDetector2(0x0)
221,fVectorFit(0)
222,fVectorFit2(0)
223{
224 //
225 // Copy constructor
226 //
227
228 if(c.fCalibraMode) fCalibraMode = new AliTRDCalibraMode(*c.fCalibraMode);
229
230 //Current variables initialised
231 for (Int_t k = 0; k < 2; k++) {
232 fCurrentCoef[k] = 0.0;
233 fCurrentCoef2[k] = 0.0;
234 }
235 for (Int_t i = 0; i < 3; i++) {
236 fPhd[i] = 0.0;
237 fDet[i] = 0;
238 }
239 if(c.fCalDet) fCalDet = new AliTRDCalDet(*c.fCalDet);
240 if(c.fCalDet2) fCalDet2 = new AliTRDCalDet(*c.fCalDet2);
241
242 if(c.fCalROC) fCalROC = new AliTRDCalROC(*c.fCalROC);
243 if(c.fCalROC2) fCalROC = new AliTRDCalROC(*c.fCalROC2);
244
245 fVectorFit.SetName(c.fVectorFit.GetName());
246 for(Int_t k = 0; k < c.fVectorFit.GetEntriesFast(); k++){
247 AliTRDFitInfo *fitInfo = new AliTRDFitInfo();
248 Int_t detector = ((AliTRDFitInfo *)c.fVectorFit.UncheckedAt(k))->GetDetector();
249 Int_t ntotal = 1;
250 if (GetChamber(detector) == 2) {
251 ntotal = 1728;
252 }
253 else {
254 ntotal = 2304;
255 }
256 Float_t *coef = new Float_t[ntotal];
257 for (Int_t i = 0; i < ntotal; i++) {
258 coef[i] = ((AliTRDFitInfo *)c.fVectorFit.UncheckedAt(k))->GetCoef()[i];
259 }
260 fitInfo->SetCoef(coef);
261 fitInfo->SetDetector(detector);
262 fVectorFit.Add((TObject *) fitInfo);
263 }
264 fVectorFit.SetName(c.fVectorFit.GetName());
265 for(Int_t k = 0; k < c.fVectorFit2.GetEntriesFast(); k++){
266 AliTRDFitInfo *fitInfo = new AliTRDFitInfo();
267 Int_t detector = ((AliTRDFitInfo *)c.fVectorFit2.UncheckedAt(k))->GetDetector();
268 Int_t ntotal = 1;
269 if (GetChamber(detector) == 2) {
270 ntotal = 1728;
271 }
272 else {
273 ntotal = 2304;
274 }
275 Float_t *coef = new Float_t[ntotal];
276 for (Int_t i = 0; i < ntotal; i++) {
277 coef[i] = ((AliTRDFitInfo *)c.fVectorFit2.UncheckedAt(k))->GetCoef()[i];
278 }
279 fitInfo->SetCoef(coef);
280 fitInfo->SetDetector(detector);
281 fVectorFit2.Add((TObject *) fitInfo);
282 }
283 if (fGeo) {
284 delete fGeo;
285 }
286 fGeo = new AliTRDgeometry();
287
288}
289//____________________________________________________________________________________
290AliTRDCalibraFit::~AliTRDCalibraFit()
291{
292 //
293 // AliTRDCalibraFit destructor
294 //
295 if ( fDebugStreamer ) delete fDebugStreamer;
296 if ( fCalDet ) delete fCalDet;
297 if ( fCalDet2 ) delete fCalDet2;
298 if ( fCalROC ) delete fCalROC;
299 if ( fCalROC2 ) delete fCalROC2;
300 fVectorFit.Delete();
301 fVectorFit2.Delete();
302 if (fGeo) {
303 delete fGeo;
304 }
305
306}
307//_____________________________________________________________________________
308void AliTRDCalibraFit::Destroy()
309{
310 //
311 // Delete instance
312 //
313
314 if (fgInstance) {
315 delete fgInstance;
316 fgInstance = 0x0;
317 }
318
319}
320//__________________________________________________________________________________
321void AliTRDCalibraFit::RangeChargeIntegration(Float_t vdrift, Float_t t0, Int_t &begin, Int_t &peak, Int_t &end)
322{
323 //
324 // From the drift velocity and t0
325 // return the position of the peak and maximum negative slope
326 //
327
328 const Float_t kDrWidth = AliTRDgeometry::DrThick(); // drift region
329 Double_t widbins = 0.1; // 0.1 mus
330
331 //peak and maxnegslope in mus
332 Double_t begind = t0*widbins + fT0Shift0;
333 Double_t peakd = t0*widbins + fT0Shift1;
334 Double_t maxnegslope = (kDrWidth + vdrift*peakd)/vdrift;
335
336 // peak and maxnegslope in timebin
337 begin = TMath::Nint(begind*widbins);
338 peak = TMath::Nint(peakd*widbins);
339 end = TMath::Nint(maxnegslope*widbins);
340
341}
342//____________Functions fit Online CH2d________________________________________
343Bool_t AliTRDCalibraFit::AnalyseCH(TH2I *ch)
344{
345 //
346 // Fit the 1D histos, projections of the 2D ch on the Xaxis, for each
347 // calibration group normalized the resulted coefficients (to 1 normally)
348 //
349
350 // Set the calibration mode
351 const char *name = ch->GetTitle();
352 SetModeCalibration(name,0);
353
354 // Number of Ybins (detectors or groups of pads)
355 Int_t nbins = ch->GetNbinsX();// charge
356 Int_t nybins = ch->GetNbinsY();// groups number
357 if (!InitFit(nybins,0)) {
358 return kFALSE;
359 }
360 if (!InitFitCH()) {
361 return kFALSE;
362 }
363 fStatisticMean = 0.0;
364 fNumberFit = 0;
365 fNumberFitSuccess = 0;
366 fNumberEnt = 0;
367 // Init fCountDet and fCount
368 InitfCountDetAndfCount(0);
369 // Beginning of the loop betwwen dect1 and dect2
370 for (Int_t idect = fDect1; idect < fDect2; idect++) {
371 // Determination of fNnZ, fNnRphi, fNfragZ and fNfragRphi...
372 UpdatefCountDetAndfCount(idect,0);
373 ReconstructFitRowMinRowMax(idect, 0);
374 // Take the histo
375 TH1I *projch = (TH1I *) ch->ProjectionX("projch",idect+1,idect+1,(Option_t *)"e");
376 projch->SetDirectory(0);
377 // Number of entries for this calibration group
378 Double_t nentries = 0.0;
379 Double_t mean = 0.0;
380 for (Int_t k = 0; k < nbins; k++) {
381 Int_t binnb = (nbins+2)*(idect+1)+(k+1);
382 nentries += ch->GetBinContent(binnb);
383 mean += projch->GetBinCenter(k+1)*projch->GetBinContent(k+1);
384 projch->SetBinError(k+1,TMath::Sqrt(projch->GetBinContent(k+1)));
385 }
386 projch->SetEntries(nentries);
387 //printf("The number of entries for the group %d is %f\n",idect,nentries);
388 if (nentries > 0) {
389 fNumberEnt++;
390 mean /= nentries;
391 }
392 // Rebin and statistic stuff
393 if (fRebin > 1) {
394 projch = ReBin((TH1I *) projch);
395 }
396 // This detector has not enough statistics or was off
397 if (nentries <= fMinEntries) {
398 NotEnoughStatisticCH(idect);
399 if (fDebugLevel != 1) {
400 delete projch;
401 }
402 continue;
403 }
404 // Statistics of the group fitted
405 fStatisticMean += nentries;
406 fNumberFit++;
407 //Method choosen
408 switch(fMethod)
409 {
410 case 0: FitMeanW((TH1 *) projch, nentries); break;
411 case 1: FitMean((TH1 *) projch, nentries, mean); break;
412 case 2: FitCH((TH1 *) projch, mean); break;
413 case 3: FitBisCH((TH1 *) projch, mean); break;
414 default: return kFALSE;
415 }
416 // Fill Infos Fit
417 FillInfosFitCH(idect);
418 // Memory!!!
419 if (fDebugLevel != 1) {
420 delete projch;
421 }
422 } // Boucle object
423 // Normierungcharge
424 if (fDebugLevel != 1) {
425 NormierungCharge();
426 }
427 // Mean Statistic
428 if (fNumberFit > 0) {
429 AliInfo(Form("There are %d with at least one entries. %d fits have been proceeded (sucessfully or not...). There is a mean statistic of: %d over these fitted histograms and %d successfulled fits",fNumberEnt, fNumberFit, (Int_t) fStatisticMean/fNumberFit, fNumberFitSuccess));
430 fStatisticMean = fStatisticMean / fNumberFit;
431 }
432 else {
433 AliInfo(Form("There are %d with at least one entries. There is no fit!",fNumberEnt));
434 }
435 delete fDebugStreamer;
436 fDebugStreamer = 0x0;
437
438 return kTRUE;
439}
440//____________Functions fit Online CH2d________________________________________
441Bool_t AliTRDCalibraFit::AnalyseCH(AliTRDCalibraVector *calvect)
442{
443 //
444 // Reconstruct a 1D histo from the vectorCH for each calibration group,
445 // fit the histo, normalized the resulted coefficients (to 1 normally)
446 //
447
448 // Set the calibraMode
449 const char *name = calvect->GetNameCH();
450 SetModeCalibration(name,0);
451
452 // Number of Xbins (detectors or groups of pads)
453 if (!InitFit((432*calvect->GetDetCha0(0)+108*calvect->GetDetCha2(0)),0)) {
454 return kFALSE;
455 }
456 if (!InitFitCH()) {
457 return kFALSE;
458 }
459 fStatisticMean = 0.0;
460 fNumberFit = 0;
461 fNumberFitSuccess = 0;
462 fNumberEnt = 0;
463 // Init fCountDet and fCount
464 InitfCountDetAndfCount(0);
465 // Beginning of the loop between dect1 and dect2
466 for (Int_t idect = fDect1; idect < fDect2; idect++) {
467 // Determination of fNnZ, fNnRphi, fNfragZ and fNfragRphi...........
468 UpdatefCountDetAndfCount(idect,0);
469 ReconstructFitRowMinRowMax(idect,0);
470 // Take the histo
471 Double_t nentries = 0.0;
472 Double_t mean = 0.0;
473 TH1F *projch = 0x0;
474 Bool_t something = kTRUE;
475 if(!calvect->GetCHEntries(fCountDet)) something = kFALSE;
476 if(something){
477 TString name("CH");
478 name += idect;
479 projch = calvect->ConvertVectorCHHisto(fCountDet,(idect-(fCount-(fCalibraMode->GetNfragZ(0)*fCalibraMode->GetNfragRphi(0)))),(const char *) name);
480 projch->SetDirectory(0);
481 for (Int_t k = 0; k < calvect->GetNumberBinCharge(); k++) {
482 nentries += projch->GetBinContent(k+1);
483 mean += projch->GetBinCenter(k+1)*projch->GetBinContent(k+1);
484 }
485 if (nentries > 0) {
486 fNumberEnt++;
487 mean /= nentries;
488 }
489 //printf("The number of entries for the group %d is %f\n",idect,nentries);
490 // Rebin
491 if (fRebin > 1) {
492 projch = ReBin((TH1F *) projch);
493 }
494 }
495 // This detector has not enough statistics or was not found in VectorCH
496 if (nentries <= fMinEntries) {
497 NotEnoughStatisticCH(idect);
498 if (fDebugLevel != 1) {
499 if(projch) delete projch;
500 }
501 continue;
502 }
503 // Statistic of the histos fitted
504 fStatisticMean += nentries;
505 fNumberFit++;
506 //Method choosen
507 switch(fMethod)
508 {
509 case 0: FitMeanW((TH1 *) projch, nentries); break;
510 case 1: FitMean((TH1 *) projch, nentries, mean); break;
511 case 2: FitCH((TH1 *) projch, mean); break;
512 case 3: FitBisCH((TH1 *) projch, mean); break;
513 default: return kFALSE;
514 }
515 // Fill Infos Fit
516 FillInfosFitCH(idect);
517 // Memory!!!
518 if (fDebugLevel != 1) {
519 delete projch;
520 }
521 } // Boucle object
522 // Normierungcharge
523 if (fDebugLevel != 1) {
524 NormierungCharge();
525 }
526 // Mean Statistics
527 if (fNumberFit > 0) {
528 AliInfo(Form("There are %d with at least one entries. %d fits have been proceeded (sucessfully or not...). There is a mean statistic of: %d over these fitted histograms and %d successfulled fits",fNumberEnt, fNumberFit, (Int_t) fStatisticMean/fNumberFit, fNumberFitSuccess));
529 fStatisticMean = fStatisticMean / fNumberFit;
530 }
531 else {
532 AliInfo(Form("There are %d with at least one entries. There is no fit!",fNumberEnt));
533 }
534 delete fDebugStreamer;
535 fDebugStreamer = 0x0;
536 return kTRUE;
537}
538//________________functions fit Online PH2d____________________________________
539Bool_t AliTRDCalibraFit::AnalysePH(TProfile2D *ph)
540{
541 //
542 // Take the 1D profiles (average pulse height), projections of the 2D PH
543 // on the Xaxis, for each calibration group
544 // Reconstruct a drift velocity
545 // A first calibration of T0 is also made using the same method
546 //
547
548 // Set the calibration mode
549 const char *name = ph->GetTitle();
550 SetModeCalibration(name,1);
551
552 // Number of Xbins (detectors or groups of pads)
553 Int_t nbins = ph->GetNbinsX();// time
554 Int_t nybins = ph->GetNbinsY();// calibration group
555 if (!InitFit(nybins,1)) {
556 return kFALSE;
557 }
558 if (!InitFitPH()) {
559 return kFALSE;
560 }
561 fStatisticMean = 0.0;
562 fNumberFit = 0;
563 fNumberFitSuccess = 0;
564 fNumberEnt = 0;
565 // Init fCountDet and fCount
566 InitfCountDetAndfCount(1);
567 // Beginning of the loop
568 for (Int_t idect = fDect1; idect < fDect2; idect++) {
569 // Determination of fNnZ, fNnRphi, fNfragZ and fNfragRphi.......
570 UpdatefCountDetAndfCount(idect,1);
571 ReconstructFitRowMinRowMax(idect,1);
572 // Take the histo
573 TH1D *projph = (TH1D *) ph->ProjectionX("projph",idect+1,idect+1,(Option_t *) "e");
574 projph->SetDirectory(0);
575 // Number of entries for this calibration group
576 Double_t nentries = 0;
577 for (Int_t k = 0; k < nbins; k++) {
578 Int_t binnb = (nbins+2)*(idect+1)+(k+1);
579 nentries += ph->GetBinEntries(binnb);
580 }
581 if (nentries > 0) {
582 fNumberEnt++;
583 }
584 //printf("The number of entries for the group %d is %f\n",idect,nentries);
585 // This detector has not enough statistics or was off
586 if (nentries <= fMinEntries) {
587 //printf("Not enough statistic!\n");
588 NotEnoughStatisticPH(idect);
589 if (fDebugLevel != 1) {
590 delete projph;
591 }
592 continue;
593 }
594 // Statistics of the histos fitted
595 fNumberFit++;
596 fStatisticMean += nentries;
597 // Calcul of "real" coef
598 CalculVdriftCoefMean();
599 CalculT0CoefMean();
600 //Method choosen
601 switch(fMethod)
602 {
603 case 0: FitLagrangePoly((TH1 *) projph); break;
604 case 1: FitPente((TH1 *) projph); break;
605 case 2: FitPH((TH1 *) projph,(Int_t) (idect - fDect1)); break;
606 default: return kFALSE;
607 }
608 // Fill the tree if end of a detector or only the pointer to the branch!!!
609 FillInfosFitPH(idect);
610 // Memory!!!
611 if (fDebugLevel != 1) {
612 delete projph;
613 }
614 } // Boucle object
615 // Mean Statistic
616 if (fNumberFit > 0) {
617 AliInfo(Form("There are %d with at least one entries. %d fits have been proceeded (sucessfully or not...). There is a mean statistic of: %d over these fitted histograms and %d successfulled fits",fNumberEnt, fNumberFit, (Int_t) fStatisticMean/fNumberFit,fNumberFitSuccess));
618 fStatisticMean = fStatisticMean / fNumberFit;
619 }
620 else {
621 AliInfo(Form("There are %d with at least one entries. There is no fit!",fNumberEnt));
622 }
623 delete fDebugStreamer;
624 fDebugStreamer = 0x0;
625 return kTRUE;
626}
627//____________Functions fit Online PH2d________________________________________
628Bool_t AliTRDCalibraFit::AnalysePH(AliTRDCalibraVector *calvect)
629{
630 //
631 // Reconstruct the average pulse height from the vectorPH for each
632 // calibration group
633 // Reconstruct a drift velocity
634 // A first calibration of T0 is also made using the same method (slope method)
635 //
636
637 // Set the calibration mode
638 const char *name = calvect->GetNamePH();
639 SetModeCalibration(name,1);
640
641 // Number of Xbins (detectors or groups of pads)
642 if (!InitFit((432*calvect->GetDetCha0(1)+108*calvect->GetDetCha2(1)),1)) {
643 return kFALSE;
644 }
645 if (!InitFitPH()) {
646 return kFALSE;
647 }
648 fStatisticMean = 0.0;
649 fNumberFit = 0;
650 fNumberFitSuccess = 0;
651 fNumberEnt = 0;
652 // Init fCountDet and fCount
653 InitfCountDetAndfCount(1);
654 // Beginning of the loop
655 for (Int_t idect = fDect1; idect < fDect2; idect++) {
656 // Determination of fNnZ, fNnRphi, fNfragZ and fNfragRphi...........
657 UpdatefCountDetAndfCount(idect,1);
658 ReconstructFitRowMinRowMax(idect,1);
659 // Take the histo
660 TH1F *projph = 0x0;
661 fEntriesCurrent = 0;
662 Bool_t something = kTRUE;
663 if(!calvect->GetPHEntries(fCountDet)) something = kFALSE;
664 if(something){
665 TString name("PH");
666 name += idect;
667 projph = CorrectTheError((TGraphErrors *) (calvect->ConvertVectorPHTGraphErrors(fCountDet,(idect-(fCount-(fCalibraMode->GetNfragZ(1)*fCalibraMode->GetNfragRphi(1)))),(const char *) name)));
668 projph->SetDirectory(0);
669 }
670 //printf("The number of entries for the group %d is %d\n",idect,fEntriesCurrent);
671 // This detector has not enough statistics or was off
672 if (fEntriesCurrent <= fMinEntries) {
673 //printf("Not enough stat!\n");
674 NotEnoughStatisticPH(idect);
675 if (fDebugLevel != 1) {
676 if(projph) delete projph;
677 }
678 continue;
679 }
680 // Statistic of the histos fitted
681 fNumberFit++;
682 fStatisticMean += fEntriesCurrent;
683 // Calcul of "real" coef
684 CalculVdriftCoefMean();
685 CalculT0CoefMean();
686 //Method choosen
687 switch(fMethod)
688 {
689 case 0: FitLagrangePoly((TH1 *) projph); break;
690 case 1: FitPente((TH1 *) projph); break;
691 case 2: FitPH((TH1 *) projph,(Int_t) (idect - fDect1)); break;
692 default: return kFALSE;
693 }
694 // Fill the tree if end of a detector or only the pointer to the branch!!!
695 FillInfosFitPH(idect);
696 // Memory!!!
697 if (fDebugLevel != 1) {
698 delete projph;
699 }
700 } // Boucle object
701
702 // Mean Statistic
703 if (fNumberFit > 0) {
704 AliInfo(Form("There are %d with at least one entries. %d fits have been proceeded (sucessfully or not...). There is a mean statistic of: %d over these fitted histograms and %d successfulled fits",fNumberEnt, fNumberFit, (Int_t) fStatisticMean/fNumberFit,fNumberFitSuccess));
705 fStatisticMean = fStatisticMean / fNumberFit;
706 }
707 else {
708 AliInfo(Form("There are %d with at least one entries. There is no fit!",fNumberEnt));
709 }
710 delete fDebugStreamer;
711 fDebugStreamer = 0x0;
712 return kTRUE;
713}
714//____________Functions fit Online PRF2d_______________________________________
715Bool_t AliTRDCalibraFit::AnalysePRF(TProfile2D *prf)
716{
717 //
718 // Take the 1D profiles (pad response function), projections of the 2D PRF
719 // on the Xaxis, for each calibration group
720 // Fit with a gaussian to reconstruct the sigma of the pad response function
721 //
722
723 // Set the calibration mode
724 const char *name = prf->GetTitle();
725 SetModeCalibration(name,2);
726
727 // Number of Ybins (detectors or groups of pads)
728 Int_t nybins = prf->GetNbinsY();// calibration groups
729 Int_t nbins = prf->GetNbinsX();// bins
730 Int_t nbg = GetNumberOfGroupsPRF((const char *)prf->GetTitle());
731 if((nbg > 0) || (nbg == -1)) return kFALSE;
732 if (!InitFit(nybins,2)) {
733 return kFALSE;
734 }
735 if (!InitFitPRF()) {
736 return kFALSE;
737 }
738 fStatisticMean = 0.0;
739 fNumberFit = 0;
740 fNumberFitSuccess = 0;
741 fNumberEnt = 0;
742 // Init fCountDet and fCount
743 InitfCountDetAndfCount(2);
744 // Beginning of the loop
745 for (Int_t idect = fDect1; idect < fDect2; idect++) {
746 // Determination of fNnZ, fNnRphi, fNfragZ and fNfragRphi......
747 UpdatefCountDetAndfCount(idect,2);
748 ReconstructFitRowMinRowMax(idect,2);
749 // Take the histo
750 TH1D *projprf = (TH1D *) prf->ProjectionX("projprf",idect+1,idect+1,(Option_t *) "e");
751 projprf->SetDirectory(0);
752 // Number of entries for this calibration group
753 Double_t nentries = 0;
754 for (Int_t k = 0; k < nbins; k++) {
755 Int_t binnb = (nbins+2)*(idect+1)+(k+1);
756 nentries += prf->GetBinEntries(binnb);
757 }
758 if(nentries > 0) fNumberEnt++;
759 // This detector has not enough statistics or was off
760 if (nentries <= fMinEntries) {
761 NotEnoughStatisticPRF(idect);
762 if (fDebugLevel != 1) {
763 delete projprf;
764 }
765 continue;
766 }
767 // Statistics of the histos fitted
768 fNumberFit++;
769 fStatisticMean += nentries;
770 // Calcul of "real" coef
771 CalculPRFCoefMean();
772 //Method choosen
773 switch(fMethod)
774 {
775 case 0: FitPRF((TH1 *) projprf); break;
776 case 1: RmsPRF((TH1 *) projprf); break;
777 default: return kFALSE;
778 }
779 // Fill the tree if end of a detector or only the pointer to the branch!!!
780 FillInfosFitPRF(idect);
781 // Memory!!!
782 if (fDebugLevel != 1) {
783 delete projprf;
784 }
785 } // Boucle object
786 // Mean Statistic
787 if (fNumberFit > 0) {
788 AliInfo(Form("There are %d with at least one entries.",fNumberEnt));
789 AliInfo(Form("%d fits have been proceeded (sucessfully or not...).",fNumberFit));
790 AliInfo(Form("There is a mean statistic of: %d over these fitted histograms and %d successfulled fits"
791 ,(Int_t) fStatisticMean/fNumberFit,fNumberFitSuccess));
792 fStatisticMean = fStatisticMean / fNumberFit;
793 }
794 else {
795 AliInfo(Form("There are %d with at least one entries. There is no fit!",fNumberEnt));
796 }
797 delete fDebugStreamer;
798 fDebugStreamer = 0x0;
799 return kTRUE;
800}
801//____________Functions fit Online PRF2d_______________________________________
802Bool_t AliTRDCalibraFit::AnalysePRFMarianFit(TProfile2D *prf)
803{
804 //
805 // Take the 1D profiles (pad response function), projections of the 2D PRF
806 // on the Xaxis, for each calibration group
807 // Fit with a gaussian to reconstruct the sigma of the pad response function
808 //
809
810 // Set the calibration mode
811 const char *name = prf->GetTitle();
812 SetModeCalibration(name,2);
813
814 // Number of Ybins (detectors or groups of pads)
815 TAxis *xprf = prf->GetXaxis();
816 TAxis *yprf = prf->GetYaxis();
817 Int_t nybins = yprf->GetNbins();// calibration groups
818 Int_t nbins = xprf->GetNbins();// bins
819 Float_t lowedge = (Float_t) xprf->GetBinLowEdge(1);//lowedge in bins
820 Float_t upedge = (Float_t) xprf->GetBinUpEdge(nbins);//upedge in bins
821 Int_t nbg = GetNumberOfGroupsPRF((const char *)name);
822 if(nbg == -1) return kFALSE;
823 if(nbg > 0) fMethod = 1;
824 else fMethod = 0;
825 if (!InitFit(nybins,2)) {
826 return kFALSE;
827 }
828 if (!InitFitPRF()) {
829 return kFALSE;
830 }
831 fStatisticMean = 0.0;
832 fNumberFit = 0;
833 fNumberFitSuccess = 0;
834 fNumberEnt = 0;
835 // Init fCountDet and fCount
836 InitfCountDetAndfCount(2);
837 // Beginning of the loop
838 for (Int_t idect = fDect1; idect < fDect2; idect++) {
839 // Determination of fNnZ, fNnRphi, fNfragZ and fNfragRphi.......
840 UpdatefCountDetAndfCount(idect,2);
841 ReconstructFitRowMinRowMax(idect,2);
842 // Build the array of entries and sum
843 TArrayD arraye = TArrayD(nbins);
844 TArrayD arraym = TArrayD(nbins);
845 TArrayD arrayme = TArrayD(nbins);
846 Double_t nentries = 0;
847 //printf("nbins %d\n",nbins);
848 for (Int_t k = 0; k < nbins; k++) {
849 Int_t binnb = (nbins+2)*(idect+1)+(k+1);
850 Double_t entries = (Double_t)prf->GetBinEntries(binnb);
851 Double_t mean = (Double_t)prf->GetBinContent(binnb);
852 Double_t error = (Double_t)prf->GetBinError(binnb);
853 //printf("for %d we have %f\n",k,entries);
854 nentries += entries;
855 arraye.AddAt(entries,k);
856 arraym.AddAt(mean,k);
857 arrayme.AddAt(error,k);
858 }
859 if(nentries > 0) fNumberEnt++;
860 //printf("The number of entries for the group %d is %f\n",idect,nentries);
861 // This detector has not enough statistics or was off
862 if (nentries <= fMinEntries) {
863 NotEnoughStatisticPRF(idect);
864 continue;
865 }
866 // Statistics of the histos fitted
867 fNumberFit++;
868 fStatisticMean += nentries;
869 // Calcul of "real" coef
870 CalculPRFCoefMean();
871 //Method choosen
872 switch(fMethod)
873 {
874 case 0: FitPRFGausMI( arraye.GetArray(), arraym.GetArray(), arrayme.GetArray(), nbins, lowedge, upedge); break;
875 case 1: FitTnpRange( arraye.GetArray(), arraym.GetArray(), arrayme.GetArray(), nbg, nbins); break;
876 default: return kFALSE;
877 }
878 // Fill the tree if end of a detector or only the pointer to the branch!!!
879 FillInfosFitPRF(idect);
880 } // Boucle object
881 // Mean Statistic
882 if (fNumberFit > 0) {
883 AliInfo(Form("There are %d with at least one entries.",fNumberEnt));
884 AliInfo(Form("%d fits have been proceeded (sucessfully or not...).",fNumberFit));
885 AliInfo(Form("There is a mean statistic of: %d over these fitted histograms and %d successfulled fits"
886 ,(Int_t) fStatisticMean/fNumberFit,fNumberFitSuccess));
887 fStatisticMean = fStatisticMean / fNumberFit;
888 }
889 else {
890 AliInfo(Form("There are %d with at least one entries. There is no fit!",fNumberEnt));
891 }
892 delete fDebugStreamer;
893 fDebugStreamer = 0x0;
894 return kTRUE;
895}
896//____________Functions fit Online PRF2d_______________________________________
897Bool_t AliTRDCalibraFit::AnalysePRF(AliTRDCalibraVector *calvect)
898{
899 //
900 // Reconstruct the 1D histo (pad response function) from the vectorPRD for
901 // each calibration group
902 // Fit with a gaussian to reconstruct the sigma of the pad response function
903 //
904
905 // Set the calibra mode
906 const char *name = calvect->GetNamePRF();
907 SetModeCalibration(name,2);
908 //printf("test0 %s\n",name);
909
910 // Number of Xbins (detectors or groups of pads)
911 if (!InitFit((432*calvect->GetDetCha0(2)+108*calvect->GetDetCha2(2)),2)) {
912 //printf("test1\n");
913 return kFALSE;
914 }
915 if (!InitFitPRF()) {
916 ///printf("test2\n");
917 return kFALSE;
918 }
919 fStatisticMean = 0.0;
920 fNumberFit = 0;
921 fNumberFitSuccess = 0;
922 fNumberEnt = 0;
923 // Init fCountDet and fCount
924 InitfCountDetAndfCount(2);
925 // Beginning of the loop
926 for (Int_t idect = fDect1; idect < fDect2; idect++) {
927 // Determination of fNnZ, fNnRphi, fNfragZ and fNfragRphi........
928 UpdatefCountDetAndfCount(idect,2);
929 ReconstructFitRowMinRowMax(idect,2);
930 // Take the histo
931 TH1F *projprf = 0x0;
932 fEntriesCurrent = 0;
933 Bool_t something = kTRUE;
934 if(!calvect->GetPRFEntries(fCountDet)) something = kFALSE;
935 if(something){
936 TString name("PRF");
937 name += idect;
938 projprf = CorrectTheError((TGraphErrors *) (calvect->ConvertVectorPRFTGraphErrors(fCountDet,(idect-(fCount-(fCalibraMode->GetNfragZ(1)*fCalibraMode->GetNfragRphi(1)))),(const char *) name)));
939 projprf->SetDirectory(0);
940 }
941 // This detector has not enough statistics or was off
942 if (fEntriesCurrent <= fMinEntries) {
943 NotEnoughStatisticPRF(idect);
944 if (fDebugLevel != 1) {
945 if(projprf) delete projprf;
946 }
947 continue;
948 }
949 // Statistic of the histos fitted
950 fNumberFit++;
951 fStatisticMean += fEntriesCurrent;
952 // Calcul of "real" coef
953 CalculPRFCoefMean();
954 //Method choosen
955 switch(fMethod)
956 {
957 case 1: FitPRF((TH1 *) projprf); break;
958 case 2: RmsPRF((TH1 *) projprf); break;
959 default: return kFALSE;
960 }
961 // Fill the tree if end of a detector or only the pointer to the branch!!!
962 FillInfosFitPRF(idect);
963 // Memory!!!
964 if (fDebugLevel != 1) {
965 delete projprf;
966 }
967 } // Boucle object
968 // Mean Statistics
969 if (fNumberFit > 0) {
970 AliInfo(Form("There are %d with at least one entries. %d fits have been proceeded (sucessfully or not...). There is a mean statistic of: %d over these fitted histograms and %d successfulled fits",fNumberEnt, fNumberFit, (Int_t) fStatisticMean/fNumberFit,fNumberFitSuccess));
971 }
972 else {
973 AliInfo(Form("There are %d with at least one entries. There is no fit!",fNumberEnt));
974 }
975 delete fDebugStreamer;
976 fDebugStreamer = 0x0;
977 return kTRUE;
978}
979//____________Functions fit Online PRF2d_______________________________________
980Bool_t AliTRDCalibraFit::AnalysePRFMarianFit(AliTRDCalibraVector *calvect)
981{
982 //
983 // Reconstruct the 1D histo (pad response function) from the vectorPRD for
984 // each calibration group
985 // Fit with a gaussian to reconstruct the sigma of the pad response function
986 //
987
988 // Set the calibra mode
989 const char *name = calvect->GetNamePRF();
990 SetModeCalibration(name,2);
991 //printf("test0 %s\n",name);
992 Int_t nbg = GetNumberOfGroupsPRF((const char *)name);
993 printf("test1 %d\n",nbg);
994 if(nbg == -1) return kFALSE;
995 if(nbg > 0) fMethod = 1;
996 else fMethod = 0;
997 // Number of Xbins (detectors or groups of pads)
998 if (!InitFit((432*calvect->GetDetCha0(2)+108*calvect->GetDetCha2(2)),2)) {
999 //printf("test2\n");
1000 return kFALSE;
1001 }
1002 if (!InitFitPRF()) {
1003 //printf("test3\n");
1004 return kFALSE;
1005 }
1006 fStatisticMean = 0.0;
1007 fNumberFit = 0;
1008 fNumberFitSuccess = 0;
1009 fNumberEnt = 0;
1010 // Variables
1011 Int_t nbins = 0;
1012 Double_t *arrayx = 0;
1013 Double_t *arraye = 0;
1014 Double_t *arraym = 0;
1015 Double_t *arrayme = 0;
1016 Float_t lowedge = 0.0;
1017 Float_t upedge = 0.0;
1018 // Init fCountDet and fCount
1019 InitfCountDetAndfCount(2);
1020 // Beginning of the loop
1021 for (Int_t idect = fDect1; idect < fDect2; idect++) {
1022 // Determination of fNnZ, fNnRphi, fNfragZ and fNfragRphi......
1023 UpdatefCountDetAndfCount(idect,2);
1024 ReconstructFitRowMinRowMax(idect,2);
1025 // Take the histo
1026 TGraphErrors *projprftree = 0x0;
1027 fEntriesCurrent = 0;
1028 Bool_t something = kTRUE;
1029 if(!calvect->GetPRFEntries(fCountDet)) something = kFALSE;
1030 if(something){
1031 TString name("PRF");
1032 name += idect;
1033 projprftree = calvect->ConvertVectorPRFTGraphErrors(fCountDet,(idect-(fCount-(fCalibraMode->GetNfragZ(1)*fCalibraMode->GetNfragRphi(1)))),(const char *) name);
1034 nbins = projprftree->GetN();
1035 arrayx = (Double_t *)projprftree->GetX();
1036 arraye = (Double_t *)projprftree->GetEX();
1037 arraym = (Double_t *)projprftree->GetY();
1038 arrayme = (Double_t *)projprftree->GetEY();
1039 Float_t step = arrayx[1]-arrayx[0];
1040 lowedge = arrayx[0] - step/2.0;
1041 upedge = arrayx[(nbins-1)] + step/2.0;
1042 //printf("nbins est %d\n",nbins);
1043 for(Int_t k = 0; k < nbins; k++){
1044 fEntriesCurrent += (Int_t)arraye[k];
1045 //printf("for %d we have %f, %f\n",k,arraye[k],((projprftree->GetEX())[k]));
1046 if(arraye[k]>0.0) arrayme[k] = TMath::Sqrt(TMath::Abs(arrayme[k]-arraym[k]*arraym[k])/arraye[k]);
1047 }
1048 if(fEntriesCurrent > 0) fNumberEnt++;
1049 }
1050 //printf("The number of entries for the group %d is %d\n",idect,fEntriesCurrent);
1051 // This detector has not enough statistics or was off
1052 if (fEntriesCurrent <= fMinEntries) {
1053 NotEnoughStatisticPRF(idect);
1054 if(projprftree) delete projprftree;
1055 continue;
1056 }
1057 // Statistic of the histos fitted
1058 fNumberFit++;
1059 fStatisticMean += fEntriesCurrent;
1060 // Calcul of "real" coef
1061 CalculPRFCoefMean();
1062 //Method choosen
1063 switch(fMethod)
1064 {
1065 case 0: FitPRFGausMI(arraye,arraym,arrayme,nbins,lowedge,upedge); break;
1066 case 1: FitTnpRange(arraye,arraym,arrayme,nbg,nbins); break;
1067 default: return kFALSE;
1068 }
1069 // Fill the tree if end of a detector or only the pointer to the branch!!!
1070 FillInfosFitPRF(idect);
1071 // Memory!!!
1072 if (fDebugLevel != 1) {
1073 delete projprftree;
1074 }
1075 } // Boucle object
1076 // Mean Statistics
1077 if (fNumberFit > 0) {
1078 AliInfo(Form("There are %d with at least one entries. %d fits have been proceeded (sucessfully or not...). There is a mean statistic of: %d over these fitted histograms and %d successfulled fits",fNumberEnt, fNumberFit, (Int_t) fStatisticMean/fNumberFit,fNumberFitSuccess));
1079 }
1080 else {
1081 AliInfo(Form("There are %d with at least one entries. There is no fit!",fNumberEnt));
1082 }
1083 delete fDebugStreamer;
1084 fDebugStreamer = 0x0;
1085 return kTRUE;
1086}
1087//____________Functions fit Online CH2d________________________________________
1088Bool_t AliTRDCalibraFit::AnalyseLinearFitters(AliTRDCalibraVdriftLinearFit *calivdli)
1089{
1090 //
1091 // The linear method
1092 //
1093
1094 fStatisticMean = 0.0;
1095 fNumberFit = 0;
1096 fNumberFitSuccess = 0;
1097 fNumberEnt = 0;
1098 if(!InitFitLinearFitter()) return kFALSE;
1099
1100
1101 for(Int_t idet = 0; idet < 540; idet++){
1102
1103
1104 //printf("detector number %d\n",idet);
1105
1106 // Take the result
1107 TVectorD param(2);
1108 TVectorD error(3);
1109 fEntriesCurrent = 0;
1110 fCountDet = idet;
1111 Bool_t here = calivdli->GetParam(idet,&param);
1112 Bool_t heree = calivdli->GetError(idet,&error);
1113 //printf("here %d and heree %d\n",here, heree);
1114 if(heree) {
1115 fEntriesCurrent = (Int_t) error[2];
1116 fNumberEnt++;
1117 }
1118 //printf("Number of entries %d\n",fEntriesCurrent);
1119 // Nothing found or not enough statistic
1120 if((!heree) || (!here) || (fEntriesCurrent <= fMinEntries)) {
1121 NotEnoughStatisticLinearFitter();
1122 continue;
1123 }
1124 //param.Print();
1125 //error.Print();
1126 //Statistics
1127 fNumberFit++;
1128 fStatisticMean += fEntriesCurrent;
1129
1130 // Check the fit
1131 if((-(param[1])) <= 0.0) {
1132 NotEnoughStatisticLinearFitter();
1133 continue;
1134 }
1135
1136 // CalculDatabaseVdriftandTan
1137 CalculVdriftLorentzCoef();
1138
1139 // Statistics
1140 fNumberFitSuccess ++;
1141
1142 // Put the fCurrentCoef
1143 fCurrentCoef[0] = -param[1];
1144 // here the database must be the one of the reconstruction for the lorentz angle....
1145 fCurrentCoef2[0] = (param[0]+fCurrentCoef[1]*fCurrentCoef2[1])/fCurrentCoef[0];
1146 fCurrentCoefE = error[1];
1147 fCurrentCoefE2 = error[0];
1148 if((fCurrentCoef2[0] != 0.0) && (param[0] != 0.0)){
1149 fCurrentCoefE2 = (fCurrentCoefE2/param[0]+fCurrentCoefE/fCurrentCoef[0])*fCurrentCoef2[0];
1150 }
1151
1152 // Fill
1153 FillInfosFitLinearFitter();
1154
1155
1156 }
1157 // Mean Statistics
1158 if (fNumberFit > 0) {
1159 AliInfo(Form("There are %d with at least one entries. %d fits have been proceeded (sucessfully or not...). There is a mean statistic of: %d over these fitted histograms and %d successfulled fits",fNumberEnt, fNumberFit, (Int_t) fStatisticMean/fNumberFit,fNumberFitSuccess));
1160 }
1161 else {
1162 AliInfo(Form("There are %d with at least one entries. There is no fit!",fNumberEnt));
1163 }
1164 delete fDebugStreamer;
1165 fDebugStreamer = 0x0;
1166 return kTRUE;
1167
1168}
1169//____________Functions for seeing if the pad is really okey___________________
1170//_____________________________________________________________________________
1171Int_t AliTRDCalibraFit::GetNumberOfGroupsPRF(const char* nametitle)
1172{
1173 //
1174 // Get numberofgroupsprf
1175 //
1176
1177 // Some patterns
1178 const Char_t *pattern0 = "Ngp0";
1179 const Char_t *pattern1 = "Ngp1";
1180 const Char_t *pattern2 = "Ngp2";
1181 const Char_t *pattern3 = "Ngp3";
1182 const Char_t *pattern4 = "Ngp4";
1183 const Char_t *pattern5 = "Ngp5";
1184 const Char_t *pattern6 = "Ngp6";
1185
1186 // Nrphi mode
1187 if (strstr(nametitle,pattern0)) {
1188 return 0;
1189 }
1190 if (strstr(nametitle,pattern1)) {
1191 return 1;
1192 }
1193 if (strstr(nametitle,pattern2)) {
1194 return 2;
1195 }
1196 if (strstr(nametitle,pattern3)) {
1197 return 3;
1198 }
1199 if (strstr(nametitle,pattern4)) {
1200 return 4;
1201 }
1202 if (strstr(nametitle,pattern5)) {
1203 return 5;
1204 }
1205 if (strstr(nametitle,pattern6)){
1206 return 6;
1207 }
1208 else return -1;
1209
1210
1211}
1212//_____________________________________________________________________________
1213Bool_t AliTRDCalibraFit::SetModeCalibration(const char *name, Int_t i)
1214{
1215 //
1216 // Set fNz[i] and fNrphi[i] of the AliTRDCalibraFit::Instance()
1217 // corresponding to the given name
1218 //
1219
1220 if(!SetNzFromTObject(name,i)) return kFALSE;
1221 if(!SetNrphiFromTObject(name,i)) return kFALSE;
1222
1223 return kTRUE;
1224
1225}
1226//_____________________________________________________________________________
1227Bool_t AliTRDCalibraFit::SetNrphiFromTObject(const char *name, Int_t i)
1228{
1229 //
1230 // Set fNrphi[i] of the AliTRDCalibraFit::Instance()
1231 // corresponding to the given TObject
1232 //
1233
1234 // Some patterns
1235 const Char_t *patternrphi0 = "Nrphi0";
1236 const Char_t *patternrphi1 = "Nrphi1";
1237 const Char_t *patternrphi2 = "Nrphi2";
1238 const Char_t *patternrphi3 = "Nrphi3";
1239 const Char_t *patternrphi4 = "Nrphi4";
1240 const Char_t *patternrphi5 = "Nrphi5";
1241 const Char_t *patternrphi6 = "Nrphi6";
1242
1243 // Nrphi mode
1244 if (strstr(name,patternrphi0)) {
1245 fCalibraMode->SetNrphi(i ,0);
1246 return kTRUE;
1247 }
1248 if (strstr(name,patternrphi1)) {
1249 fCalibraMode->SetNrphi(i, 1);
1250 return kTRUE;
1251 }
1252 if (strstr(name,patternrphi2)) {
1253 fCalibraMode->SetNrphi(i, 2);
1254 return kTRUE;
1255 }
1256 if (strstr(name,patternrphi3)) {
1257 fCalibraMode->SetNrphi(i, 3);
1258 return kTRUE;
1259 }
1260 if (strstr(name,patternrphi4)) {
1261 fCalibraMode->SetNrphi(i, 4);
1262 return kTRUE;
1263 }
1264 if (strstr(name,patternrphi5)) {
1265 fCalibraMode->SetNrphi(i, 5);
1266 return kTRUE;
1267 }
1268 if (strstr(name,patternrphi6)) {
1269 fCalibraMode->SetNrphi(i, 6);
1270 return kTRUE;
1271 }
1272
1273 fCalibraMode->SetNrphi(i ,0);
1274 return kFALSE;
1275
1276}
1277//_____________________________________________________________________________
1278Bool_t AliTRDCalibraFit::SetNzFromTObject(const char *name, Int_t i)
1279{
1280 //
1281 // Set fNz[i] of the AliTRDCalibraFit::Instance()
1282 // corresponding to the given TObject
1283 //
1284
1285 // Some patterns
1286 const Char_t *patternz0 = "Nz0";
1287 const Char_t *patternz1 = "Nz1";
1288 const Char_t *patternz2 = "Nz2";
1289 const Char_t *patternz3 = "Nz3";
1290 const Char_t *patternz4 = "Nz4";
1291
1292 if (strstr(name,patternz0)) {
1293 fCalibraMode->SetNz(i, 0);
1294 return kTRUE;
1295 }
1296 if (strstr(name,patternz1)) {
1297 fCalibraMode->SetNz(i ,1);
1298 return kTRUE;
1299 }
1300 if (strstr(name,patternz2)) {
1301 fCalibraMode->SetNz(i ,2);
1302 return kTRUE;
1303 }
1304 if (strstr(name,patternz3)) {
1305 fCalibraMode->SetNz(i ,3);
1306 return kTRUE;
1307 }
1308 if (strstr(name,patternz4)) {
1309 fCalibraMode->SetNz(i ,4);
1310 return kTRUE;
1311 }
1312
1313 fCalibraMode->SetNz(i ,0);
1314 return kFALSE;
1315}
1316//_____________________________________________________________________________
1317AliTRDCalDet *AliTRDCalibraFit::CreateDetObjectVdrift(TObjArray *vectorFit, Bool_t perdetector)
1318{
1319 //
1320 // It creates the AliTRDCalDet object from the AliTRDFitInfo
1321 // It takes the mean value of the coefficients per detector
1322 // This object has to be written in the database
1323 //
1324
1325 // Create the DetObject
1326 AliTRDCalDet *object = new AliTRDCalDet("ChamberVdrift","TRD drift velocities (detector value)");
1327
1328 Int_t loop = (Int_t) vectorFit->GetEntriesFast();
1329 if(loop != 540) AliInfo("The Vector Fit is not complete!");
1330 Int_t detector = -1;
1331 Float_t value = 0.0;
1332
1333 for (Int_t k = 0; k < loop; k++) {
1334 detector = ((AliTRDFitInfo *) vectorFit->At(k))->GetDetector();
1335 Float_t mean = 0.0;
1336 if(perdetector){
1337 mean = TMath::Abs(((AliTRDFitInfo *) vectorFit->At(k))->GetCoef()[0]);
1338 }
1339 else {
1340 Int_t count = 0;
1341 Int_t rowMax = fGeo->GetRowMax(GetPlane(detector),GetChamber(detector),GetSector(detector));
1342 Int_t colMax = fGeo->GetColMax(GetPlane(detector));
1343 for (Int_t row = 0; row < rowMax; row++) {
1344 for (Int_t col = 0; col < colMax; col++) {
1345 value = ((AliTRDFitInfo *) vectorFit->At(k))->GetCoef()[(Int_t)(col*rowMax+row)];
1346 mean += TMath::Abs(value);
1347 count++;
1348 } // Col
1349 } // Row
1350 if(count > 0) mean = mean/count;
1351 }
1352 object->SetValue(detector,mean);
1353 }
1354
1355 return object;
1356}
1357//_____________________________________________________________________________
1358AliTRDCalDet *AliTRDCalibraFit::CreateDetObjectGain(TObjArray *vectorFit, Double_t scaleFitFactor, Bool_t perdetector)
1359{
1360 //
1361 // It creates the AliTRDCalDet object from the AliTRDFitInfo
1362 // It takes the mean value of the coefficients per detector
1363 // This object has to be written in the database
1364 //
1365
1366 // Create the DetObject
1367 AliTRDCalDet *object = new AliTRDCalDet("ChamberGainFactor","GainFactor (detector value)");
1368
1369
1370 Int_t loop = (Int_t) vectorFit->GetEntriesFast();
1371 if(loop != 540) AliInfo("The Vector Fit is not complete!");
1372 Int_t detector = -1;
1373 Float_t value = 0.0;
1374
1375 for (Int_t k = 0; k < loop; k++) {
1376 detector = ((AliTRDFitInfo *) vectorFit->At(k))->GetDetector();
1377 Float_t mean = 0.0;
1378 if(perdetector){
1379 value = ((AliTRDFitInfo *) vectorFit->At(k))->GetCoef()[0];
1380 if(value > 0) value = value*scaleFitFactor;
1381 mean = TMath::Abs(value);
1382 }
1383 else{
1384 Int_t count = 0;
1385 Int_t rowMax = fGeo->GetRowMax(GetPlane(detector),GetChamber(detector),GetSector(detector));
1386 Int_t colMax = fGeo->GetColMax(GetPlane(detector));
1387 for (Int_t row = 0; row < rowMax; row++) {
1388 for (Int_t col = 0; col < colMax; col++) {
1389 value = ((AliTRDFitInfo *) vectorFit->At(k))->GetCoef()[(Int_t)(col*rowMax+row)];
1390 if(value > 0) value = value*scaleFitFactor;
1391 mean += TMath::Abs(value);
1392 count++;
1393 } // Col
1394 } // Row
1395 if(count > 0) mean = mean/count;
1396 }
1397 object->SetValue(detector,mean);
1398 }
1399
1400 return object;
1401}
1402//_____________________________________________________________________________
1403AliTRDCalDet *AliTRDCalibraFit::CreateDetObjectT0(TObjArray *vectorFit, Bool_t perdetector)
1404{
1405 //
1406 // It creates the AliTRDCalDet object from the AliTRDFitInfo2
1407 // It takes the min value of the coefficients per detector
1408 // This object has to be written in the database
1409 //
1410
1411 // Create the DetObject
1412 AliTRDCalDet *object = new AliTRDCalDet("ChamberT0","T0 (detector value)");
1413
1414 Int_t loop = (Int_t) vectorFit->GetEntriesFast();
1415 if(loop != 540) AliInfo("The Vector Fit is not complete!");
1416 Int_t detector = -1;
1417 Float_t value = 0.0;
1418
1419 for (Int_t k = 0; k < loop; k++) {
1420 detector = ((AliTRDFitInfo *) vectorFit->At(k))->GetDetector();
1421 Float_t min = 100.0;
1422 if(perdetector){
1423 min = ((AliTRDFitInfo *) vectorFit->At(k))->GetCoef()[0];
1424 }
1425 else{
1426 Int_t rowMax = fGeo->GetRowMax(GetPlane(detector),GetChamber(detector),GetSector(detector));
1427 Int_t colMax = fGeo->GetColMax(GetPlane(detector));
1428 for (Int_t row = 0; row < rowMax; row++) {
1429 for (Int_t col = 0; col < colMax; col++) {
1430 value = ((AliTRDFitInfo *) vectorFit->At(k))->GetCoef()[(Int_t)(col*rowMax+row)];
1431 if(min > value) min = value;
1432 } // Col
1433 } // Row
1434 }
1435 object->SetValue(detector,min);
1436 }
1437
1438 return object;
1439
1440}
1441//_____________________________________________________________________________
1442AliTRDCalDet *AliTRDCalibraFit::CreateDetObjectLorentzAngle(TObjArray *vectorFit)
1443{
1444 //
1445 // It creates the AliTRDCalDet object from the AliTRDFitInfo2
1446 // It takes the min value of the coefficients per detector
1447 // This object has to be written in the database
1448 //
1449
1450 // Create the DetObject
1451 AliTRDCalDet *object = new AliTRDCalDet("tan(lorentzangle)","tan(lorentzangle) (detector value)");
1452
1453
1454 Int_t loop = (Int_t) vectorFit->GetEntriesFast();
1455 if(loop != 540) AliInfo("The Vector Fit is not complete!");
1456 Int_t detector = -1;
1457 Float_t value = 0.0;
1458
1459 for (Int_t k = 0; k < loop; k++) {
1460 detector = ((AliTRDFitInfo *) vectorFit->At(k))->GetDetector();
1461 /*
1462 Int_t rowMax = fGeo->GetRowMax(GetPlane(detector),GetChamber(detector),GetSector(detector));
1463 Int_t colMax = fGeo->GetColMax(GetPlane(detector));
1464 Float_t min = 100.0;
1465 for (Int_t row = 0; row < rowMax; row++) {
1466 for (Int_t col = 0; col < colMax; col++) {
1467 value = ((AliTRDFitInfo *) fVectorFit2.At(k))->GetCoef()[(Int_t)(col*rowMax+row)];
1468 mean += -TMath::Abs(value);
1469 count++;
1470 } // Col
1471 } // Row
1472 if(count > 0) mean = mean/count;
1473 */
1474 value = ((AliTRDFitInfo *) vectorFit->At(k))->GetCoef()[0];
1475 object->SetValue(detector,-TMath::Abs(value));
1476 }
1477
1478 return object;
1479
1480}
1481//_____________________________________________________________________________
1482TObject *AliTRDCalibraFit::CreatePadObjectGain(TObjArray *vectorFit, Double_t scaleFitFactor, AliTRDCalDet *detobject)
1483{
1484 //
1485 // It Creates the AliTRDCalPad object from AliTRDFitInfo
1486 // You need first to create the object for the detectors,
1487 // where the mean value is put.
1488 // This object has to be written in the database
1489 //
1490
1491 // Create the DetObject
1492 AliTRDCalPad *object = new AliTRDCalPad("GainFactor","GainFactor (local variations)");
1493
1494 if(!vectorFit){
1495 for(Int_t k = 0; k < 540; k++){
1496 AliTRDCalROC *calROC = object->GetCalROC(k);
1497 Int_t nchannels = calROC->GetNchannels();
1498 for(Int_t ch = 0; ch < nchannels; ch++){
1499 calROC->SetValue(ch,1.0);
1500 }
1501 }
1502 }
1503 else{
1504
1505 Int_t loop = (Int_t) vectorFit->GetEntriesFast();
1506 if(loop != 540) AliInfo("The Vector Fit is not complete!");
1507 Int_t detector = -1;
1508 Float_t value = 0.0;
1509
1510 for (Int_t k = 0; k < loop; k++) {
1511 detector = ((AliTRDFitInfo *) vectorFit->At(k))->GetDetector();
1512 AliTRDCalROC *calROC = object->GetCalROC(detector);
1513 Float_t mean = detobject->GetValue(detector);
1514 if(mean == 0) continue;
1515 Int_t rowMax = calROC->GetNrows();
1516 Int_t colMax = calROC->GetNcols();
1517 for (Int_t row = 0; row < rowMax; row++) {
1518 for (Int_t col = 0; col < colMax; col++) {
1519 value = ((AliTRDFitInfo *) vectorFit->At(k))->GetCoef()[(Int_t)(col*rowMax+row)];
1520 if(value > 0) value = value*scaleFitFactor;
1521 calROC->SetValue(col,row,TMath::Abs(value)/mean);
1522 } // Col
1523 } // Row
1524 }
1525 }
1526
1527 return object;
1528}
1529//_____________________________________________________________________________
1530TObject *AliTRDCalibraFit::CreatePadObjectVdrift(TObjArray *vectorFit, AliTRDCalDet *detobject)
1531{
1532 //
1533 // It Creates the AliTRDCalPad object from AliTRDFitInfo
1534 // You need first to create the object for the detectors,
1535 // where the mean value is put.
1536 // This object has to be written in the database
1537 //
1538
1539 // Create the DetObject
1540 AliTRDCalPad *object = new AliTRDCalPad("LocalVdrift","TRD drift velocities (local variations)");
1541
1542 if(!vectorFit){
1543 for(Int_t k = 0; k < 540; k++){
1544 AliTRDCalROC *calROC = object->GetCalROC(k);
1545 Int_t nchannels = calROC->GetNchannels();
1546 for(Int_t ch = 0; ch < nchannels; ch++){
1547 calROC->SetValue(ch,1.0);
1548 }
1549 }
1550 }
1551 else {
1552
1553 Int_t loop = (Int_t) vectorFit->GetEntriesFast();
1554 if(loop != 540) AliInfo("The Vector Fit is not complete!");
1555 Int_t detector = -1;
1556 Float_t value = 0.0;
1557
1558 for (Int_t k = 0; k < loop; k++) {
1559 detector = ((AliTRDFitInfo *) vectorFit->At(k))->GetDetector();
1560 AliTRDCalROC *calROC = object->GetCalROC(detector);
1561 Float_t mean = detobject->GetValue(detector);
1562 if(mean == 0) continue;
1563 Int_t rowMax = calROC->GetNrows();
1564 Int_t colMax = calROC->GetNcols();
1565 for (Int_t row = 0; row < rowMax; row++) {
1566 for (Int_t col = 0; col < colMax; col++) {
1567 value = ((AliTRDFitInfo *) vectorFit->At(k))->GetCoef()[(Int_t)(col*rowMax+row)];
1568 calROC->SetValue(col,row,TMath::Abs(value)/mean);
1569 } // Col
1570 } // Row
1571 }
1572 }
1573 return object;
1574
1575}
1576//_____________________________________________________________________________
1577TObject *AliTRDCalibraFit::CreatePadObjectT0(TObjArray *vectorFit, AliTRDCalDet *detobject)
1578{
1579 //
1580 // It Creates the AliTRDCalPad object from AliTRDFitInfo2
1581 // You need first to create the object for the detectors,
1582 // where the mean value is put.
1583 // This object has to be written in the database
1584 //
1585
1586 // Create the DetObject
1587 AliTRDCalPad *object = new AliTRDCalPad("LocalT0","T0 (local variations)");
1588
1589 if(!vectorFit){
1590 for(Int_t k = 0; k < 540; k++){
1591 AliTRDCalROC *calROC = object->GetCalROC(k);
1592 Int_t nchannels = calROC->GetNchannels();
1593 for(Int_t ch = 0; ch < nchannels; ch++){
1594 calROC->SetValue(ch,0.0);
1595 }
1596 }
1597 }
1598 else {
1599
1600 Int_t loop = (Int_t) vectorFit->GetEntriesFast();
1601 if(loop != 540) AliInfo("The Vector Fit is not complete!");
1602 Int_t detector = -1;
1603 Float_t value = 0.0;
1604
1605 for (Int_t k = 0; k < loop; k++) {
1606 detector = ((AliTRDFitInfo *) vectorFit->At(k))->GetDetector();
1607 AliTRDCalROC *calROC = object->GetCalROC(detector);
1608 Float_t min = detobject->GetValue(detector);
1609 Int_t rowMax = calROC->GetNrows();
1610 Int_t colMax = calROC->GetNcols();
1611 for (Int_t row = 0; row < rowMax; row++) {
1612 for (Int_t col = 0; col < colMax; col++) {
1613 value = ((AliTRDFitInfo *) vectorFit->At(k))->GetCoef()[(Int_t)(col*rowMax+row)];
1614 calROC->SetValue(col,row,value-min);
1615 } // Col
1616 } // Row
1617 }
1618 }
1619 return object;
1620
1621}
1622//_____________________________________________________________________________
1623TObject *AliTRDCalibraFit::CreatePadObjectPRF(TObjArray *vectorFit)
1624{
1625 //
1626 // It Creates the AliTRDCalPad object from AliTRDFitInfo
1627 // This object has to be written in the database
1628 //
1629
1630 // Create the DetObject
1631 AliTRDCalPad *object = new AliTRDCalPad("PRFWidth","PRFWidth");
1632
1633 Int_t loop = (Int_t) vectorFit->GetEntriesFast();
1634 if(loop != 540) AliInfo("The Vector Fit is not complete!");
1635 Int_t detector = -1;
1636 Float_t value = 0.0;
1637
1638 for (Int_t k = 0; k < loop; k++) {
1639 detector = ((AliTRDFitInfo *) vectorFit->At(k))->GetDetector();
1640 AliTRDCalROC *calROC = object->GetCalROC(detector);
1641 Int_t rowMax = calROC->GetNrows();
1642 Int_t colMax = calROC->GetNcols();
1643 for (Int_t row = 0; row < rowMax; row++) {
1644 for (Int_t col = 0; col < colMax; col++) {
1645 value = ((AliTRDFitInfo *) vectorFit->At(k))->GetCoef()[(Int_t)(col*rowMax+row)];
1646 calROC->SetValue(col,row,TMath::Abs(value));
1647 } // Col
1648 } // Row
1649 }
1650
1651 return object;
1652
1653}
1654//_____________________________________________________________________________
1655AliTRDCalDet *AliTRDCalibraFit::MakeOutliersStatDet(TObjArray *vectorFit, const char *name, Double_t &mean)
1656{
1657 //
1658 // It Creates the AliTRDCalDet object from AliTRDFitInfo
1659 // 0 successful fit 1 not successful fit
1660 // mean is the mean value over the successful fit
1661 // do not use it for t0: no meaning
1662 //
1663
1664 // Create the CalObject
1665 AliTRDCalDet *object = new AliTRDCalDet(name,name);
1666 mean = 0.0;
1667 Int_t count = 0;
1668
1669 Int_t loop = (Int_t) vectorFit->GetEntriesFast();
1670 if(loop != 540) {
1671 AliInfo("The Vector Fit is not complete! We initialise all outliers");
1672 for(Int_t k = 0; k < 540; k++){
1673 object->SetValue(k,1.0);
1674 }
1675 }
1676 Int_t detector = -1;
1677 Float_t value = 0.0;
1678
1679 for (Int_t k = 0; k < loop; k++) {
1680 detector = ((AliTRDFitInfo *) vectorFit->At(k))->GetDetector();
1681 value = ((AliTRDFitInfo *) vectorFit->At(k))->GetCoef()[0];
1682 if(value <= 0) object->SetValue(detector,1.0);
1683 else {
1684 object->SetValue(detector,0.0);
1685 mean += value;
1686 count++;
1687 }
1688 }
1689 if(count > 0) mean /= count;
1690 return object;
1691}
1692//_____________________________________________________________________________
1693TObject *AliTRDCalibraFit::MakeOutliersStatPad(TObjArray *vectorFit, const char *name, Double_t &mean)
1694{
1695 //
1696 // It Creates the AliTRDCalPad object from AliTRDFitInfo
1697 // 0 not successful fit 1 successful fit
1698 // mean mean value over the successful fit
1699 //
1700
1701 // Create the CalObject
1702 AliTRDCalPad *object = new AliTRDCalPad(name,name);
1703 mean = 0.0;
1704 Int_t count = 0;
1705
1706 Int_t loop = (Int_t) vectorFit->GetEntriesFast();
1707 if(loop != 540) {
1708 AliInfo("The Vector Fit is not complete! We initialise all outliers");
1709 for(Int_t k = 0; k < 540; k++){
1710 AliTRDCalROC *calROC = object->GetCalROC(k);
1711 Int_t nchannels = calROC->GetNchannels();
1712 for(Int_t ch = 0; ch < nchannels; ch++){
1713 calROC->SetValue(ch,1.0);
1714 }
1715 }
1716 }
1717 Int_t detector = -1;
1718 Float_t value = 0.0;
1719
1720 for (Int_t k = 0; k < loop; k++) {
1721 detector = ((AliTRDFitInfo *) vectorFit->At(k))->GetDetector();
1722 AliTRDCalROC *calROC = object->GetCalROC(detector);
1723 Int_t nchannels = calROC->GetNchannels();
1724 for (Int_t ch = 0; ch < nchannels; ch++) {
1725 value = ((AliTRDFitInfo *) vectorFit->At(k))->GetCoef()[ch];
1726 if(value <= 0) calROC->SetValue(ch,1.0);
1727 else {
1728 calROC->SetValue(ch,0.0);
1729 mean += value;
1730 count++;
1731 }
1732 } // channels
1733 }
1734 if(count > 0) mean /= count;
1735 return object;
1736}
1737//_____________________________________________________________________________
1738void AliTRDCalibraFit::SetPeriodeFitPH(Int_t periodeFitPH)
1739{
1740 //
1741 // Set FitPH if 1 then each detector will be fitted
1742 //
1743
1744 if (periodeFitPH > 0) {
1745 fFitPHPeriode = periodeFitPH;
1746 }
1747 else {
1748 AliInfo("periodeFitPH must be higher than 0!");
1749 }
1750
1751}
1752//_____________________________________________________________________________
1753void AliTRDCalibraFit::SetBeginFitCharge(Float_t beginFitCharge)
1754{
1755 //
1756 // The fit of the deposited charge distribution begins at
1757 // histo->Mean()/beginFitCharge
1758 // You can here set beginFitCharge
1759 //
1760
1761 if (beginFitCharge > 0) {
1762 fBeginFitCharge = beginFitCharge;
1763 }
1764 else {
1765 AliInfo("beginFitCharge must be strict positif!");
1766 }
1767
1768}
1769
1770//_____________________________________________________________________________
1771void AliTRDCalibraFit::SetT0Shift0(Float_t t0Shift)
1772{
1773 //
1774 // The t0 calculated with the maximum positif slope is shift from t0Shift0
1775 // You can here set t0Shift0
1776 //
1777
1778 if (t0Shift > 0) {
1779 fT0Shift0 = t0Shift;
1780 }
1781 else {
1782 AliInfo("t0Shift0 must be strict positif!");
1783 }
1784
1785}
1786
1787//_____________________________________________________________________________
1788void AliTRDCalibraFit::SetT0Shift1(Float_t t0Shift)
1789{
1790 //
1791 // The t0 calculated with the maximum of the amplification region is shift from t0Shift1
1792 // You can here set t0Shift1
1793 //
1794
1795 if (t0Shift > 0) {
1796 fT0Shift1 = t0Shift;
1797 }
1798 else {
1799 AliInfo("t0Shift must be strict positif!");
1800 }
1801
1802}
1803
1804//_____________________________________________________________________________
1805void AliTRDCalibraFit::SetRangeFitPRF(Float_t rangeFitPRF)
1806{
1807 //
1808 // The fit of the PRF is from -rangeFitPRF to rangeFitPRF
1809 // You can here set rangeFitPRF
1810 //
1811
1812 if ((rangeFitPRF > 0) &&
1813 (rangeFitPRF <= 1.5)) {
1814 fRangeFitPRF = rangeFitPRF;
1815 }
1816 else {
1817 AliInfo("rangeFitPRF must be between 0 and 1.0");
1818 }
1819
1820}
1821
1822//_____________________________________________________________________________
1823void AliTRDCalibraFit::SetMinEntries(Int_t minEntries)
1824{
1825 //
1826 // Minimum entries for fitting
1827 //
1828
1829 if (minEntries > 0) {
1830 fMinEntries = minEntries;
1831 }
1832 else {
1833 AliInfo("fMinEntries must be >= 0.");
1834 }
1835
1836}
1837
1838//_____________________________________________________________________________
1839void AliTRDCalibraFit::SetRebin(Short_t rebin)
1840{
1841 //
1842 // Rebin with rebin time less bins the Ch histo
1843 // You can set here rebin that should divide the number of bins of CH histo
1844 //
1845
1846 if (rebin > 0) {
1847 fRebin = rebin;
1848 AliInfo("You have to be sure that fRebin divides fNumberBinCharge used!");
1849 }
1850 else {
1851 AliInfo("You have to choose a positiv value!");
1852 }
1853
1854}
1855//_____________________________________________________________________________
1856Bool_t AliTRDCalibraFit::FillVectorFit()
1857{
1858 //
1859 // For the Fit functions fill the vector Fit
1860 //
1861
1862 AliTRDFitInfo *fitInfo = new AliTRDFitInfo();
1863
1864 Int_t ntotal = 1;
1865 if (GetChamber(fCountDet) == 2) {
1866 ntotal = 1728;
1867 }
1868 else {
1869 ntotal = 2304;
1870 }
1871
1872 //printf("For the detector %d , ntotal %d and fCoefCH[0] %f\n",countdet,ntotal,fCoefCH[0]);
1873 Float_t *coef = new Float_t[ntotal];
1874 for (Int_t i = 0; i < ntotal; i++) {
1875 coef[i] = fCurrentCoefDetector[i];
1876 }
1877
1878 Int_t detector = fCountDet;
1879 // Set
1880 fitInfo->SetCoef(coef);
1881 fitInfo->SetDetector(detector);
1882 fVectorFit.Add((TObject *) fitInfo);
1883
1884 return kTRUE;
1885
1886}
1887//_____________________________________________________________________________
1888Bool_t AliTRDCalibraFit::FillVectorFit2()
1889{
1890 //
1891 // For the Fit functions fill the vector Fit
1892 //
1893
1894 AliTRDFitInfo *fitInfo = new AliTRDFitInfo();
1895
1896 Int_t ntotal = 1;
1897 if (GetChamber(fCountDet) == 2) {
1898 ntotal = 1728;
1899 }
1900 else {
1901 ntotal = 2304;
1902 }
1903
1904 //printf("For the detector %d , ntotal %d and fCoefCH[0] %f\n",countdet,ntotal,fCoefCH[0]);
1905 Float_t *coef = new Float_t[ntotal];
1906 for (Int_t i = 0; i < ntotal; i++) {
1907 coef[i] = fCurrentCoefDetector2[i];
1908 }
1909
1910 Int_t detector = fCountDet;
1911 // Set
1912 fitInfo->SetCoef(coef);
1913 fitInfo->SetDetector(detector);
1914 fVectorFit2.Add((TObject *) fitInfo);
1915
1916 return kTRUE;
1917
1918}
1919//____________Functions for initialising the AliTRDCalibraFit in the code_________
1920Bool_t AliTRDCalibraFit::InitFit(Int_t nbins, Int_t i)
1921{
1922 //
1923 // Init the number of expected bins and fDect1[i] fDect2[i]
1924 //
1925
1926 gStyle->SetPalette(1);
1927 gStyle->SetOptStat(1111);
1928 gStyle->SetPadBorderMode(0);
1929 gStyle->SetCanvasColor(10);
1930 gStyle->SetPadLeftMargin(0.13);
1931 gStyle->SetPadRightMargin(0.01);
1932
1933 // Mode groups of pads: the total number of bins!
1934 CalculNumberOfBinsExpected(i);
1935
1936 // Quick verification that we have the good pad calibration mode!
1937 if (fNumberOfBinsExpected != nbins) {
1938 AliInfo("It doesn't correspond to the mode of pad group calibration!");
1939 return kFALSE;
1940 }
1941
1942 // Security for fDebug 3 and 4
1943 if ((fDebugLevel >= 3) &&
1944 ((fDet[0] > 5) ||
1945 (fDet[1] > 4) ||
1946 (fDet[2] > 17))) {
1947 AliInfo("This detector doesn't exit!");
1948 return kFALSE;
1949 }
1950
1951 // Determine fDet1 and fDet2 and set the fNfragZ and fNfragRphi for debug 3 and 4
1952 CalculDect1Dect2(i);
1953
1954
1955 return kTRUE;
1956}
1957//____________Functions for initialising the AliTRDCalibraFit in the code_________
1958Bool_t AliTRDCalibraFit::InitFitCH()
1959{
1960 //
1961 // Init the fVectorFitCH for normalisation
1962 // Init the histo for debugging
1963 //
1964
1965 gDirectory = gROOT;
1966
1967 fScaleFitFactor = 0.0;
1968 fCurrentCoefDetector = new Float_t[2304];
1969 for (Int_t k = 0; k < 2304; k++) {
1970 fCurrentCoefDetector[k] = 0.0;
1971 }
1972 fVectorFit.SetName("gainfactorscoefficients");
1973
1974 // fDebug == 0 nothing
1975 // fDebug == 1 and fFitVoir no histo
1976 if (fDebugLevel == 1) {
1977 if(!CheckFitVoir()) return kFALSE;
1978 }
1979 //Get the CalDet object
1980 if(fAccCDB){
1981 AliTRDcalibDB *cal = AliTRDcalibDB::Instance();
1982 if (!cal) {
1983 AliInfo("Could not get calibDB");
1984 return kFALSE;
1985 }
1986 if(fCalDet) delete fCalDet;
1987 fCalDet = new AliTRDCalDet(*(cal->GetGainFactorDet()));
1988 }
1989 else{
1990 Float_t devalue = 1.0;
1991 if(fCalDet) delete fCalDet;
1992 fCalDet = new AliTRDCalDet("ChamberGainFactor","GainFactor (detector value)");
1993 for(Int_t k = 0; k < 540; k++){
1994 fCalDet->SetValue(k,devalue);
1995 }
1996 }
1997 return kTRUE;
1998
1999}
2000//____________Functions for initialising the AliTRDCalibraFit in the code_________
2001Bool_t AliTRDCalibraFit::InitFitPH()
2002{
2003 //
2004 // Init the arrays of results
2005 // Init the histos for debugging
2006 //
2007
2008 gDirectory = gROOT;
2009 fVectorFit.SetName("driftvelocitycoefficients");
2010 fVectorFit2.SetName("t0coefficients");
2011
2012 fCurrentCoefDetector = new Float_t[2304];
2013 for (Int_t k = 0; k < 2304; k++) {
2014 fCurrentCoefDetector[k] = 0.0;
2015 }
2016
2017 fCurrentCoefDetector2 = new Float_t[2304];
2018 for (Int_t k = 0; k < 2304; k++) {
2019 fCurrentCoefDetector2[k] = 0.0;
2020 }
2021
2022 //fDebug == 0 nothing
2023 // fDebug == 1 and fFitVoir no histo
2024 if (fDebugLevel == 1) {
2025 if(!CheckFitVoir()) return kFALSE;
2026 }
2027 //Get the CalDet object
2028 if(fAccCDB){
2029 AliTRDcalibDB *cal = AliTRDcalibDB::Instance();
2030 if (!cal) {
2031 AliInfo("Could not get calibDB");
2032 return kFALSE;
2033 }
2034 if(fCalDet) delete fCalDet;
2035 if(fCalDet2) delete fCalDet2;
2036 fCalDet = new AliTRDCalDet(*(cal->GetVdriftDet()));
2037 fCalDet2 = new AliTRDCalDet(*(cal->GetT0Det()));
2038 }
2039 else{
2040 Float_t devalue = 1.5;
2041 Float_t devalue2 = 0.0;
2042 if(fCalDet) delete fCalDet;
2043 if(fCalDet2) delete fCalDet2;
2044 fCalDet = new AliTRDCalDet("ChamberVdrift","TRD drift velocities (detector value)");
2045 fCalDet2 = new AliTRDCalDet("ChamberT0","T0 (detector value)");
2046 for(Int_t k = 0; k < 540; k++){
2047 fCalDet->SetValue(k,devalue);
2048 fCalDet2->SetValue(k,devalue2);
2049 }
2050 }
2051 return kTRUE;
2052}
2053//____________Functions for initialising the AliTRDCalibraFit in the code_________
2054Bool_t AliTRDCalibraFit::InitFitPRF()
2055{
2056 //
2057 // Init the calibration mode (Nz, Nrphi), the histograms for
2058 // debugging the fit methods if fDebug > 0,
2059 //
2060
2061 gDirectory = gROOT;
2062 fVectorFit.SetName("prfwidthcoefficients");
2063
2064 fCurrentCoefDetector = new Float_t[2304];
2065 for (Int_t k = 0; k < 2304; k++) {
2066 fCurrentCoefDetector[k] = 0.0;
2067 }
2068
2069 // fDebug == 0 nothing
2070 // fDebug == 1 and fFitVoir no histo
2071 if (fDebugLevel == 1) {
2072 if(!CheckFitVoir()) return kFALSE;
2073 }
2074 return kTRUE;
2075}
2076//____________Functions for initialising the AliTRDCalibraFit in the code_________
2077Bool_t AliTRDCalibraFit::InitFitLinearFitter()
2078{
2079 //
2080 // Init the fCalDet, fVectorFit fCurrentCoefDetector
2081 //
2082
2083 gDirectory = gROOT;
2084
2085 fCurrentCoefDetector = new Float_t[2304];
2086 fCurrentCoefDetector2 = new Float_t[2304];
2087 for (Int_t k = 0; k < 2304; k++) {
2088 fCurrentCoefDetector[k] = 0.0;
2089 fCurrentCoefDetector2[k] = 0.0;
2090 }
2091
2092 //printf("test0\n");
2093
2094 AliTRDcalibDB *cal = AliTRDcalibDB::Instance();
2095 if (!cal) {
2096 AliInfo("Could not get calibDB");
2097 return kFALSE;
2098 }
2099
2100 //Get the CalDet object
2101 if(fAccCDB){
2102 if(fCalDet) delete fCalDet;
2103 if(fCalDet2) delete fCalDet2;
2104 fCalDet = new AliTRDCalDet(*(cal->GetVdriftDet()));
2105 //printf("test1\n");
2106 fCalDet2 = new AliTRDCalDet("lorentz angle tan","lorentz angle tan (detector value)");
2107 //printf("test2\n");
2108 for(Int_t k = 0; k < 540; k++){
2109 fCalDet2->SetValue(k,cal->GetOmegaTau(fCalDet->GetValue(k),-fMagneticField));
2110 }
2111 //printf("test3\n");
2112 }
2113 else{
2114 Float_t devalue = 1.5;
2115 Float_t devalue2 = cal->GetOmegaTau(1.5,-fMagneticField);
2116 if(fCalDet) delete fCalDet;
2117 if(fCalDet2) delete fCalDet2;
2118 //printf("test1\n");
2119 fCalDet = new AliTRDCalDet("ChamberVdrift","TRD drift velocities (detector value)");
2120 fCalDet2 = new AliTRDCalDet("lorentz angle tan","lorentz angle tan (detector value)");
2121 //printf("test2\n");
2122 for(Int_t k = 0; k < 540; k++){
2123 fCalDet->SetValue(k,devalue);
2124 fCalDet2->SetValue(k,devalue2);
2125 }
2126 //printf("test3\n");
2127 }
2128 return kTRUE;
2129}
2130
2131//____________Functions for initialising the AliTRDCalibraFit in the code_________
2132void AliTRDCalibraFit::InitfCountDetAndfCount(Int_t i)
2133{
2134 //
2135 // Init the current detector where we are fCountDet and the
2136 // next fCount for the functions Fit...
2137 //
2138
2139 // Loop on the Xbins of ch!!
2140 fCountDet = -1; // Current detector
2141 fCount = 0; // To find the next detector
2142
2143 // If fDebug >= 3
2144 if (fDebugLevel >= 3) {
2145 // Set countdet to the detector
2146 fCountDet = AliTRDgeometry::GetDetector(fDet[0],fDet[1],fDet[2]);
2147 // Set counter to write at the end of the detector
2148 fCount = fDect2;
2149 // Get the right calib objects
2150 SetCalROC(i);
2151 }
2152 if(fDebugLevel == 1) {
2153 fCountDet = 0;
2154 fCalibraMode->CalculXBins(fCountDet,i);
2155 while(fCalibraMode->GetXbins(i) <=fFitVoir){
2156 fCountDet++;
2157 fCalibraMode->CalculXBins(fCountDet,i);
2158 }
2159 fCount = fCalibraMode->GetXbins(i);
2160 fCountDet--;
2161 // Determination of fNnZ, fNnRphi, fNfragZ and fNfragRphi
2162 fCalibraMode->ModePadCalibration((Int_t) GetChamber(fCountDet),i);
2163 fCalibraMode->ModePadFragmentation((Int_t) GetPlane(fCountDet)
2164 ,(Int_t) GetChamber(fCountDet)
2165 ,(Int_t) GetSector(fCountDet),i);
2166 }
2167}
2168//_______________________________________________________________________________
2169void AliTRDCalibraFit::CalculNumberOfBinsExpected(Int_t i)
2170{
2171 //
2172 // Calculate the number of bins expected (calibration groups)
2173 //
2174
2175 fNumberOfBinsExpected = 0;
2176 fCalibraMode->ModePadCalibration(2,i);
2177 fCalibraMode->ModePadFragmentation(0,2,0,i);
2178 fCalibraMode->SetDetChamb2(i);
2179 if (fDebugLevel > 1) {
2180 AliInfo(Form("For the chamber 2: %d",fCalibraMode->GetDetChamb2(i)));
2181 }
2182 fNumberOfBinsExpected += 6 * 18 * fCalibraMode->GetDetChamb2(i);
2183 fCalibraMode->ModePadCalibration(0,i);
2184 fCalibraMode->ModePadFragmentation(0,0,0,i);
2185 fCalibraMode->SetDetChamb0(i);
2186 if (fDebugLevel > 1) {
2187 AliInfo(Form("For the other chamber 0: %d",fCalibraMode->GetDetChamb0(i)));
2188 }
2189 fNumberOfBinsExpected += 6 * 4 * 18 * fCalibraMode->GetDetChamb0(i);
2190
2191}
2192//_______________________________________________________________________________
2193void AliTRDCalibraFit::CalculDect1Dect2(Int_t i)
2194{
2195 //
2196 // Calculate the range of fits
2197 //
2198
2199 fDect1 = -1;
2200 fDect2 = -1;
2201 if (fDebugLevel == 1) {
2202 fDect1 = fFitVoir;
2203 fDect2 = fDect1 +1;
2204 }
2205 if ((fDebugLevel == 2) || (fDebugLevel == 0)) {
2206 fDect1 = 0;
2207 fDect2 = fNumberOfBinsExpected;
2208 }
2209 if (fDebugLevel >= 3) {
2210 fCountDet = AliTRDgeometry::GetDetector(fDet[0],fDet[1],fDet[2]);
2211 fCalibraMode->CalculXBins(fCountDet,i);
2212 fDect1 = fCalibraMode->GetXbins(i);
2213 // Determination of fNnZ, fNnRphi, fNfragZ and fNfragRphi
2214 fCalibraMode->ModePadCalibration((Int_t) GetChamber(fCountDet),i);
2215 fCalibraMode->ModePadFragmentation((Int_t) GetPlane(fCountDet)
2216 ,(Int_t) GetChamber(fCountDet)
2217 ,(Int_t) GetSector(fCountDet),i);
2218 // Set for the next detector
2219 fDect2 = fDect1 + fCalibraMode->GetNfragZ(i)*fCalibraMode->GetNfragRphi(i);
2220 }
2221}
2222//_______________________________________________________________________________
2223Bool_t AliTRDCalibraFit::CheckFitVoir()
2224{
2225 //
2226 // Check if fFitVoir is in the range
2227 //
2228
2229 if (fFitVoir < fNumberOfBinsExpected) {
2230 AliInfo(Form("We will see the fit of the object %d",fFitVoir));
2231 }
2232 else {
2233 AliInfo("fFitVoir is out of range of the histo!");
2234 return kFALSE;
2235 }
2236 return kTRUE;
2237}
2238//____________Functions for initialising the AliTRDCalibraFit in the code_________
2239void AliTRDCalibraFit::UpdatefCountDetAndfCount(Int_t idect, Int_t i)
2240{
2241 //
2242 // See if we are in a new detector and update the
2243 // variables fNfragZ and fNfragRphi if yes
2244 // Will never happen for only one detector (3 and 4)
2245 // Doesn't matter for 2
2246 //
2247 if (fCount == idect) {
2248 // On en est au detector
2249 fCountDet += 1;
2250 // Determination of fNnZ, fNnRphi, fNfragZ and fNfragRphi
2251 fCalibraMode->ModePadCalibration((Int_t) GetChamber(fCountDet),i);
2252 fCalibraMode->ModePadFragmentation((Int_t) GetPlane(fCountDet)
2253 ,(Int_t) GetChamber(fCountDet)
2254 ,(Int_t) GetSector(fCountDet),i);
2255 // Set for the next detector
2256 fCount += fCalibraMode->GetNfragZ(i)*fCalibraMode->GetNfragRphi(i);
2257 // calib objects
2258 SetCalROC(i);
2259 }
2260}
2261//____________Functions for initialising the AliTRDCalibraFit in the code_________
2262void AliTRDCalibraFit::ReconstructFitRowMinRowMax(Int_t idect, Int_t i)
2263{
2264 //
2265 // Reconstruct the min pad row, max pad row, min pad col and
2266 // max pad col of the calibration group for the Fit functions
2267 //
2268 if (fDebugLevel != 1) {
2269 fCalibraMode->ReconstructionRowPadGroup((Int_t) (idect-(fCount-(fCalibraMode->GetNfragZ(i)*fCalibraMode->GetNfragRphi(i)))),i);
2270 }
2271}
2272//____________Functions for initialising the AliTRDCalibraFit in the code_________
2273Bool_t AliTRDCalibraFit::NotEnoughStatisticCH(Int_t idect)
2274{
2275 //
2276 // For the case where there are not enough entries in the histograms
2277 // of the calibration group, the value present in the choosen database
2278 // will be put. A negativ sign enables to know that a fit was not possible.
2279 //
2280
2281 if (fDebugLevel == 1) {
2282 AliInfo("The element has not enough statistic to be fitted");
2283 }
2284
2285 else {
2286
2287 AliInfo(Form("The element %d in this detector %d has not enough statistic to be fitted"
2288 ,idect-(fCount-(fCalibraMode->GetNfragZ(0)*fCalibraMode->GetNfragRphi(0))),fCountDet));
2289
2290 // Calcul the coef from the database choosen
2291 CalculChargeCoefMean(kFALSE);
2292
2293 //chamber 2, not chamber 2
2294 Int_t factor = 0;
2295 if(GetChamber(fCountDet) == 2) factor = 12;
2296 else factor = 16;
2297
2298 // Fill the fCurrentCoefDetector with negative value to say: not fitted
2299 for (Int_t k = fCalibraMode->GetRowMin(0); k < fCalibraMode->GetRowMax(0); k++) {
2300 for (Int_t j = fCalibraMode->GetColMin(0); j < fCalibraMode->GetColMax(0); j++) {
2301 fCurrentCoefDetector[(Int_t)(j*factor+k)] = -TMath::Abs(fCurrentCoef[1]);
2302 }
2303 }
2304
2305 //Put default value negative
2306 fCurrentCoef[0] = -TMath::Abs(fCurrentCoef[1]);
2307 fCurrentCoefE = 0.0;
2308
2309 FillFillCH(idect);
2310 }
2311
2312 return kTRUE;
2313}
2314
2315
2316//____________Functions for initialising the AliTRDCalibraFit in the code_________
2317Bool_t AliTRDCalibraFit::NotEnoughStatisticPH(Int_t idect)
2318{
2319 //
2320 // For the case where there are not enough entries in the histograms
2321 // of the calibration group, the value present in the choosen database
2322 // will be put. A negativ sign enables to know that a fit was not possible.
2323 //
2324 if (fDebugLevel == 1) {
2325 AliInfo("The element has not enough statistic to be fitted");
2326 }
2327 else {
2328
2329 AliInfo(Form("The element %d in this detector %d has not enough statistic to be fitted"
2330 ,idect-(fCount-(fCalibraMode->GetNfragZ(1)*fCalibraMode->GetNfragRphi(1))),fCountDet));
2331
2332 CalculVdriftCoefMean();
2333 CalculT0CoefMean();
2334
2335 //chamber 2 and not chamber 2
2336 Int_t factor = 0;
2337 if(GetChamber(fCountDet) == 2) factor = 12;
2338 else factor = 16;
2339
2340
2341 // Fill the fCurrentCoefDetector 2
2342 for (Int_t k = fCalibraMode->GetRowMin(1); k < fCalibraMode->GetRowMax(1); k++) {
2343 for (Int_t j = fCalibraMode->GetColMin(1); j < fCalibraMode->GetColMax(1); j++) {
2344 fCurrentCoefDetector[(Int_t)(j*factor+k)] = -TMath::Abs(fCurrentCoef[1]);
2345 fCurrentCoefDetector2[(Int_t)(j*factor+k)] = fCurrentCoef2[1];
2346 }
2347 }
2348
2349 // Put the default value
2350 fCurrentCoef[0] = -TMath::Abs(fCurrentCoef[1]);
2351 fCurrentCoefE = 0.0;
2352 fCurrentCoef2[0] = fCurrentCoef2[1];
2353 fCurrentCoefE2 = 0.0;
2354
2355 FillFillPH(idect);
2356
2357 }
2358
2359 return kTRUE;
2360
2361}
2362
2363
2364//____________Functions for initialising the AliTRDCalibraFit in the code_________
2365Bool_t AliTRDCalibraFit::NotEnoughStatisticPRF(Int_t idect)
2366{
2367 //
2368 // For the case where there are not enough entries in the histograms
2369 // of the calibration group, the value present in the choosen database
2370 // will be put. A negativ sign enables to know that a fit was not possible.
2371 //
2372
2373 if (fDebugLevel == 1) {
2374 AliInfo("The element has not enough statistic to be fitted");
2375 }
2376 else {
2377
2378 AliInfo(Form("The element %d in this detector %d has not enough statistic to be fitted"
2379 ,idect-(fCount-(fCalibraMode->GetNfragZ(2)*fCalibraMode->GetNfragRphi(2))),fCountDet));
2380
2381 CalculPRFCoefMean();
2382
2383 // chamber 2 and not chamber 2
2384 Int_t factor = 0;
2385 if(GetChamber(fCountDet) == 2) factor = 12;
2386 else factor = 16;
2387
2388
2389 // Fill the fCurrentCoefDetector
2390 for (Int_t k = fCalibraMode->GetRowMin(2); k < fCalibraMode->GetRowMax(2); k++) {
2391 for (Int_t j = fCalibraMode->GetColMin(2); j < fCalibraMode->GetColMax(2); j++) {
2392 fCurrentCoefDetector[(Int_t)(j*factor+k)] = -fCurrentCoef[1];
2393 }
2394 }
2395
2396 // Put the default value
2397 fCurrentCoef[0] = -fCurrentCoef[1];
2398 fCurrentCoefE = 0.0;
2399
2400 FillFillPRF(idect);
2401 }
2402
2403 return kTRUE;
2404
2405}
2406//____________Functions for initialising the AliTRDCalibraFit in the code_________
2407Bool_t AliTRDCalibraFit::NotEnoughStatisticLinearFitter()
2408{
2409 //
2410 // For the case where there are not enough entries in the histograms
2411 // of the calibration group, the value present in the choosen database
2412 // will be put. A negativ sign enables to know that a fit was not possible.
2413 //
2414
2415 // Calcul the coef from the database choosen
2416 CalculVdriftLorentzCoef();
2417
2418 Int_t factor = 0;
2419 if(GetChamber(fCountDet) == 2) factor = 1728;
2420 else factor = 2304;
2421
2422
2423 // Fill the fCurrentCoefDetector
2424 for (Int_t k = 0; k < factor; k++) {
2425 fCurrentCoefDetector[k] = -TMath::Abs(fCurrentCoef[1]);
2426 // should be negative
2427 fCurrentCoefDetector2[k] = +TMath::Abs(fCurrentCoef2[1]);
2428 }
2429
2430
2431 //Put default opposite sign
2432 fCurrentCoef[0] = -TMath::Abs(fCurrentCoef[1]);
2433 fCurrentCoefE = 0.0;
2434 fCurrentCoef2[0] = +TMath::Abs(fCurrentCoef2[1]);
2435 fCurrentCoefE2 = 0.0;
2436
2437 FillFillLinearFitter();
2438
2439 return kTRUE;
2440}
2441
2442//____________Functions for initialising the AliTRDCalibraFit in the code_________
2443Bool_t AliTRDCalibraFit::FillInfosFitCH(Int_t idect)
2444{
2445 //
2446 // Fill the coefficients found with the fits or other
2447 // methods from the Fit functions
2448 //
2449
2450 if (fDebugLevel != 1) {
2451
2452 Int_t factor = 0;
2453 if(GetChamber(fCountDet) == 2) factor = 12;
2454 else factor = 16;
2455
2456 for (Int_t k = fCalibraMode->GetRowMin(0); k < fCalibraMode->GetRowMax(0); k++) {
2457 for (Int_t j = fCalibraMode->GetColMin(0); j < fCalibraMode->GetColMax(0); j++) {
2458 fCurrentCoefDetector[(Int_t)(j*factor+k)] = fCurrentCoef[0];
2459 }
2460 }
2461
2462 FillFillCH(idect);
2463
2464 }
2465
2466 return kTRUE;
2467
2468}
2469//____________Functions for initialising the AliTRDCalibraFit in the code_________
2470Bool_t AliTRDCalibraFit::FillInfosFitPH(Int_t idect)
2471{
2472 //
2473 // Fill the coefficients found with the fits or other
2474 // methods from the Fit functions
2475 //
2476
2477 if (fDebugLevel != 1) {
2478
2479 Int_t factor = 0;
2480 if(GetChamber(fCountDet) == 2) factor = 12;
2481 else factor = 16;
2482
2483 for (Int_t k = fCalibraMode->GetRowMin(1); k < fCalibraMode->GetRowMax(1); k++) {
2484 for (Int_t j = fCalibraMode->GetColMin(1); j < fCalibraMode->GetColMax(1); j++) {
2485 fCurrentCoefDetector[(Int_t)(j*factor+k)] = fCurrentCoef[0];
2486 fCurrentCoefDetector2[(Int_t)(j*factor+k)] = fCurrentCoef2[0];
2487 }
2488 }
2489 FillFillPH(idect);
2490 }
2491 return kTRUE;
2492}
2493//____________Functions for initialising the AliTRDCalibraFit in the code_________
2494Bool_t AliTRDCalibraFit::FillInfosFitPRF(Int_t idect)
2495{
2496 //
2497 // Fill the coefficients found with the fits or other
2498 // methods from the Fit functions
2499 //
2500
2501 if (fDebugLevel != 1) {
2502
2503 Int_t factor = 0;
2504 if(GetChamber(fCountDet) == 2) factor = 12;
2505 else factor = 16;
2506
2507 // Pointer to the branch
2508 for (Int_t k = fCalibraMode->GetRowMin(2); k < fCalibraMode->GetRowMax(2); k++) {
2509 for (Int_t j = fCalibraMode->GetColMin(2); j < fCalibraMode->GetColMax(2); j++) {
2510 fCurrentCoefDetector[(Int_t)(j*factor+k)] = fCurrentCoef[0];
2511 }
2512 }
2513 FillFillPRF(idect);
2514 }
2515
2516 return kTRUE;
2517
2518}
2519//____________Functions for initialising the AliTRDCalibraFit in the code_________
2520Bool_t AliTRDCalibraFit::FillInfosFitLinearFitter()
2521{
2522 //
2523 // Fill the coefficients found with the fits or other
2524 // methods from the Fit functions
2525 //
2526
2527 Int_t factor = 0;
2528 if(GetChamber(fCountDet) == 2) factor = 1728;
2529 else factor = 2304;
2530
2531 // Pointer to the branch
2532 for (Int_t k = 0; k < factor; k++) {
2533 fCurrentCoefDetector[k] = fCurrentCoef[0];
2534 fCurrentCoefDetector2[k] = fCurrentCoef2[0];
2535 }
2536
2537 FillFillLinearFitter();
2538
2539 return kTRUE;
2540
2541}
2542//________________________________________________________________________________
2543void AliTRDCalibraFit::FillFillCH(Int_t idect)
2544{
2545 //
2546 // DebugStream and fVectorFit
2547 //
2548
2549 // End of one detector
2550 if ((idect == (fCount-1))) {
2551 FillVectorFit();
2552 // Reset
2553 for (Int_t k = 0; k < 2304; k++) {
2554 fCurrentCoefDetector[k] = 0.0;
2555 }
2556 }
2557
2558 if(fDebugLevel > 1){
2559
2560 if ( !fDebugStreamer ) {
2561 //debug stream
2562 TDirectory *backup = gDirectory;
2563 fDebugStreamer = new TTreeSRedirector("TRDDebugFit.root");
2564 if ( backup ) backup->cd(); //we don't want to be cd'd to the debug streamer
2565 }
2566
2567 Int_t detector = fCountDet;
2568 Int_t caligroup = idect;
2569 Short_t rowmin = fCalibraMode->GetRowMin(0);
2570 Short_t rowmax = fCalibraMode->GetRowMax(0);
2571 Short_t colmin = fCalibraMode->GetColMin(0);
2572 Short_t colmax = fCalibraMode->GetColMax(0);
2573 Float_t gf = fCurrentCoef[0];
2574 Float_t gfs = fCurrentCoef[1];
2575 Float_t gfE = fCurrentCoefE;
2576
2577 (*fDebugStreamer) << "FillFillCH" <<
2578 "detector=" << detector <<
2579 "caligroup=" << caligroup <<
2580 "rowmin=" << rowmin <<
2581 "rowmax=" << rowmax <<
2582 "colmin=" << colmin <<
2583 "colmax=" << colmax <<
2584 "gf=" << gf <<
2585 "gfs=" << gfs <<
2586 "gfE=" << gfE <<
2587 "\n";
2588
2589 }
2590}
2591//________________________________________________________________________________
2592void AliTRDCalibraFit::FillFillPH(Int_t idect)
2593{
2594 //
2595 // DebugStream and fVectorFit and fVectorFit2
2596 //
2597
2598 // End of one detector
2599 if ((idect == (fCount-1))) {
2600 FillVectorFit();
2601 FillVectorFit2();
2602 // Reset
2603 for (Int_t k = 0; k < 2304; k++) {
2604 fCurrentCoefDetector[k] = 0.0;
2605 fCurrentCoefDetector2[k] = 0.0;
2606 }
2607 }
2608
2609 if(fDebugLevel > 1){
2610
2611 if ( !fDebugStreamer ) {
2612 //debug stream
2613 TDirectory *backup = gDirectory;
2614 fDebugStreamer = new TTreeSRedirector("TRDDebugFit.root");
2615 if ( backup ) backup->cd(); //we don't want to be cd'd to the debug streamer
2616 }
2617
2618
2619 Int_t detector = fCountDet;
2620 Int_t caligroup = idect;
2621 Short_t rowmin = fCalibraMode->GetRowMin(1);
2622 Short_t rowmax = fCalibraMode->GetRowMax(1);
2623 Short_t colmin = fCalibraMode->GetColMin(1);
2624 Short_t colmax = fCalibraMode->GetColMax(1);
2625 Float_t vf = fCurrentCoef[0];
2626 Float_t vs = fCurrentCoef[1];
2627 Float_t vfE = fCurrentCoefE;
2628 Float_t t0f = fCurrentCoef2[0];
2629 Float_t t0s = fCurrentCoef2[1];
2630 Float_t t0E = fCurrentCoefE2;
2631
2632
2633
2634 (* fDebugStreamer) << "FillFillPH"<<
2635 "detector="<<detector<<
2636 "caligroup="<<caligroup<<
2637 "rowmin="<<rowmin<<
2638 "rowmax="<<rowmax<<
2639 "colmin="<<colmin<<
2640 "colmax="<<colmax<<
2641 "vf="<<vf<<
2642 "vs="<<vs<<
2643 "vfE="<<vfE<<
2644 "t0f="<<t0f<<
2645 "t0s="<<t0s<<
2646 "t0E="<<t0E<<
2647 "\n";
2648 }
2649
2650}
2651//________________________________________________________________________________
2652void AliTRDCalibraFit::FillFillPRF(Int_t idect)
2653{
2654 //
2655 // DebugStream and fVectorFit
2656 //
2657
2658 // End of one detector
2659 if ((idect == (fCount-1))) {
2660 FillVectorFit();
2661 // Reset
2662 for (Int_t k = 0; k < 2304; k++) {
2663 fCurrentCoefDetector[k] = 0.0;
2664 }
2665 }
2666
2667
2668 if(fDebugLevel > 1){
2669
2670 if ( !fDebugStreamer ) {
2671 //debug stream
2672 TDirectory *backup = gDirectory;
2673 fDebugStreamer = new TTreeSRedirector("TRDDebugFit.root");
2674 if ( backup ) backup->cd(); //we don't want to be cd'd to the debug streamer
2675 }
2676
2677 Int_t detector = fCountDet;
2678 Int_t plane = GetPlane(fCountDet);
2679 Int_t caligroup = idect;
2680 Short_t rowmin = fCalibraMode->GetRowMin(2);
2681 Short_t rowmax = fCalibraMode->GetRowMax(2);
2682 Short_t colmin = fCalibraMode->GetColMin(2);
2683 Short_t colmax = fCalibraMode->GetColMax(2);
2684 Float_t widf = fCurrentCoef[0];
2685 Float_t wids = fCurrentCoef[1];
2686 Float_t widfE = fCurrentCoefE;
2687
2688 (* fDebugStreamer) << "FillFillPRF"<<
2689 "detector="<<detector<<
2690 "plane="<<plane<<
2691 "caligroup="<<caligroup<<
2692 "rowmin="<<rowmin<<
2693 "rowmax="<<rowmax<<
2694 "colmin="<<colmin<<
2695 "colmax="<<colmax<<
2696 "widf="<<widf<<
2697 "wids="<<wids<<
2698 "widfE="<<widfE<<
2699 "\n";
2700 }
2701
2702}
2703//________________________________________________________________________________
2704void AliTRDCalibraFit::FillFillLinearFitter()
2705{
2706 //
2707 // DebugStream and fVectorFit
2708 //
2709
2710 // End of one detector
2711 FillVectorFit();
2712 FillVectorFit2();
2713
2714
2715 // Reset
2716 for (Int_t k = 0; k < 2304; k++) {
2717 fCurrentCoefDetector[k] = 0.0;
2718 fCurrentCoefDetector2[k] = 0.0;
2719 }
2720
2721
2722 if(fDebugLevel > 1){
2723
2724 if ( !fDebugStreamer ) {
2725 //debug stream
2726 TDirectory *backup = gDirectory;
2727 fDebugStreamer = new TTreeSRedirector("TRDDebugFit.root");
2728 if ( backup ) backup->cd(); //we don't want to be cd'd to the debug streamer
2729 }
2730
2731 //Debug: comparaison of the different methods (okey for first time but not for iterative procedure)
2732 AliTRDpadPlane *padplane = fGeo->GetPadPlane(GetPlane(fCountDet),GetChamber(fCountDet));
2733 Float_t rowmd = (padplane->GetRow0()+padplane->GetRowEnd())/2.;
2734 Float_t r = AliTRDgeometry::GetTime0(GetPlane(fCountDet));
2735 Float_t tiltangle = padplane->GetTiltingAngle();
2736 Int_t detector = fCountDet;
2737 Int_t chamber = GetChamber(fCountDet);
2738 Int_t plane = GetChamber(fCountDet);
2739 Float_t vf = fCurrentCoef[0];
2740 Float_t vs = fCurrentCoef[1];
2741 Float_t vfE = fCurrentCoefE;
2742 Float_t lorentzangler = fCurrentCoef2[0];
2743 Float_t elorentzangler = fCurrentCoefE2;
2744 Float_t lorentzangles = fCurrentCoef2[1];
2745
2746 (* fDebugStreamer) << "FillFillLinearFitter"<<
2747 "detector="<<detector<<
2748 "chamber="<<chamber<<
2749 "plane="<<plane<<
2750 "rowmd="<<rowmd<<
2751 "r="<<r<<
2752 "tiltangle="<<tiltangle<<
2753 "vf="<<vf<<
2754 "vs="<<vs<<
2755 "vfE="<<vfE<<
2756 "lorentzangler="<<lorentzangler<<
2757 "Elorentzangler="<<elorentzangler<<
2758 "lorentzangles="<<lorentzangles<<
2759 "\n";
2760 }
2761
2762}
2763//
2764//____________Calcul Coef Mean_________________________________________________
2765//
2766//_____________________________________________________________________________
2767Bool_t AliTRDCalibraFit::CalculT0CoefMean()
2768{
2769 //
2770 // For the detector Dect calcul the mean time 0
2771 // for the calibration group idect from the choosen database
2772 //
2773
2774 fCurrentCoef2[1] = 0.0;
2775 if(fDebugLevel != 1){
2776 if((fCalibraMode->GetNz(1) > 0) ||
2777 (fCalibraMode->GetNrphi(1) > 0)) {
2778 for (Int_t row = fCalibraMode->GetRowMin(1); row < fCalibraMode->GetRowMax(1); row++) {
2779 for (Int_t col = fCalibraMode->GetColMin(1); col < fCalibraMode->GetColMax(1); col++) {
2780 fCurrentCoef2[1] += (Float_t) (fCalROC2->GetValue(col,row)+fCalDet2->GetValue(fCountDet));
2781 }
2782 }
2783 fCurrentCoef2[1] = fCurrentCoef2[1] / ((fCalibraMode->GetColMax(1)-fCalibraMode->GetColMin(1))*(fCalibraMode->GetRowMax(1)-fCalibraMode->GetRowMin(1)));
2784 }
2785 else {
2786 if(!fAccCDB){
2787 fCurrentCoef2[1] = fCalDet2->GetValue(fCountDet);
2788 }
2789 else{
2790 for(Int_t row = 0; row < fGeo->GetRowMax(GetPlane(fCountDet),GetChamber(fCountDet),GetSector(fCountDet)); row++){
2791 for(Int_t col = 0; col < fGeo->GetColMax(GetPlane(fCountDet)); col++){
2792 fCurrentCoef2[1] += (Float_t) (fCalROC2->GetValue(col,row)+fCalDet2->GetValue(fCountDet));
2793 }
2794 }
2795 fCurrentCoef2[1] = fCurrentCoef2[1] / ((fGeo->GetRowMax(GetPlane(fCountDet),GetChamber(fCountDet),GetSector(fCountDet)))*(fGeo->GetColMax(GetPlane(fCountDet))));
2796 }
2797 }
2798 }
2799 return kTRUE;
2800}
2801
2802//_____________________________________________________________________________
2803Bool_t AliTRDCalibraFit::CalculChargeCoefMean(Bool_t vrai)
2804{
2805 //
2806 // For the detector Dect calcul the mean gain factor
2807 // for the calibration group idect from the choosen database
2808 //
2809
2810 fCurrentCoef[1] = 0.0;
2811 if(fDebugLevel != 1){
2812 if ((fCalibraMode->GetNz(0) > 0) ||
2813 (fCalibraMode->GetNrphi(0) > 0)) {
2814 for (Int_t row = fCalibraMode->GetRowMin(0); row < fCalibraMode->GetRowMax(0); row++) {
2815 for (Int_t col = fCalibraMode->GetColMin(0); col < fCalibraMode->GetColMax(0); col++) {
2816 fCurrentCoef[1] += (Float_t) (fCalROC->GetValue(col,row)*fCalDet->GetValue(fCountDet));
2817 if (vrai) fScaleFitFactor += (Float_t) (fCalROC->GetValue(col,row)*fCalDet->GetValue(fCountDet));
2818 }
2819 }
2820 fCurrentCoef[1] = fCurrentCoef[1] / ((fCalibraMode->GetColMax(0)-fCalibraMode->GetColMin(0))*(fCalibraMode->GetRowMax(0)-fCalibraMode->GetRowMin(0)));
2821 }
2822 else {
2823 //Per detectors
2824 fCurrentCoef[1] = (Float_t) fCalDet->GetValue(fCountDet);
2825 if (vrai) fScaleFitFactor += ((Float_t) fCalDet->GetValue(fCountDet))*(fCalibraMode->GetColMax(0)-fCalibraMode->GetColMin(0))*(fCalibraMode->GetRowMax(0)-fCalibraMode->GetRowMin(0));
2826 }
2827 }
2828 return kTRUE;
2829}
2830//_____________________________________________________________________________
2831Bool_t AliTRDCalibraFit::CalculPRFCoefMean()
2832{
2833 //
2834 // For the detector Dect calcul the mean sigma of pad response
2835 // function for the calibration group idect from the choosen database
2836 //
2837
2838 fCurrentCoef[1] = 0.0;
2839 if(fDebugLevel != 1){
2840 for (Int_t row = fCalibraMode->GetRowMin(2); row < fCalibraMode->GetRowMax(2); row++) {
2841 for (Int_t col = fCalibraMode->GetColMin(2); col < fCalibraMode->GetColMax(2); col++) {
2842 fCurrentCoef[1] += (Float_t) fCalROC->GetValue(col,row);
2843 }
2844 }
2845 fCurrentCoef[1] = fCurrentCoef[1] / ((fCalibraMode->GetColMax(2)-fCalibraMode->GetColMin(2))*(fCalibraMode->GetRowMax(2)-fCalibraMode->GetRowMin(2)));
2846 }
2847 return kTRUE;
2848}
2849//_____________________________________________________________________________
2850Bool_t AliTRDCalibraFit::CalculVdriftCoefMean()
2851{
2852 //
2853 // For the detector dect calcul the mean drift velocity for the
2854 // calibration group idect from the choosen database
2855 //
2856
2857 fCurrentCoef[1] = 0.0;
2858 if(fDebugLevel != 1){
2859 if ((fCalibraMode->GetNz(1) > 0) ||
2860 (fCalibraMode->GetNrphi(1) > 0)) {
2861 for (Int_t row = fCalibraMode->GetRowMin(1); row < fCalibraMode->GetRowMax(1); row++) {
2862 for (Int_t col = fCalibraMode->GetColMin(1); col < fCalibraMode->GetColMax(1); col++) {
2863 fCurrentCoef[1] += (Float_t) (fCalROC->GetValue(col,row)*fCalDet->GetValue(fCountDet));
2864 }
2865 }
2866 fCurrentCoef[1] = fCurrentCoef[1] / ((fCalibraMode->GetColMax(1)-fCalibraMode->GetColMin(1))*(fCalibraMode->GetRowMax(1)-fCalibraMode->GetRowMin(1)));
2867 }
2868 else {
2869 //per detectors
2870 fCurrentCoef[1] = (Float_t) fCalDet->GetValue(fCountDet);
2871 }
2872 }
2873 return kTRUE;
2874}
2875//_____________________________________________________________________________
2876Bool_t AliTRDCalibraFit::CalculVdriftLorentzCoef()
2877{
2878 //
2879 // For the detector fCountDet, mean drift velocity and tan lorentzangle
2880 //
2881
2882 fCurrentCoef[1] = fCalDet->GetValue(fCountDet);
2883 fCurrentCoef2[1] = fCalDet2->GetValue(fCountDet);
2884
2885 return kTRUE;
2886}
2887//_____________________________________________________________________________
2888Float_t AliTRDCalibraFit::GetPRFDefault(Int_t plane) const
2889{
2890 //
2891 // Default width of the PRF if there is no database as reference
2892 //
2893 switch(plane)
2894 {
2895 // default database
2896 //case 0: return 0.515;
2897 //case 1: return 0.502;
2898 //case 2: return 0.491;
2899 //case 3: return 0.481;
2900 //case 4: return 0.471;
2901 //case 5: return 0.463;
2902 //default: return 0.0;
2903
2904 // fit database
2905 case 0: return 0.538429;
2906 case 1: return 0.524302;
2907 case 2: return 0.511591;
2908 case 3: return 0.500140;
2909 case 4: return 0.489821;
2910 case 5: return 0.480524;
2911 default: return 0.0;
2912 }
2913}
2914//________________________________________________________________________________
2915void AliTRDCalibraFit::SetCalROC(Int_t i)
2916{
2917 //
2918 // Set the calib object for fCountDet
2919 //
2920
2921 Float_t value = 0.0;
2922
2923 //Get the CalDet object
2924 if(fAccCDB){
2925 AliTRDcalibDB *cal = AliTRDcalibDB::Instance();
2926 if (!cal) {
2927 AliInfo("Could not get calibDB");
2928 return;
2929 }
2930 switch (i)
2931 {
2932 case 0:
2933 if(fCalROC) delete fCalROC;
2934 fCalROC = new AliTRDCalROC(*(cal->GetGainFactorROC(fCountDet)));
2935 break;
2936 case 1:
2937 if(fCalROC) delete fCalROC;
2938 if(fCalROC2) delete fCalROC2;
2939 fCalROC = new AliTRDCalROC(*(cal->GetVdriftROC(fCountDet)));
2940 fCalROC2 = new AliTRDCalROC(*(cal->GetT0ROC(fCountDet)));
2941 break;
2942 case 2:
2943 if(fCalROC) delete fCalROC;
2944 fCalROC = new AliTRDCalROC(*(cal->GetPRFROC(fCountDet))); break;
2945 default: return;
2946 }
2947 }
2948 else{
2949 switch (i)
2950 {
2951 case 0:
2952 if(fCalROC) delete fCalROC;
2953 fCalROC = new AliTRDCalROC(GetPlane(fCountDet),GetChamber(fCountDet));
2954 for(Int_t k = 0; k < fCalROC->GetNchannels(); k++){
2955 fCalROC->SetValue(k,1.0);
2956 }
2957 break;
2958 case 1:
2959 if(fCalROC) delete fCalROC;
2960 if(fCalROC2) delete fCalROC2;
2961 fCalROC = new AliTRDCalROC(GetPlane(fCountDet),GetChamber(fCountDet));
2962 fCalROC2 = new AliTRDCalROC(GetPlane(fCountDet),GetChamber(fCountDet));
2963 for(Int_t k = 0; k < fCalROC->GetNchannels(); k++){
2964 fCalROC->SetValue(k,1.0);
2965 fCalROC2->SetValue(k,0.0);
2966 }
2967 break;
2968 case 2:
2969 if(fCalROC) delete fCalROC;
2970 value = GetPRFDefault(GetPlane(fCountDet));
2971 fCalROC = new AliTRDCalROC(GetPlane(fCountDet),GetChamber(fCountDet));
2972 for(Int_t k = 0; k < fCalROC->GetNchannels(); k++){
2973 fCalROC->SetValue(k,value);
2974 }
2975 break;
2976 default: return;
2977 }
2978 }
2979
2980}
2981//____________Fit Methods______________________________________________________
2982
2983//_____________________________________________________________________________
2984void AliTRDCalibraFit::FitPente(TH1* projPH)
2985{
2986 //
2987 // Slope methode for the drift velocity
2988 //
2989
2990 // Constants
2991 const Float_t kDrWidth = AliTRDgeometry::DrThick();
2992 Int_t binmax = 0;
2993 Int_t binmin = 0;
2994 fPhd[0] = 0.0;
2995 fPhd[1] = 0.0;
2996 fPhd[2] = 0.0;
2997 Int_t ju = 0;
2998 fCurrentCoefE = 0.0;
2999 fCurrentCoefE2 = 0.0;
3000 fCurrentCoef[0] = 0.0;
3001 fCurrentCoef2[0] = 0.0;
3002 TLine *line = new TLine();
3003
3004 // Some variables
3005 TAxis *xpph = projPH->GetXaxis();
3006 Int_t nbins = xpph->GetNbins();
3007 Double_t lowedge = xpph->GetBinLowEdge(1);
3008 Double_t upedge = xpph->GetBinUpEdge(xpph->GetNbins());
3009 Double_t widbins = (upedge - lowedge) / nbins;
3010 Double_t limit = upedge + 0.5 * widbins;
3011 Bool_t put = kTRUE;
3012
3013 // Beginning of the signal
3014 TH1D *pentea = new TH1D("pentea","pentea",projPH->GetNbinsX(),0,(Float_t) limit);
3015 for (Int_t k = 1; k < projPH->GetNbinsX(); k++) {
3016 pentea->SetBinContent(k,(Double_t) (projPH->GetBinContent(k+1) - projPH->GetBinContent(k)));
3017 }
3018 binmax = (Int_t) pentea->GetMaximumBin();
3019 if (binmax <= 1) {
3020 binmax = 2;
3021 AliInfo("Put the binmax from 1 to 2 to enable the fit");
3022 }
3023 if (binmax >= nbins) {
3024 binmax = nbins-1;
3025 put = kFALSE;
3026 AliInfo("Put the binmax from nbins-1 to nbins-2 to enable the fit");
3027 }
3028 pentea->Fit("pol2","0MR","",TMath::Max(pentea->GetBinCenter(binmax-1),0.0),pentea->GetBinCenter(binmax+1));
3029 Float_t l3P1am = pentea->GetFunction("pol2")->GetParameter(1);
3030 Float_t l3P2am = pentea->GetFunction("pol2")->GetParameter(2);
3031 Float_t l3P1amE = pentea->GetFunction("pol2")->GetParError(1);
3032 Float_t l3P2amE = pentea->GetFunction("pol2")->GetParError(2);
3033 if (l3P2am != 0) {
3034 fPhd[0] = -(l3P1am / (2 * l3P2am));
3035 }
3036 if(!fTakeTheMaxPH){
3037 if((l3P1am != 0.0) && (l3P2am != 0.0)){
3038 fCurrentCoefE2 = (l3P1amE/l3P1am + l3P2amE/l3P2am)*fPhd[0];
3039 }
3040 }
3041 // Amplification region
3042 binmax = 0;
3043 ju = 0;
3044 for (Int_t kbin = 1; kbin < projPH->GetNbinsX(); kbin ++) {
3045 if (((projPH->GetBinContent(kbin+1) - projPH->GetBinContent(kbin)) <= 0.0) && (ju == 0) && (kbin > (fPhd[0]/widbins))) {
3046 binmax = kbin;
3047 ju = 1;
3048 }
3049 }
3050 if (binmax <= 1) {
3051 binmax = 2;
3052 AliInfo("Put the binmax from 1 to 2 to enable the fit");
3053 }
3054 if (binmax >= nbins) {
3055 binmax = nbins-1;
3056 put = kFALSE;
3057 AliInfo("Put the binmax from nbins-1 to nbins-2 to enable the fit");
3058 }
3059 projPH->Fit("pol2","0MR","",TMath::Max(projPH->GetBinCenter(binmax-1),0.0),projPH->GetBinCenter(binmax+1));
3060 Float_t l3P1amf = projPH->GetFunction("pol2")->GetParameter(1);
3061 Float_t l3P2amf = projPH->GetFunction("pol2")->GetParameter(2);
3062 Float_t l3P1amfE = projPH->GetFunction("pol2")->GetParError(1);
3063 Float_t l3P2amfE = projPH->GetFunction("pol2")->GetParError(2);
3064 if (l3P2amf != 0) {
3065 fPhd[1] = -(l3P1amf / (2 * l3P2amf));
3066 }
3067 if((l3P1amf != 0.0) && (l3P2amf != 0.0)){
3068 fCurrentCoefE = (l3P1amfE/l3P1amf + l3P2amfE/l3P2amf)*fPhd[1];
3069 }
3070 if(fTakeTheMaxPH){
3071 fCurrentCoefE2 = fCurrentCoefE;
3072 }
3073 // Drift region
3074 TH1D *pente = new TH1D("pente","pente",projPH->GetNbinsX(),0,(Float_t) limit);
3075 for (Int_t k = TMath::Min(binmax+4,projPH->GetNbinsX()); k < projPH->GetNbinsX(); k++) {
3076 pente->SetBinContent(k,(Double_t) (projPH->GetBinContent(k+1) - projPH->GetBinContent(k)));
3077 }
3078 binmin = 0;
3079 if(pente->GetEntries() > 0) binmin = (Int_t) pente->GetMinimumBin();
3080 if (binmin <= 1) {
3081 binmin = 2;
3082 AliInfo("Put the binmax from 1 to 2 to enable the fit");
3083 }
3084 if (binmin >= nbins) {
3085 binmin = nbins-1;
3086 put = kFALSE;
3087 AliInfo("Put the binmax from nbins-1 to nbins-2 to enable the fit");
3088 }
3089 pente->Fit("pol2"
3090 ,"0MR"
3091 ,""
3092 ,TMath::Max(pente->GetBinCenter(binmin-1), 0.0)
3093 ,TMath::Min(pente->GetBinCenter(binmin+1),(Double_t) limit));
3094 Float_t l3P1dr = pente->GetFunction("pol2")->GetParameter(1);
3095 Float_t l3P2dr = pente->GetFunction("pol2")->GetParameter(2);
3096 Float_t l3P1drE = pente->GetFunction("pol2")->GetParError(1);
3097 Float_t l3P2drE = pente->GetFunction("pol2")->GetParError(2);
3098 if (l3P2dr != 0) {
3099 fPhd[2] = -(l3P1dr / (2 * l3P2dr));
3100 }
3101 if((l3P1dr != 0.0) && (l3P2dr != 0.0)){
3102 fCurrentCoefE += (l3P1drE/l3P1dr + l3P2drE/l3P2dr)*fPhd[2];
3103 }
3104 Float_t fPhdt0 = 0.0;
3105 Float_t t0Shift = 0.0;
3106 if(fTakeTheMaxPH) {
3107 fPhdt0 = fPhd[1];
3108 t0Shift = fT0Shift1;
3109 }
3110 else {
3111 fPhdt0 = fPhd[0];
3112 t0Shift = fT0Shift0;
3113 }
3114
3115 if ((fPhd[2] > fPhd[0]) &&
3116 (fPhd[2] > fPhd[1]) &&
3117 (fPhd[1] > fPhd[0]) &&
3118 (put)) {
3119 fCurrentCoef[0] = (kDrWidth) / (fPhd[2]-fPhd[1]);
3120 fNumberFitSuccess++;
3121
3122 if (fPhdt0 >= 0.0) {
3123 fCurrentCoef2[0] = (fPhdt0 - t0Shift) / widbins;
3124 if (fCurrentCoef2[0] < -1.0) {
3125 fCurrentCoef2[0] = fCurrentCoef2[1];
3126 }
3127 }
3128 else {
3129 fCurrentCoef2[0] = fCurrentCoef2[1];
3130 }
3131
3132 }
3133 else {
3134 fCurrentCoef[0] = -TMath::Abs(fCurrentCoef[1]);
3135 fCurrentCoef2[0] = fCurrentCoef2[1];
3136 }
3137
3138 if (fDebugLevel == 1) {
3139 TCanvas *cpentei = new TCanvas("cpentei","cpentei",50,50,600,800);
3140 cpentei->cd();
3141 projPH->Draw();
3142 line->SetLineColor(2);
3143 line->DrawLine(fPhd[0],0,fPhd[0],projPH->GetMaximum());
3144 line->DrawLine(fPhd[1],0,fPhd[1],projPH->GetMaximum());
3145 line->DrawLine(fPhd[2],0,fPhd[2],projPH->GetMaximum());
3146 AliInfo(Form("fPhd[0] (beginning of the signal): %f" ,(Float_t) fPhd[0]));
3147 AliInfo(Form("fPhd[1] (end of the amplification region): %f" ,(Float_t) fPhd[1]));
3148 AliInfo(Form("fPhd[2] (end of the drift region): %f" ,(Float_t) fPhd[2]));
3149 AliInfo(Form("fVriftCoef[1] (with only the drift region(default)): %f",(Float_t) fCurrentCoef[0]));
3150 TCanvas *cpentei2 = new TCanvas("cpentei2","cpentei2",50,50,600,800);
3151 cpentei2->cd();
3152 pentea->Draw();
3153 TCanvas *cpentei3 = new TCanvas("cpentei3","cpentei3",50,50,600,800);
3154 cpentei3->cd();
3155 pente->Draw();
3156 }
3157 else {
3158 delete pentea;
3159 delete pente;
3160 }
3161}
3162//_____________________________________________________________________________
3163void AliTRDCalibraFit::FitLagrangePoly(TH1* projPH)
3164{
3165 //
3166 // Slope methode but with polynomes de Lagrange
3167 //
3168
3169 // Constants
3170 const Float_t kDrWidth = AliTRDgeometry::DrThick();
3171 Int_t binmax = 0;
3172 Int_t binmin = 0;
3173 Double_t *x = new Double_t[5];
3174 Double_t *y = new Double_t[5];
3175 x[0] = 0.0;
3176 x[1] = 0.0;
3177 x[2] = 0.0;
3178 x[3] = 0.0;
3179 x[4] = 0.0;
3180 y[0] = 0.0;
3181 y[1] = 0.0;
3182 y[2] = 0.0;
3183 y[3] = 0.0;
3184 y[4] = 0.0;
3185 fPhd[0] = 0.0;
3186 fPhd[1] = 0.0;
3187 fPhd[2] = 0.0;
3188 Int_t ju = 0;
3189 fCurrentCoefE = 0.0;
3190 fCurrentCoefE2 = 1.0;
3191 fCurrentCoef[0] = 0.0;
3192 fCurrentCoef2[0] = 0.0;
3193 TLine *line = new TLine();
3194 TF1 * polynome = 0x0;
3195 TF1 * polynomea = 0x0;
3196 TF1 * polynomeb = 0x0;
3197 Double_t *c = 0x0;
3198
3199 // Some variables
3200 TAxis *xpph = projPH->GetXaxis();
3201 Int_t nbins = xpph->GetNbins();
3202 Double_t lowedge = xpph->GetBinLowEdge(1);
3203 Double_t upedge = xpph->GetBinUpEdge(xpph->GetNbins());
3204 Double_t widbins = (upedge - lowedge) / nbins;
3205 Double_t limit = upedge + 0.5 * widbins;
3206
3207
3208 Bool_t put = kTRUE;
3209
3210 // Beginning of the signal
3211 TH1D *pentea = new TH1D("pentea","pentea",projPH->GetNbinsX(),0,(Float_t) limit);
3212 for (Int_t k = 1; k < projPH->GetNbinsX(); k++) {
3213 pentea->SetBinContent(k,(Double_t) (projPH->GetBinContent(k+1) - projPH->GetBinContent(k)));
3214 }
3215
3216 binmax = (Int_t) pentea->GetMaximumBin();
3217
3218 Double_t minnn = 0.0;
3219 Double_t maxxx = 0.0;
3220
3221 Int_t kase = nbins-binmax;
3222
3223 switch(kase)
3224 {
3225 case 0:
3226 put = kFALSE;
3227 break;
3228 case 1:
3229 minnn = pentea->GetBinCenter(binmax-2);
3230 maxxx = pentea->GetBinCenter(binmax);
3231 x[0] = pentea->GetBinCenter(binmax-2);
3232 x[1] = pentea->GetBinCenter(binmax-1);
3233 x[2] = pentea->GetBinCenter(binmax);
3234 y[0] = pentea->GetBinContent(binmax-2);
3235 y[1] = pentea->GetBinContent(binmax-1);
3236 y[2] = pentea->GetBinContent(binmax);
3237 c = CalculPolynomeLagrange2(x,y);
3238 AliInfo("At the limit for beginning!");
3239 break;
3240 case 2:
3241 minnn = pentea->GetBinCenter(binmax-2);
3242 maxxx = pentea->GetBinCenter(binmax+1);
3243 x[0] = pentea->GetBinCenter(binmax-2);
3244 x[1] = pentea->GetBinCenter(binmax-1);
3245 x[2] = pentea->GetBinCenter(binmax);
3246 x[3] = pentea->GetBinCenter(binmax+1);
3247 y[0] = pentea->GetBinContent(binmax-2);
3248 y[1] = pentea->GetBinContent(binmax-1);
3249 y[2] = pentea->GetBinContent(binmax);
3250 y[3] = pentea->GetBinContent(binmax+1);
3251 c = CalculPolynomeLagrange3(x,y);
3252 break;
3253 default:
3254 switch(binmax){
3255 case 0:
3256 put = kFALSE;
3257 break;
3258 case 1:
3259 minnn = pentea->GetBinCenter(binmax);
3260 maxxx = pentea->GetBinCenter(binmax+2);
3261 x[0] = pentea->GetBinCenter(binmax);
3262 x[1] = pentea->GetBinCenter(binmax+1);
3263 x[2] = pentea->GetBinCenter(binmax+2);
3264 y[0] = pentea->GetBinContent(binmax);
3265 y[1] = pentea->GetBinContent(binmax+1);
3266 y[2] = pentea->GetBinContent(binmax+2);
3267 c = CalculPolynomeLagrange2(x,y);
3268 break;
3269 case 2:
3270 minnn = pentea->GetBinCenter(binmax-1);
3271 maxxx = pentea->GetBinCenter(binmax+2);
3272 x[0] = pentea->GetBinCenter(binmax-1);
3273 x[1] = pentea->GetBinCenter(binmax);
3274 x[2] = pentea->GetBinCenter(binmax+1);
3275 x[3] = pentea->GetBinCenter(binmax+2);
3276 y[0] = pentea->GetBinContent(binmax-1);
3277 y[1] = pentea->GetBinContent(binmax);
3278 y[2] = pentea->GetBinContent(binmax+1);
3279 y[3] = pentea->GetBinContent(binmax+2);
3280 c = CalculPolynomeLagrange3(x,y);
3281 break;
3282 default:
3283 minnn = pentea->GetBinCenter(binmax-2);
3284 maxxx = pentea->GetBinCenter(binmax+2);
3285 x[0] = pentea->GetBinCenter(binmax-2);
3286 x[1] = pentea->GetBinCenter(binmax-1);
3287 x[2] = pentea->GetBinCenter(binmax);
3288 x[3] = pentea->GetBinCenter(binmax+1);
3289 x[4] = pentea->GetBinCenter(binmax+2);
3290 y[0] = pentea->GetBinContent(binmax-2);
3291 y[1] = pentea->GetBinContent(binmax-1);
3292 y[2] = pentea->GetBinContent(binmax);
3293 y[3] = pentea->GetBinContent(binmax+1);
3294 y[4] = pentea->GetBinContent(binmax+2);
3295 c = CalculPolynomeLagrange4(x,y);
3296 break;
3297 }
3298 break;
3299 }
3300
3301
3302 if(put) {
3303 polynomeb = new TF1("polb","[0]+[1]*x+[2]*x*x+[3]*x*x*x+[4]*x*x*x*x",minnn,maxxx);
3304 polynomeb->SetParameters(c[0],c[1],c[2],c[3],c[4]);
3305
3306 Double_t step = (maxxx-minnn)/10000;
3307 Double_t l = minnn;
3308 Double_t maxvalue = 0.0;
3309 Double_t placemaximum = minnn;
3310 for(Int_t o = 0; o < 10000; o++){
3311 if(o == 0) maxvalue = polynomeb->Eval(l);
3312 if(maxvalue < (polynomeb->Eval(l))){
3313 maxvalue = polynomeb->Eval(l);
3314 placemaximum = l;
3315 }
3316 l += step;
3317 }
3318 fPhd[0] = placemaximum;
3319 }
3320
3321 // Amplification region
3322 binmax = 0;
3323 ju = 0;
3324 for (Int_t kbin = 1; kbin < projPH->GetNbinsX(); kbin ++) {
3325 if (((projPH->GetBinContent(kbin+1) - projPH->GetBinContent(kbin)) <= 0.0) && (ju == 0) && (kbin > (fPhd[0]/widbins))) {
3326 binmax = kbin;
3327 ju = 1;
3328 }
3329 }
3330
3331 Double_t minn = 0.0;
3332 Double_t maxx = 0.0;
3333 x[0] = 0.0;
3334 x[1] = 0.0;
3335 x[2] = 0.0;
3336 x[3] = 0.0;
3337 x[4] = 0.0;
3338 y[0] = 0.0;
3339 y[1] = 0.0;
3340 y[2] = 0.0;
3341 y[3] = 0.0;
3342 y[4] = 0.0;
3343
3344 Int_t kase1 = nbins - binmax;
3345
3346 //Determination of minn and maxx
3347 //case binmax = nbins
3348 //pol2
3349 switch(kase1)
3350 {
3351 case 0:
3352 minn = projPH->GetBinCenter(binmax-2);
3353 maxx = projPH->GetBinCenter(binmax);
3354 x[0] = projPH->GetBinCenter(binmax-2);
3355 x[1] = projPH->GetBinCenter(binmax-1);
3356 x[2] = projPH->GetBinCenter(binmax);
3357 y[0] = projPH->GetBinContent(binmax-2);
3358 y[1] = projPH->GetBinContent(binmax-1);
3359 y[2] = projPH->GetBinContent(binmax);
3360 c = CalculPolynomeLagrange2(x,y);
3361 //AliInfo("At the limit for the drift!");
3362 break;
3363 case 1:
3364 minn = projPH->GetBinCenter(binmax-2);
3365 maxx = projPH->GetBinCenter(binmax+1);
3366 x[0] = projPH->GetBinCenter(binmax-2);
3367 x[1] = projPH->GetBinCenter(binmax-1);
3368 x[2] = projPH->GetBinCenter(binmax);
3369 x[3] = projPH->GetBinCenter(binmax+1);
3370 y[0] = projPH->GetBinContent(binmax-2);
3371 y[1] = projPH->GetBinContent(binmax-1);
3372 y[2] = projPH->GetBinContent(binmax);
3373 y[3] = projPH->GetBinContent(binmax+1);
3374 c = CalculPolynomeLagrange3(x,y);
3375 break;
3376 default:
3377 switch(binmax)
3378 {
3379 case 0:
3380 put = kFALSE;
3381 break;
3382 case 1:
3383 minn = projPH->GetBinCenter(binmax);
3384 maxx = projPH->GetBinCenter(binmax+2);
3385 x[0] = projPH->GetBinCenter(binmax);
3386 x[1] = projPH->GetBinCenter(binmax+1);
3387 x[2] = projPH->GetBinCenter(binmax+2);
3388 y[0] = projPH->GetBinContent(binmax);
3389 y[1] = projPH->GetBinContent(binmax+1);
3390 y[2] = projPH->GetBinContent(binmax+2);
3391 c = CalculPolynomeLagrange2(x,y);
3392 break;
3393 case 2:
3394 minn = projPH->GetBinCenter(binmax-1);
3395 maxx = projPH->GetBinCenter(binmax+2);
3396 x[0] = projPH->GetBinCenter(binmax-1);
3397 x[1] = projPH->GetBinCenter(binmax);
3398 x[2] = projPH->GetBinCenter(binmax+1);
3399 x[3] = projPH->GetBinCenter(binmax+2);
3400 y[0] = projPH->GetBinContent(binmax-1);
3401 y[1] = projPH->GetBinContent(binmax);
3402 y[2] = projPH->GetBinContent(binmax+1);
3403 y[3] = projPH->GetBinContent(binmax+2);
3404 c = CalculPolynomeLagrange3(x,y);
3405 break;
3406 default:
3407 minn = projPH->GetBinCenter(binmax-2);
3408 maxx = projPH->GetBinCenter(binmax+2);
3409 x[0] = projPH->GetBinCenter(binmax-2);
3410 x[1] = projPH->GetBinCenter(binmax-1);
3411 x[2] = projPH->GetBinCenter(binmax);
3412 x[3] = projPH->GetBinCenter(binmax+1);
3413 x[4] = projPH->GetBinCenter(binmax+2);
3414 y[0] = projPH->GetBinContent(binmax-2);
3415 y[1] = projPH->GetBinContent(binmax-1);
3416 y[2] = projPH->GetBinContent(binmax);
3417 y[3] = projPH->GetBinContent(binmax+1);
3418 y[4] = projPH->GetBinContent(binmax+2);
3419 c = CalculPolynomeLagrange4(x,y);
3420 break;
3421 }
3422 break;
3423 }
3424
3425 if(put) {
3426 polynomea = new TF1("pola","[0]+[1]*x+[2]*x*x+[3]*x*x*x+[4]*x*x*x*x",minn,maxx);
3427 polynomea->SetParameters(c[0],c[1],c[2],c[3],c[4]);
3428
3429 Double_t step = (maxx-minn)/1000;
3430 Double_t l = minn;
3431 Double_t maxvalue = 0.0;
3432 Double_t placemaximum = minn;
3433 for(Int_t o = 0; o < 1000; o++){
3434 if(o == 0) maxvalue = polynomea->Eval(l);
3435 if(maxvalue < (polynomea->Eval(l))){
3436 maxvalue = polynomea->Eval(l);
3437 placemaximum = l;
3438 }
3439 l += step;
3440 }
3441 fPhd[1] = placemaximum;
3442 }
3443
3444 // Drift region
3445 TH1D *pente = new TH1D("pente","pente", projPH->GetNbinsX(),0,(Float_t) limit);
3446 for (Int_t k = TMath::Min(binmax+4, projPH->GetNbinsX()); k < projPH->GetNbinsX(); k++) {
3447 pente->SetBinContent(k,(Double_t) (projPH->GetBinContent(k+1) - projPH->GetBinContent(k)));
3448 }
3449 binmin = 0;
3450 if(pente->GetEntries() > 0) binmin = (Int_t) pente->GetMinimumBin();
3451
3452 //should not happen
3453 if (binmin <= 1) {
3454 binmin = 2;
3455 put = 1;
3456 AliInfo("Put the binmax from 1 to 2 to enable the fit");
3457 }
3458
3459 //check
3460 if((projPH->GetBinContent(binmin)-projPH->GetBinError(binmin)) < (projPH->GetBinContent(binmin+1))) put = kFALSE;
3461 if((projPH->GetBinContent(binmin)+projPH->GetBinError(binmin)) > (projPH->GetBinContent(binmin-1))) put = kFALSE;
3462
3463 x[0] = 0.0;
3464 x[1] = 0.0;
3465 x[2] = 0.0;
3466 x[3] = 0.0;
3467 x[4] = 0.0;
3468 y[0] = 0.0;
3469 y[1] = 0.0;
3470 y[2] = 0.0;
3471 y[3] = 0.0;
3472 y[4] = 0.0;
3473 Double_t min = 0.0;
3474 Double_t max = 0.0;
3475 Bool_t case1 = kFALSE;
3476 Bool_t case2 = kFALSE;
3477 Bool_t case4 = kFALSE;
3478
3479 //Determination of min and max
3480 //case binmin <= nbins-3
3481 //pol4 case 3
3482 if((binmin <= (nbins-3)) && ((binmin-2) >= TMath::Min(binmax+4, projPH->GetNbinsX()))){
3483 min = pente->GetBinCenter(binmin-2);
3484 max = pente->GetBinCenter(binmin+2);
3485 x[0] = pente->GetBinCenter(binmin-2);
3486 x[1] = pente->GetBinCenter(binmin-1);
3487 x[2] = pente->GetBinCenter(binmin);
3488 x[3] = pente->GetBinCenter(binmin+1);
3489 x[4] = pente->GetBinCenter(binmin+2);
3490 y[0] = pente->GetBinContent(binmin-2);
3491 y[1] = pente->GetBinContent(binmin-1);
3492 y[2] = pente->GetBinContent(binmin);
3493 y[3] = pente->GetBinContent(binmin+1);
3494 y[4] = pente->GetBinContent(binmin+2);
3495 //Calcul the polynome de Lagrange
3496 c = CalculPolynomeLagrange4(x,y);
3497 //richtung +/-
3498 if((pente->GetBinContent(binmin+2) <= pente->GetBinContent(binmin+1)) &&
3499 (pente->GetBinContent(binmin-2) <= pente->GetBinContent(binmin-1))) put = kFALSE;
3500 if(((binmin+3) <= (nbins-1)) &&
3501 (pente->GetBinContent(binmin+3) <= pente->GetBinContent(binmin+2)) &&
3502 ((binmin-3) >= TMath::Min(binmax+4, projPH->GetNbinsX())) &&
3503 (pente->GetBinContent(binmin-3) <= pente->GetBinContent(binmin-2))) put = kFALSE;
3504 if((pente->GetBinContent(binmin+2) <= pente->GetBinContent(binmin+1)) &&
3505 (pente->GetBinContent(binmin-2) > pente->GetBinContent(binmin-1))) case1 = kTRUE;
3506 if((pente->GetBinContent(binmin+2) > pente->GetBinContent(binmin+1)) &&
3507 (pente->GetBinContent(binmin-2) <= pente->GetBinContent(binmin-1))) case4 = kTRUE;
3508 }
3509 //case binmin = nbins-2
3510 //pol3 case 1
3511 if(((binmin == (nbins-2)) && ((binmin-2) >= TMath::Min(binmax+4, projPH->GetNbinsX()))) ||
3512 (case1)){
3513 min = pente->GetBinCenter(binmin-2);
3514 max = pente->GetBinCenter(binmin+1);
3515 x[0] = pente->GetBinCenter(binmin-2);
3516 x[1] = pente->GetBinCenter(binmin-1);
3517 x[2] = pente->GetBinCenter(binmin);
3518 x[3] = pente->GetBinCenter(binmin+1);
3519 y[0] = pente->GetBinContent(binmin-2);
3520 y[1] = pente->GetBinContent(binmin-1);
3521 y[2] = pente->GetBinContent(binmin);
3522 y[3] = pente->GetBinContent(binmin+1);
3523 //Calcul the polynome de Lagrange
3524 c = CalculPolynomeLagrange3(x,y);
3525 //richtung +: nothing
3526 //richtung -
3527 if((pente->GetBinContent(binmin-2) <= pente->GetBinContent(binmin-1))) case2 = kTRUE;
3528 }
3529 //pol3 case 4
3530 if(((binmin <= (nbins-3)) && ((binmin-1) == TMath::Min(binmax+4, projPH->GetNbinsX()))) ||
3531 (case4)){
3532 min = pente->GetBinCenter(binmin-1);
3533 max = pente->GetBinCenter(binmin+2);
3534 x[0] = pente->GetBinCenter(binmin-1);
3535 x[1] = pente->GetBinCenter(binmin);
3536 x[2] = pente->GetBinCenter(binmin+1);
3537 x[3] = pente->GetBinCenter(binmin+2);
3538 y[0] = pente->GetBinContent(binmin-1);
3539 y[1] = pente->GetBinContent(binmin);
3540 y[2] = pente->GetBinContent(binmin+1);
3541 y[3] = pente->GetBinContent(binmin+2);
3542 //Calcul the polynome de Lagrange
3543 c = CalculPolynomeLagrange3(x,y);
3544 //richtung +
3545 if((pente->GetBinContent(binmin+2) <= pente->GetBinContent(binmin+1))) case2 = kTRUE;
3546 }
3547 //pol2 case 5
3548 if((binmin <= (nbins-3)) && (binmin == TMath::Min(binmax+4, projPH->GetNbinsX()))){
3549 min = pente->GetBinCenter(binmin);
3550 max = pente->GetBinCenter(binmin+2);
3551 x[0] = pente->GetBinCenter(binmin);
3552 x[1] = pente->GetBinCenter(binmin+1);
3553 x[2] = pente->GetBinCenter(binmin+2);
3554 y[0] = pente->GetBinContent(binmin);
3555 y[1] = pente->GetBinContent(binmin+1);
3556 y[2] = pente->GetBinContent(binmin+2);
3557 //Calcul the polynome de Lagrange
3558 c = CalculPolynomeLagrange2(x,y);
3559 //richtung +
3560 if((pente->GetBinContent(binmin+2) <= pente->GetBinContent(binmin+1))) put = kFALSE;
3561 }
3562 //pol2 case 2
3563 if(((binmin == (nbins-2)) && ((binmin-1) == TMath::Min(binmax+4, projPH->GetNbinsX()))) ||
3564 (case2)){
3565 min = pente->GetBinCenter(binmin-1);
3566 max = pente->GetBinCenter(binmin+1);
3567 x[0] = pente->GetBinCenter(binmin-1);
3568 x[1] = pente->GetBinCenter(binmin);
3569 x[2] = pente->GetBinCenter(binmin+1);
3570 y[0] = pente->GetBinContent(binmin-1);
3571 y[1] = pente->GetBinContent(binmin);
3572 y[2] = pente->GetBinContent(binmin+1);
3573 //Calcul the polynome de Lagrange
3574 c = CalculPolynomeLagrange2(x,y);
3575 //richtung +: nothing
3576 //richtung -: nothing
3577 }
3578 //case binmin = nbins-1
3579 //pol2 case 0
3580 if((binmin == (nbins-1)) && ((binmin-2) >= TMath::Min(binmax+4, projPH->GetNbinsX()))){
3581 min = pente->GetBinCenter(binmin-2);
3582 max = pente->GetBinCenter(binmin);
3583 x[0] = pente->GetBinCenter(binmin-2);
3584 x[1] = pente->GetBinCenter(binmin-1);
3585 x[2] = pente->GetBinCenter(binmin);
3586 y[0] = pente->GetBinContent(binmin-2);
3587 y[1] = pente->GetBinContent(binmin-1);
3588 y[2] = pente->GetBinContent(binmin);
3589 //Calcul the polynome de Lagrange
3590 c = CalculPolynomeLagrange2(x,y);
3591 //AliInfo("At the limit for the drift!");
3592 //fluctuation too big!
3593 //richtung +: nothing
3594 //richtung -
3595 if((pente->GetBinContent(binmin-2) <= pente->GetBinContent(binmin-1))) put = kFALSE;
3596 }
3597 if((binmin == (nbins-1)) && ((binmin-2) < TMath::Min(binmax+4, projPH->GetNbinsX()))) {
3598 put = kFALSE;
3599 AliInfo("At the limit for the drift and not usable!");
3600 }
3601
3602 //pass
3603 if((binmin == (nbins-2)) && ((binmin-1) < TMath::Min(binmax+4, projPH->GetNbinsX()))){
3604 put = kFALSE;
3605 AliInfo("For the drift...problem!");
3606 }
3607 //pass but should not happen
3608 if((binmin <= (nbins-3)) && (binmin < TMath::Min(binmax+4, projPH->GetNbinsX()))){
3609 put = kFALSE;
3610 AliInfo("For the drift...problem!");
3611 }
3612
3613 if(put) {
3614 polynome = new TF1("pol","[0]+[1]*x+[2]*x*x+[3]*x*x*x+[4]*x*x*x*x",min,max);
3615 polynome->SetParameters(c[0],c[1],c[2],c[3],c[4]);
3616 //AliInfo(Form("GetMinimum of the function %f",polynome->GetMinimumX()));
3617 Double_t step = (max-min)/1000;
3618 Double_t l = min;
3619 Double_t minvalue = 0.0;
3620 Double_t placeminimum = min;
3621 for(Int_t o = 0; o < 1000; o++){
3622 if(o == 0) minvalue = polynome->Eval(l);
3623 if(minvalue > (polynome->Eval(l))){
3624 minvalue = polynome->Eval(l);
3625 placeminimum = l;
3626 }
3627 l += step;
3628 }
3629 fPhd[2] = placeminimum;
3630 }
3631
3632 Float_t fPhdt0 = 0.0;
3633 Float_t t0Shift = 0.0;
3634 if(fTakeTheMaxPH) {
3635 fPhdt0 = fPhd[1];
3636 t0Shift = fT0Shift1;
3637 }
3638 else {
3639 fPhdt0 = fPhd[0];
3640 t0Shift = fT0Shift0;
3641 }
3642
3643 if ((fPhd[2] > fPhd[0]) &&
3644 (fPhd[2] > fPhd[1]) &&
3645 (fPhd[1] > fPhd[0]) &&
3646 (put)) {
3647 fCurrentCoef[0] = (kDrWidth) / (fPhd[2]-fPhd[1]);
3648 fNumberFitSuccess++;
3649 if (fPhdt0 >= 0.0) {
3650 fCurrentCoef2[0] = (fPhdt0 - t0Shift) / widbins;
3651 if (fCurrentCoef2[0] < -1.0) {
3652 fCurrentCoef2[0] = fCurrentCoef2[1];
3653 }
3654 }
3655 else {
3656 fCurrentCoef2[0] = fCurrentCoef2[1];
3657 }
3658 }
3659 else {
3660 fCurrentCoef[0] = -TMath::Abs(fCurrentCoef[1]);
3661 fCurrentCoef2[0] = fCurrentCoef2[1];
3662 //printf("Fit failed!\n");
3663 }
3664
3665 if (fDebugLevel == 1) {
3666 TCanvas *cpentei = new TCanvas("cpentei","cpentei",50,50,600,800);
3667 cpentei->cd();
3668 projPH->Draw();
3669 line->SetLineColor(2);
3670 line->DrawLine(fPhd[0],0,fPhd[0],projPH->GetMaximum());
3671 line->DrawLine(fPhd[1],0,fPhd[1],projPH->GetMaximum());
3672 line->DrawLine(fPhd[2],0,fPhd[2],projPH->GetMaximum());
3673 AliInfo(Form("fPhd[0] (beginning of the signal): %f" ,(Float_t) fPhd[0]));
3674 AliInfo(Form("fPhd[1] (end of the amplification region): %f" ,(Float_t) fPhd[1]));
3675 AliInfo(Form("fPhd[2] (end of the drift region): %f" ,(Float_t) fPhd[2]));
3676 AliInfo(Form("fVriftCoef[3] (with only the drift region(default)): %f",(Float_t) fCurrentCoef[0]));
3677 TCanvas *cpentei2 = new TCanvas("cpentei2","cpentei2",50,50,600,800);
3678 cpentei2->cd();
3679 pentea->Draw();
3680 TCanvas *cpentei3 = new TCanvas("cpentei3","cpentei3",50,50,600,800);
3681 cpentei3->cd();
3682 pente->Draw();
3683 }
3684 else {
3685 delete pentea;
3686 delete pente;
3687 delete polynome;
3688 delete polynomea;
3689 delete polynomeb;
3690 }
3691
3692 projPH->SetDirectory(0);
3693
3694}
3695
3696//_____________________________________________________________________________
3697void AliTRDCalibraFit::FitPH(TH1* projPH, Int_t idect)
3698{
3699 //
3700 // Fit methode for the drift velocity
3701 //
3702
3703 // Constants
3704 const Float_t kDrWidth = AliTRDgeometry::DrThick();
3705
3706 // Some variables
3707 TAxis *xpph = projPH->GetXaxis();
3708 Double_t upedge = xpph->GetBinUpEdge(xpph->GetNbins());
3709
3710 TF1 *fPH = new TF1("fPH",AliTRDCalibraFit::PH,-0.05,3.2,6);
3711 fPH->SetParameter(0,0.469); // Scaling
3712 fPH->SetParameter(1,0.18); // Start
3713 fPH->SetParameter(2,0.0857325); // AR
3714 fPH->SetParameter(3,1.89); // DR
3715 fPH->SetParameter(4,0.08); // QA/QD
3716 fPH->SetParameter(5,0.0); // Baseline
3717
3718 TLine *line = new TLine();
3719
3720 fCurrentCoef[0] = 0.0;
3721 fCurrentCoef2[0] = 0.0;
3722 fCurrentCoefE = 0.0;
3723 fCurrentCoefE2 = 0.0;
3724
3725 if (idect%fFitPHPeriode == 0) {
3726
3727 AliInfo(Form("The detector %d will be fitted",idect));
3728 fPH->SetParameter(0,(projPH->Integral()-(projPH->GetBinContent(1)*projPH->GetNbinsX())) * 0.00028); // Scaling
3729 fPH->SetParameter(1,fPhd[0] - 0.1); // Start
3730 fPH->SetParameter(2,fPhd[1] - fPhd[0]); // AR
3731 fPH->SetParameter(3,fPhd[2] - fPhd[1]); // DR
3732 fPH->SetParameter(4,0.225); // QA/QD
3733 fPH->SetParameter(5,(Float_t) projPH->GetBinContent(1));
3734
3735 if (fDebugLevel != 1) {
3736 projPH->Fit(fPH,"0M","",0.0,upedge);
3737 }
3738 else {
3739 TCanvas *cpente = new TCanvas("cpente","cpente",50,50,600,800);
3740 cpente->cd();
3741 projPH->Fit(fPH,"M+","",0.0,upedge);
3742 projPH->Draw("E0");
3743 line->SetLineColor(4);
3744 line->DrawLine(fPH->GetParameter(1)
3745 ,0
3746 ,fPH->GetParameter(1)
3747 ,projPH->GetMaximum());
3748 line->DrawLine(fPH->GetParameter(1)+fPH->GetParameter(2)
3749 ,0
3750 ,fPH->GetParameter(1)+fPH->GetParameter(2)
3751 ,projPH->GetMaximum());
3752 line->DrawLine(fPH->GetParameter(1)+fPH->GetParameter(2)+fPH->GetParameter(3)
3753 ,0
3754 ,fPH->GetParameter(1)+fPH->GetParameter(2)+fPH->GetParameter(3)
3755 ,projPH->GetMaximum());
3756 }
3757
3758 if (fPH->GetParameter(3) != 0) {
3759 fNumberFitSuccess++;
3760 fCurrentCoef[0] = kDrWidth / (fPH->GetParameter(3));
3761 fCurrentCoefE = (fPH->GetParError(3)/fPH->GetParameter(3))*fCurrentCoef[0];
3762 fCurrentCoef2[0] = fPH->GetParameter(1);
3763 fCurrentCoefE2 = fPH->GetParError(1);
3764 }
3765 else {
3766 fCurrentCoef[0] = -TMath::Abs(fCurrentCoef[1]);
3767 fCurrentCoef2[0] = fCurrentCoef2[1];
3768 }
3769
3770 }
3771 else {
3772
3773 // Put the default value
3774 fCurrentCoef[0] = -TMath::Abs(fCurrentCoef[1]);
3775 fCurrentCoef2[0] = fCurrentCoef2[1];
3776 }
3777
3778 if (fDebugLevel != 1) {
3779 delete fPH;
3780 }
3781
3782}
3783//_____________________________________________________________________________
3784Bool_t AliTRDCalibraFit::FitPRFGausMI(Double_t *arraye, Double_t *arraym, Double_t *arrayme, Int_t nBins, Float_t xMin, Float_t xMax)
3785{
3786 //
3787 // Fit methode for the sigma of the pad response function
3788 //
3789
3790 TVectorD param(3);
3791
3792 fCurrentCoef[0] = 0.0;
3793 fCurrentCoefE = 0.0;
3794
3795 Double_t ret = FitGausMI(arraye, arraym, arrayme, nBins, xMin, xMax,&param);
3796
3797 if(ret == -4){
3798 fCurrentCoef[0] = -fCurrentCoef[1];
3799 return kFALSE;
3800 }
3801 else {
3802 fNumberFitSuccess++;
3803 fCurrentCoef[0] = param[2];
3804 fCurrentCoefE = ret;
3805 return kTRUE;
3806 }
3807}
3808//_____________________________________________________________________________
3809Double_t AliTRDCalibraFit::FitGausMI(Double_t *arraye, Double_t *arraym, Double_t *arrayme, Int_t nBins, Float_t xMin, Float_t xMax, TVectorD *param, Bool_t kError)
3810{
3811 //
3812 // Fit methode for the sigma of the pad response function
3813 //
3814
3815 //We should have at least 3 points
3816 if(nBins <=3) return -4.0;
3817
3818 TLinearFitter fitter(3,"pol2");
3819 fitter.StoreData(kFALSE);
3820 fitter.ClearPoints();
3821 TVectorD par(3);
3822 Float_t binWidth = (xMax-xMin)/(Float_t)nBins;
3823 Float_t entries = 0;
3824 Int_t nbbinwithentries = 0;
3825 for (Int_t i=0; i<nBins; i++){
3826 entries+=arraye[i];
3827 if(arraye[i] > 15) nbbinwithentries++;
3828 //printf("entries for i %d: %f\n",i,arraye[i]);
3829 }
3830 if ((entries<700) || (nbbinwithentries < ((Int_t)(nBins/2)))) return -4;
3831 //printf("entries %f\n",entries);
3832 //printf("nbbinwithentries %d\n",nbbinwithentries);
3833
3834 Int_t npoints=0;
3835 Float_t errorm = 0.0;
3836 Float_t errorn = 0.0;
3837 Float_t error = 0.0;
3838
3839 //
3840 for (Int_t ibin=0;ibin<nBins; ibin++){
3841 Float_t entriesI = arraye[ibin];
3842 Float_t valueI = arraym[ibin];
3843 Double_t xcenter = 0.0;
3844 Float_t val = 0.0;
3845 if ((entriesI>15) && (valueI>0.0)){
3846 xcenter = xMin+(ibin+0.5)*binWidth;
3847 errorm = 0.0;
3848 errorn = 0.0;
3849 error = 0.0;
3850 if(!kError){
3851 if((valueI + 0.01) > 0.0) errorm = TMath::Log((valueI + 0.01)/valueI);
3852 if((valueI - 0.01) > 0.0) errorn = TMath::Log((valueI - 0.01)/valueI);
3853 error = TMath::Max(TMath::Abs(errorm),TMath::Abs(errorn));
3854 }
3855 else{
3856 if((valueI + arrayme[ibin]) > 0.0) errorm = TMath::Log((valueI + arrayme[ibin])/valueI);
3857 if((valueI - arrayme[ibin]) > 0.0) errorn = TMath::Log((valueI - arrayme[ibin])/valueI);
3858 error = TMath::Max(TMath::Abs(errorm),TMath::Abs(errorn));
3859 }
3860 if(error == 0.0) continue;
3861 val = TMath::Log(Float_t(valueI));
3862 fitter.AddPoint(&xcenter,val,error);
3863 npoints++;
3864 }
3865
3866 if(fDebugLevel > 1){
3867
3868 if ( !fDebugStreamer ) {
3869 //debug stream
3870 TDirectory *backup = gDirectory;
3871 fDebugStreamer = new TTreeSRedirector("TRDDebugFit.root");
3872 if ( backup ) backup->cd(); //we don't want to be cd'd to the debug streamer
3873 }
3874
3875 Int_t detector = fCountDet;
3876 Int_t plane = GetPlane(fCountDet);
3877 Int_t group = ibin;
3878
3879 (* fDebugStreamer) << "FitGausMIFill"<<
3880 "detector="<<detector<<
3881 "plane="<<plane<<
3882 "nbins="<<nBins<<
3883 "group="<<group<<
3884 "entriesI="<<entriesI<<
3885 "valueI="<<valueI<<
3886 "val="<<val<<
3887 "xcenter="<<xcenter<<
3888 "errorm="<<errorm<<
3889 "errorn="<<errorn<<
3890 "error="<<error<<
3891 "kError="<<kError<<
3892 "\n";
3893 }
3894
3895 }
3896
3897 if(npoints <=3) return -4.0;
3898
3899 Double_t chi2 = 0;
3900 if (npoints>3){
3901 fitter.Eval();
3902 fitter.GetParameters(par);
3903 chi2 = fitter.GetChisquare()/Float_t(npoints);
3904
3905
3906 if (!param) param = new TVectorD(3);
3907 if(par[2] == 0.0) return -4.0;
3908 Double_t x = TMath::Sqrt(TMath::Abs(-2*par[2]));
3909 Double_t deltax = (fitter.GetParError(2))/x;
3910 Double_t errorparam2 = TMath::Abs(deltax)/(x*x);
3911 chi2 = errorparam2;
3912
3913 (*param)[1] = par[1]/(-2.*par[2]);
3914 (*param)[2] = 1./TMath::Sqrt(TMath::Abs(-2.*par[2]));
3915 Double_t lnparam0 = par[0]+ par[1]* (*param)[1] + par[2]*(*param)[1]*(*param)[1];
3916 if ( lnparam0>307 ) return -4;
3917 (*param)[0] = TMath::Exp(lnparam0);
3918
3919 if(fDebugLevel > 1){
3920
3921 if ( !fDebugStreamer ) {
3922 //debug stream
3923 TDirectory *backup = gDirectory;
3924 fDebugStreamer = new TTreeSRedirector("TRDDebugFit.root");
3925 if ( backup ) backup->cd(); //we don't want to be cd'd to the debug streamer
3926 }
3927
3928 Int_t detector = fCountDet;
3929 Int_t plane = GetPlane(fCountDet);
3930
3931
3932 (* fDebugStreamer) << "FitGausMIFit"<<
3933 "detector="<<detector<<
3934 "plane="<<plane<<
3935 "nbins="<<nBins<<
3936 "errorsigma="<<chi2<<
3937 "mean="<<(*param)[1]<<
3938 "sigma="<<(*param)[2]<<
3939 "constant="<<(*param)[0]<<
3940 "\n";
3941 }
3942 }
3943
3944 if((chi2/(*param)[2]) > 0.1){
3945 if(kError){
3946 chi2 = FitGausMI(arraye,arraym,arrayme,nBins,xMin,xMax,param,kFALSE);
3947 }
3948 else return -4.0;
3949 }
3950
3951 if(fDebugLevel == 1){
3952 TString name("PRF");
3953 name += (Int_t)xMin;
3954 name += (Int_t)xMax;
3955 TCanvas *c1 = new TCanvas((const char *)name,(const char *)name,50,50,600,800);
3956 c1->cd();
3957 name += "histo";
3958 TH1F *histo = new TH1F((const char *)name,(const char *)name,nBins,xMin,xMax);
3959 for(Int_t k = 0; k < nBins; k++){
3960 histo->SetBinContent(k+1,arraym[k]);
3961 histo->SetBinError(k+1,arrayme[k]);
3962 }
3963 histo->Draw();
3964 name += "functionf";
3965 TF1 *f1= new TF1((const char*)name,"[0]*exp(-(x-[1])^2/(2*[2]*[2]))",xMin,xMax);
3966 f1->SetParameter(0, (*param)[0]);
3967 f1->SetParameter(1, (*param)[1]);
3968 f1->SetParameter(2, (*param)[2]);
3969 f1->Draw("same");
3970 }
3971
3972
3973 return chi2;
3974
3975}
3976//_____________________________________________________________________________
3977void AliTRDCalibraFit::FitPRF(TH1 *projPRF)
3978{
3979 //
3980 // Fit methode for the sigma of the pad response function
3981 //
3982
3983 fCurrentCoef[0] = 0.0;
3984 fCurrentCoefE = 0.0;
3985
3986 if (fDebugLevel != 1) {
3987 projPRF->Fit("gaus","0M","",-fRangeFitPRF,fRangeFitPRF);
3988 }
3989 else {
3990 TCanvas *cfit = new TCanvas("cfit","cfit",50,50,600,800);
3991 cfit->cd();
3992 projPRF->Fit("gaus","M+","",-fRangeFitPRF,fRangeFitPRF);
3993 projPRF->Draw();
3994 }
3995 fCurrentCoef[0] = projPRF->GetFunction("gaus")->GetParameter(2);
3996 fCurrentCoefE = projPRF->GetFunction("gaus")->GetParError(2);
3997 if(fCurrentCoef[0] <= 0.0) fCurrentCoef[0] = -fCurrentCoef[1];
3998 else {
3999 fNumberFitSuccess++;
4000 }
4001}
4002//_____________________________________________________________________________
4003void AliTRDCalibraFit::RmsPRF(TH1 *projPRF)
4004{
4005 //
4006 // Fit methode for the sigma of the pad response function
4007 //
4008 fCurrentCoef[0] = 0.0;
4009 fCurrentCoefE = 0.0;
4010 if (fDebugLevel == 1) {
4011 TCanvas *cfit = new TCanvas("cfit","cfit",50,50,600,800);
4012 cfit->cd();
4013 projPRF->Draw();
4014 }
4015 fCurrentCoef[0] = projPRF->GetRMS();
4016 if(fCurrentCoef[0] <= 0.0) fCurrentCoef[0] = -fCurrentCoef[1];
4017 else {
4018 fNumberFitSuccess++;
4019 }
4020}
4021//_____________________________________________________________________________
4022void AliTRDCalibraFit::FitTnpRange(Double_t *arraye, Double_t *arraym, Double_t *arrayme, Int_t nbg, Int_t nybins)
4023{
4024 //
4025 // Fit methode for the sigma of the pad response function with 2*nbg tan bins
4026 //
4027
4028 TLinearFitter linearfitter = TLinearFitter(3,"pol2");
4029
4030
4031 Int_t nbins = (Int_t)(nybins/(2*nbg));
4032 Float_t lowedge = -3.0*nbg;
4033 Float_t upedge = lowedge + 3.0;
4034 Int_t offset = 0;
4035 Int_t npoints = 0;
4036 Double_t xvalues = -0.2*nbg+0.1;
4037 Double_t y = 0.0;
4038 Int_t total = 2*nbg;
4039
4040
4041 for(Int_t k = 0; k < total; k++){
4042 if(FitPRFGausMI(arraye+offset, arraym+offset, arrayme+offset, nbins, lowedge, upedge)){
4043 npoints++;
4044 y = fCurrentCoef[0]*fCurrentCoef[0];
4045 linearfitter.AddPoint(&xvalues,y,2*fCurrentCoefE*fCurrentCoef[0]);
4046 }
4047
4048 if(fDebugLevel > 1){
4049
4050 if ( !fDebugStreamer ) {
4051 //debug stream
4052 TDirectory *backup = gDirectory;
4053 fDebugStreamer = new TTreeSRedirector("TRDDebugFit.root");
4054 if ( backup ) backup->cd(); //we don't want to be cd'd to the debug streamer
4055 }
4056
4057 Int_t detector = fCountDet;
4058 Int_t plane = GetPlane(fCountDet);
4059 Int_t nbtotal = total;
4060 Int_t group = k;
4061 Float_t low = lowedge;
4062 Float_t up = upedge;
4063 Float_t tnp = xvalues;
4064 Float_t wid = fCurrentCoef[0];
4065 Float_t widfE = fCurrentCoefE;
4066
4067 (* fDebugStreamer) << "FitTnpRange0"<<
4068 "detector="<<detector<<
4069 "plane="<<plane<<
4070 "nbtotal="<<nbtotal<<
4071 "group="<<group<<
4072 "low="<<low<<
4073 "up="<<up<<
4074 "offset="<<offset<<
4075 "tnp="<<tnp<<
4076 "wid="<<wid<<
4077 "widfE="<<widfE<<
4078 "\n";
4079 }
4080
4081 offset += nbins;
4082 lowedge += 3.0;
4083 upedge += 3.0;
4084 xvalues += 0.2;
4085
4086 }
4087
4088 fCurrentCoefE = 0.0;
4089 fCurrentCoef[0] = 0.0;
4090
4091 //printf("npoints\n",npoints);
4092
4093 if(npoints < 3){
4094 fCurrentCoef[0] = -TMath::Abs(fCurrentCoef[1]);
4095 }
4096 else{
4097
4098 TVectorD pars0;
4099 linearfitter.Eval();
4100 linearfitter.GetParameters(pars0);
4101 Double_t pointError0 = TMath::Sqrt(linearfitter.GetChisquare()/npoints);
4102 Double_t errorsx0 = linearfitter.GetParError(2)*pointError0;
4103 Double_t min0 = 0.0;
4104 Double_t ermin0 = 0.0;
4105 //Double_t prfe0 = 0.0;
4106 Double_t prf0 = 0.0;
4107 if((pars0[2] > 0.0) && (pars0[1] != 0.0)) {
4108 min0 = -pars0[1]/(2*pars0[2]);
4109 ermin0 = TMath::Abs(min0*(errorsx0/pars0[2]+linearfitter.GetParError(1)*pointError0/pars0[1]));
4110 prf0 = pars0[0]+pars0[1]*min0+pars0[2]*min0*min0;
4111 if(prf0 > 0.0) {
4112 /*
4113 prfe0 = linearfitter->GetParError(0)*pointError0
4114 +(linearfitter->GetParError(1)*pointError0/pars0[1]+ermin0/min0)*pars0[1]*min0
4115 +(linearfitter->GetParError(2)*pointError0/pars0[2]+2*ermin0/min0)*pars0[2]*min0*min0;
4116 prfe0 = prfe0/(2*TMath::Sqrt(prf0));
4117 fCurrentCoefE = (Float_t) prfe0;
4118 */
4119 fCurrentCoef[0] = (Float_t) TMath::Sqrt(TMath::Abs(prf0));
4120 }
4121 else{
4122 fCurrentCoef[0] = -TMath::Abs(fCurrentCoef[1]);
4123 }
4124 }
4125 else {
4126 fCurrentCoef[0] = -TMath::Abs(fCurrentCoef[1]);
4127 }
4128
4129 if(fDebugLevel > 1){
4130
4131 if ( !fDebugStreamer ) {
4132 //debug stream
4133 TDirectory *backup = gDirectory;
4134 fDebugStreamer = new TTreeSRedirector("TRDDebugFit.root");
4135 if ( backup ) backup->cd(); //we don't want to be cd'd to the debug streamer
4136 }
4137
4138 Int_t detector = fCountDet;
4139 Int_t plane = GetPlane(fCountDet);
4140 Int_t nbtotal = total;
4141 Double_t colsize[6] = {0.635,0.665,0.695,0.725,0.755,0.785};
4142 Double_t sigmax = TMath::Sqrt(TMath::Abs(pars0[2]))*10000*colsize[plane];
4143
4144 (* fDebugStreamer) << "FitTnpRange1"<<
4145 "detector="<<detector<<
4146 "plane="<<plane<<
4147 "nbtotal="<<nbtotal<<
4148 "par0="<<pars0[0]<<
4149 "par1="<<pars0[1]<<
4150 "par2="<<pars0[2]<<
4151 "npoints="<<npoints<<
4152 "sigmax="<<sigmax<<
4153 "tan="<<min0<<
4154 "sigmaprf="<<fCurrentCoef[0]<<
4155 "sigprf="<<fCurrentCoef[1]<<
4156 "\n";
4157 }
4158
4159 }
4160
4161}
4162//_____________________________________________________________________________
4163void AliTRDCalibraFit::FitMean(TH1 *projch, Double_t nentries, Double_t mean)
4164{
4165 //
4166 // Only mean methode for the gain factor
4167 //
4168
4169 fCurrentCoef[0] = mean;
4170 fCurrentCoefE = 0.0;
4171 if(nentries > 0) fCurrentCoefE = projch->GetRMS()/TMath::Sqrt(nentries);
4172 if (fDebugLevel == 1) {
4173 TCanvas *cpmean = new TCanvas("cpmean","cpmean",50,50,600,800);
4174 cpmean->cd();
4175 projch->Draw();
4176 }
4177 CalculChargeCoefMean(kTRUE);
4178 fNumberFitSuccess++;
4179}
4180//_____________________________________________________________________________
4181void AliTRDCalibraFit::FitMeanW(TH1 *projch, Double_t nentries)
4182{
4183 //
4184 // mean w methode for the gain factor
4185 //
4186
4187 //Number of bins
4188 Int_t nybins = projch->GetNbinsX();
4189
4190 //The weight function
4191 Double_t a = 0.00228515;
4192 Double_t b = -0.00231487;
4193 Double_t c = 0.00044298;
4194 Double_t d = -0.00379239;
4195 Double_t e = 0.00338349;
4196
4197// 0 |0.00228515
4198// 1 |-0.00231487
4199// 2 |0.00044298
4200// 3 |-0.00379239
4201// 4 |0.00338349
4202
4203
4204
4205 //A arbitrary error for the moment
4206 fCurrentCoefE = 0.0;
4207 fCurrentCoef[0] = 0.0;
4208
4209 //Calcul
4210 Double_t sumw = 0.0;
4211 Double_t sum = 0.0;
4212 Float_t sumAll = (Float_t) nentries;
4213 Int_t sumCurrent = 0;
4214 for(Int_t k = 0; k <nybins; k++){
4215 Double_t fraction = Float_t(sumCurrent)/Float_t(sumAll);
4216 if (fraction>0.95) break;
4217 Double_t weight = a + b*fraction + c*fraction*fraction + d *fraction*fraction*fraction+
4218 e*fraction*fraction*fraction*fraction;
4219 sumw += weight*projch->GetBinContent(k+1)*projch->GetBinCenter(k+1);
4220 sum += weight*projch->GetBinContent(k+1);
4221 sumCurrent += (Int_t) projch->GetBinContent(k+1);
4222 //printf("fraction %f, weight %f, bincontent %f\n",fraction,weight,projch->GetBinContent(k+1));
4223 }
4224 if(sum > 0.0) fCurrentCoef[0] = (sumw/sum);
4225
4226 if (fDebugLevel == 1) {
4227 TCanvas *cpmeanw = new TCanvas("cpmeanw","cpmeanw",50,50,600,800);
4228 cpmeanw->cd();
4229 projch->Draw();
4230 }
4231 fNumberFitSuccess++;
4232 CalculChargeCoefMean(kTRUE);
4233}
4234//_____________________________________________________________________________
4235void AliTRDCalibraFit::FitMeanWSm(TH1 *projch, Float_t sumAll)
4236{
4237 //
4238 // mean w methode for the gain factor
4239 //
4240
4241 //Number of bins
4242 Int_t nybins = projch->GetNbinsX();
4243
4244 //The weight function
4245 Double_t a = 0.00228515;
4246 Double_t b = -0.00231487;
4247 Double_t c = 0.00044298;
4248 Double_t d = -0.00379239;
4249 Double_t e = 0.00338349;
4250
4251// 0 |0.00228515
4252// 1 |-0.00231487
4253// 2 |0.00044298
4254// 3 |-0.00379239
4255// 4 |0.00338349
4256
4257
4258
4259 //A arbitrary error for the moment
4260 fCurrentCoefE = 0.0;
4261 fCurrentCoef[0] = 0.0;
4262
4263 //Calcul
4264 Double_t sumw = 0.0;
4265 Double_t sum = 0.0;
4266 Int_t sumCurrent = 0;
4267 for(Int_t k = 0; k <nybins; k++){
4268 Double_t fraction = Float_t(sumCurrent)/Float_t(sumAll);
4269 if (fraction>0.95) break;
4270 Double_t weight = a + b*fraction + c*fraction*fraction + d *fraction*fraction*fraction+
4271 e*fraction*fraction*fraction*fraction;
4272 sumw += weight*projch->GetBinContent(k+1)*projch->GetBinCenter(k+1);
4273 sum += weight*projch->GetBinContent(k+1);
4274 sumCurrent += (Int_t) projch->GetBinContent(k+1);
4275 //printf("fraction %f, weight %f, bincontent %f\n",fraction,weight,projch->GetBinContent(k+1));
4276 }
4277 if(sum > 0.0) fCurrentCoef[0] = (sumw/sum);
4278
4279 if (fDebugLevel == 1) {
4280 TCanvas *cpmeanw = new TCanvas("cpmeanw","cpmeanw",50,50,600,800);
4281 cpmeanw->cd();
4282 projch->Draw();
4283 }
4284 fNumberFitSuccess++;
4285}
4286//_____________________________________________________________________________
4287void AliTRDCalibraFit::FitCH(TH1 *projch, Double_t mean)
4288{
4289 //
4290 // Fit methode for the gain factor
4291 //
4292
4293 fCurrentCoef[0] = 0.0;
4294 fCurrentCoefE = 0.0;
4295 Double_t chisqrl = 0.0;
4296 Double_t chisqrg = 0.0;
4297 Double_t chisqr = 0.0;
4298 TF1 *fLandauGaus = new TF1("fLandauGaus",FuncLandauGaus,0,300,5);
4299
4300 projch->Fit("landau","0",""
4301 ,(Double_t) mean/fBeginFitCharge
4302 ,projch->GetBinCenter(projch->GetNbinsX()));
4303 Double_t l3P0 = projch->GetFunction("landau")->GetParameter(0);
4304 Double_t l3P1 = projch->GetFunction("landau")->GetParameter(1);
4305 Double_t l3P2 = projch->GetFunction("landau")->GetParameter(2);
4306 chisqrl = projch->GetFunction("landau")->GetChisquare();
4307
4308 projch->Fit("gaus","0",""
4309 ,(Double_t) mean/fBeginFitCharge
4310 ,projch->GetBinCenter(projch->GetNbinsX()));
4311 Double_t g3P0 = projch->GetFunction("gaus")->GetParameter(0);
4312 Double_t g3P2 = projch->GetFunction("gaus")->GetParameter(2);
4313 chisqrg = projch->GetFunction("gaus")->GetChisquare();
4314
4315 fLandauGaus->SetParameters(l3P0,l3P1,l3P2,g3P0,g3P2);
4316 if (fDebugLevel != 1) {
4317 projch->Fit("fLandauGaus","0",""
4318 ,(Double_t) mean/fBeginFitCharge
4319 ,projch->GetBinCenter(projch->GetNbinsX()));
4320 chisqr = projch->GetFunction("fLandauGaus")->GetChisquare();
4321 }
4322 else {
4323 TCanvas *cp = new TCanvas("cp","cp",50,50,600,800);
4324 cp->cd();
4325 projch->Fit("fLandauGaus","+",""
4326 ,(Double_t) mean/fBeginFitCharge
4327 ,projch->GetBinCenter(projch->GetNbinsX()));
4328 chisqr = projch->GetFunction("fLandauGaus")->GetChisquare();
4329 projch->Draw();
4330 fLandauGaus->Draw("same");
4331 }
4332
4333 if ((projch->GetFunction("fLandauGaus")->GetParameter(1) > 0) && (projch->GetFunction("fLandauGaus")->GetParError(1) < (0.05*projch->GetFunction("fLandauGaus")->GetParameter(1))) && (chisqr < chisqrl) && (chisqr < chisqrg)) {
4334 //if ((projch->GetFunction("fLandauGaus")->GetParameter(1) > 0) && (chisqr < chisqrl) && (chisqr < chisqrg)) {
4335 fNumberFitSuccess++;
4336 CalculChargeCoefMean(kTRUE);
4337 fCurrentCoef[0] = projch->GetFunction("fLandauGaus")->GetParameter(1);
4338 fCurrentCoefE = projch->GetFunction("fLandauGaus")->GetParError(1);
4339 }
4340 else {
4341 CalculChargeCoefMean(kFALSE);
4342 fCurrentCoef[0] = -TMath::Abs(fCurrentCoef[1]);
4343 }
4344
4345 if (fDebugLevel != 1) {
4346 delete fLandauGaus;
4347 }
4348
4349}
4350//_____________________________________________________________________________
4351void AliTRDCalibraFit::FitBisCH(TH1* projch, Double_t mean)
4352{
4353 //
4354 // Fit methode for the gain factor more time consuming
4355 //
4356
4357
4358 //Some parameters to initialise
4359 Double_t widthLandau, widthGaus, mPV, integral;
4360 Double_t chisquarel = 0.0;
4361 Double_t chisquareg = 0.0;
4362 projch->Fit("landau","0M+",""
4363 ,(Double_t) mean/6
4364 ,projch->GetBinCenter(projch->GetNbinsX()));
4365 widthLandau = projch->GetFunction("landau")->GetParameter(2);
4366 chisquarel = projch->GetFunction("landau")->GetChisquare();
4367 projch->Fit("gaus","0M+",""
4368 ,(Double_t) mean/6
4369 ,projch->GetBinCenter(projch->GetNbinsX()));
4370 widthGaus = projch->GetFunction("gaus")->GetParameter(2);
4371 chisquareg = projch->GetFunction("gaus")->GetChisquare();
4372
4373 mPV = (projch->GetFunction("landau")->GetParameter(1))/2;
4374 integral = (projch->GetFunction("gaus")->Integral(0.3*mean,3*mean)+projch->GetFunction("landau")->Integral(0.3*mean,3*mean))/2;
4375
4376 // Setting fit range and start values
4377 Double_t fr[2];
4378 //Double_t sv[4] = { l3P2, fChargeCoef[1], projch->Integral("width"), fG3P2 };
4379 //Double_t sv[4] = { fL3P2, fChargeCoef[1], fL3P0, fG3P2 };
4380 Double_t sv[4] = { widthLandau, mPV, integral, widthGaus};
4381 Double_t pllo[4] = { 0.001, 0.001, projch->Integral()/3, 0.001};
4382 Double_t plhi[4] = { 300.0, 300.0, 30*projch->Integral(), 300.0};
4383 Double_t fp[4] = { 1.0, 1.0, 1.0, 1.0 };
4384 Double_t fpe[4] = { 1.0, 1.0, 1.0, 1.0 };
4385 fr[0] = 0.3 * mean;
4386 fr[1] = 3.0 * mean;
4387 fCurrentCoef[0] = 0.0;
4388 fCurrentCoefE = 0.0;
4389
4390 Double_t chisqr;
4391 Int_t ndf;
4392 TF1 *fitsnr = LanGauFit(projch,&fr[0],&sv[0]
4393 ,&pllo[0],&plhi[0]
4394 ,&fp[0],&fpe[0]
4395 ,&chisqr,&ndf);
4396
4397 Double_t projchPeak;
4398 Double_t projchFWHM;
4399 LanGauPro(fp,projchPeak,projchFWHM);
4400
4401 if ((fp[1] > 0) && ((fpe[1] < (0.05*fp[1])) && (chisqr < chisquarel) && (chisqr < chisquareg))) {
4402 //if ((fp[1] > 0) && ((chisqr < chisquarel) && (chisqr < chisquareg))) {
4403 fNumberFitSuccess++;
4404 CalculChargeCoefMean(kTRUE);
4405 fCurrentCoef[0] = fp[1];
4406 fCurrentCoefE = fpe[1];
4407 //chargeCoefE2 = chisqr;
4408 }
4409 else {
4410 CalculChargeCoefMean(kFALSE);
4411 fCurrentCoef[0] = -TMath::Abs(fCurrentCoef[1]);
4412 }
4413 if (fDebugLevel == 1) {
4414 AliInfo(Form("fChargeCoef[0]: %f",(Float_t) fCurrentCoef[0]));
4415 TCanvas *cpy = new TCanvas("cpy","cpy",50,50,600,800);
4416 cpy->cd();
4417 projch->Draw();
4418 fitsnr->Draw("same");
4419 }
4420 else {
4421 delete fitsnr;
4422 }
4423}
4424//_____________________________________________________________________________
4425Double_t *AliTRDCalibraFit::CalculPolynomeLagrange2(Double_t *x, Double_t *y) const
4426{
4427 //
4428 // Calcul the coefficients of the polynome passant par ces trois points de degre 2
4429 //
4430 Double_t *c = new Double_t[5];
4431 Double_t x0 = y[0]/((x[0]-x[1])*(x[0]-x[2]));
4432 Double_t x1 = y[1]/((x[1]-x[0])*(x[1]-x[2]));
4433 Double_t x2 = y[2]/((x[2]-x[0])*(x[2]-x[1]));
4434
4435 c[4] = 0.0;
4436 c[3] = 0.0;
4437 c[2] = x0+x1+x2;
4438 c[1] = -(x0*(x[1]+x[2])+x1*(x[0]+x[2])+x2*(x[0]+x[1]));
4439 c[0] = x0*x[1]*x[2]+x1*x[0]*x[2]+x2*x[0]*x[1];
4440
4441 return c;
4442
4443
4444}
4445
4446//_____________________________________________________________________________
4447Double_t *AliTRDCalibraFit::CalculPolynomeLagrange3(Double_t *x, Double_t *y) const
4448{
4449 //
4450 // Calcul the coefficients of the polynome passant par ces quatre points de degre 3
4451 //
4452 Double_t *c = new Double_t[5];
4453 Double_t x0 = y[0]/((x[0]-x[1])*(x[0]-x[2])*(x[0]-x[3]));
4454 Double_t x1 = y[1]/((x[1]-x[0])*(x[1]-x[2])*(x[1]-x[3]));
4455 Double_t x2 = y[2]/((x[2]-x[0])*(x[2]-x[1])*(x[2]-x[3]));
4456 Double_t x3 = y[3]/((x[3]-x[0])*(x[3]-x[1])*(x[3]-x[2]));
4457
4458 c[4] = 0.0;
4459 c[3] = x0+x1+x2+x3;
4460 c[2] = -(x0*(x[1]+x[2]+x[3])
4461 +x1*(x[0]+x[2]+x[3])
4462 +x2*(x[0]+x[1]+x[3])
4463 +x3*(x[0]+x[1]+x[2]));
4464 c[1] = (x0*(x[1]*x[2]+x[1]*x[3]+x[2]*x[3])
4465 +x1*(x[0]*x[2]+x[0]*x[3]+x[2]*x[3])
4466 +x2*(x[0]*x[1]+x[0]*x[3]+x[1]*x[3])
4467 +x3*(x[0]*x[1]+x[0]*x[2]+x[1]*x[2]));
4468
4469 c[0] = -(x0*x[1]*x[2]*x[3]
4470 +x1*x[0]*x[2]*x[3]
4471 +x2*x[0]*x[1]*x[3]
4472 +x3*x[0]*x[1]*x[2]);
4473
4474
4475 return c;
4476
4477
4478}
4479
4480//_____________________________________________________________________________
4481Double_t *AliTRDCalibraFit::CalculPolynomeLagrange4(Double_t *x, Double_t *y) const
4482{
4483 //
4484 // Calcul the coefficients of the polynome passant par ces cinqs points de degre 4
4485 //
4486 Double_t *c = new Double_t[5];
4487 Double_t x0 = y[0]/((x[0]-x[1])*(x[0]-x[2])*(x[0]-x[3])*(x[0]-x[4]));
4488 Double_t x1 = y[1]/((x[1]-x[0])*(x[1]-x[2])*(x[1]-x[3])*(x[1]-x[4]));
4489 Double_t x2 = y[2]/((x[2]-x[0])*(x[2]-x[1])*(x[2]-x[3])*(x[2]-x[4]));
4490 Double_t x3 = y[3]/((x[3]-x[0])*(x[3]-x[1])*(x[3]-x[2])*(x[3]-x[4]));
4491 Double_t x4 = y[4]/((x[4]-x[0])*(x[4]-x[1])*(x[4]-x[2])*(x[4]-x[3]));
4492
4493
4494 c[4] = x0+x1+x2+x3+x4;
4495 c[3] = -(x0*(x[1]+x[2]+x[3]+x[4])
4496 +x1*(x[0]+x[2]+x[3]+x[4])
4497 +x2*(x[0]+x[1]+x[3]+x[4])
4498 +x3*(x[0]+x[1]+x[2]+x[4])
4499 +x4*(x[0]+x[1]+x[2]+x[3]));
4500 c[2] = (x0*(x[1]*x[2]+x[1]*x[3]+x[1]*x[4]+x[2]*x[3]+x[2]*x[4]+x[3]*x[4])
4501 +x1*(x[0]*x[2]+x[0]*x[3]+x[0]*x[4]+x[2]*x[3]+x[2]*x[4]+x[3]*x[4])
4502 +x2*(x[0]*x[1]+x[0]*x[3]+x[0]*x[4]+x[1]*x[3]+x[1]*x[4]+x[3]*x[4])
4503 +x3*(x[0]*x[1]+x[0]*x[2]+x[0]*x[4]+x[1]*x[2]+x[1]*x[4]+x[2]*x[4])
4504 +x4*(x[0]*x[1]+x[0]*x[2]+x[0]*x[3]+x[1]*x[2]+x[1]*x[3]+x[2]*x[3]));
4505
4506 c[1] = -(x0*(x[1]*x[2]*x[3]+x[1]*x[2]*x[4]+x[1]*x[3]*x[4]+x[2]*x[3]*x[4])
4507 +x1*(x[0]*x[2]*x[3]+x[0]*x[2]*x[4]+x[0]*x[3]*x[4]+x[2]*x[3]*x[4])
4508 +x2*(x[0]*x[1]*x[3]+x[0]*x[1]*x[4]+x[0]*x[3]*x[4]+x[1]*x[3]*x[4])
4509 +x3*(x[0]*x[1]*x[2]+x[0]*x[1]*x[4]+x[0]*x[2]*x[4]+x[1]*x[2]*x[4])
4510 +x4*(x[0]*x[1]*x[2]+x[0]*x[1]*x[3]+x[0]*x[2]*x[3]+x[1]*x[2]*x[3]));
4511
4512 c[0] = (x0*x[1]*x[2]*x[3]*x[4]
4513 +x1*x[0]*x[2]*x[3]*x[4]
4514 +x2*x[0]*x[1]*x[3]*x[4]
4515 +x3*x[0]*x[1]*x[2]*x[4]
4516 +x4*x[0]*x[1]*x[2]*x[3]);
4517
4518 return c;
4519
4520
4521}
4522//_____________________________________________________________________________
4523void AliTRDCalibraFit::NormierungCharge()
4524{
4525 //
4526 // Normalisation of the gain factor resulting for the fits
4527 //
4528
4529 // Calcul of the mean of choosen method by fFitChargeNDB
4530 Double_t sum = 0.0;
4531 //printf("total number of entries %d\n",fVectorFitCH->GetEntriesFast());
4532 for (Int_t k = 0; k < (Int_t) fVectorFit.GetEntriesFast(); k++) {
4533 Int_t total = 0;
4534 Int_t detector = ((AliTRDFitInfo *) fVectorFit.At(k))->GetDetector();
4535 Float_t *coef = ((AliTRDFitInfo *) fVectorFit.At(k))->GetCoef();
4536 //printf("detector %d coef[0] %f\n",detector,coef[0]);
4537 if (GetChamber(detector) == 2) {
4538 total = 1728;
4539 }
4540 if (GetChamber(detector) != 2) {
4541 total = 2304;
4542 }
4543 for (Int_t j = 0; j < total; j++) {
4544 if (coef[j] >= 0) {
4545 sum += coef[j];
4546 }
4547 }
4548 }
4549
4550 if (sum > 0) {
4551 fScaleFitFactor = fScaleFitFactor / sum;
4552 }
4553 else {
4554 fScaleFitFactor = 1.0;
4555 }
4556
4557 //methode de boeuf mais bon...
4558 Double_t scalefactor = fScaleFitFactor;
4559
4560 if(fDebugLevel > 1){
4561
4562 if ( !fDebugStreamer ) {
4563 //debug stream
4564 TDirectory *backup = gDirectory;
4565 fDebugStreamer = new TTreeSRedirector("TRDDebugFit.root");
4566 if ( backup ) backup->cd(); //we don't want to be cd'd to the debug streamer
4567 }
4568 (* fDebugStreamer) << "NormierungCharge"<<
4569 "scalefactor="<<scalefactor<<
4570 "\n";
4571 }
4572}
4573//_____________________________________________________________________________
4574TH1I *AliTRDCalibraFit::ReBin(TH1I *hist) const
4575{
4576 //
4577 // Rebin of the 1D histo for the gain calibration if needed.
4578 // you have to choose fRebin, divider of fNumberBinCharge
4579 //
4580
4581 TAxis *xhist = hist->GetXaxis();
4582 TH1I *rehist = new TH1I("projrebin","",(Int_t) xhist->GetNbins()/fRebin
4583 ,xhist->GetBinLowEdge(1)
4584 ,xhist->GetBinUpEdge(xhist->GetNbins()));
4585
4586 AliInfo(Form("fRebin: %d",fRebin));
4587 Int_t i = 1;
4588 for (Int_t k = 1; k <= (Int_t) xhist->GetNbins()/fRebin; k++) {
4589 Double_t sum = 0.0;
4590 for (Int_t ji = i; ji < i+fRebin; ji++) {
4591 sum += hist->GetBinContent(ji);
4592 }
4593 sum = sum / fRebin;
4594 rehist->SetBinContent(k,sum);
4595 i += fRebin;
4596 }
4597
4598 return rehist;
4599
4600}
4601
4602//_____________________________________________________________________________
4603TH1F *AliTRDCalibraFit::ReBin(TH1F *hist) const
4604{
4605 //
4606 // Rebin of the 1D histo for the gain calibration if needed
4607 // you have to choose fRebin divider of fNumberBinCharge
4608 //
4609
4610 TAxis *xhist = hist->GetXaxis();
4611 TH1F *rehist = new TH1F("projrebin","",(Int_t) xhist->GetNbins()/fRebin
4612 ,xhist->GetBinLowEdge(1)
4613 ,xhist->GetBinUpEdge(xhist->GetNbins()));
4614
4615 AliInfo(Form("fRebin: %d",fRebin));
4616 Int_t i = 1;
4617 for (Int_t k = 1; k <= (Int_t) xhist->GetNbins()/fRebin; k++) {
4618 Double_t sum = 0.0;
4619 for (Int_t ji = i; ji < i+fRebin; ji++) {
4620 sum += hist->GetBinContent(ji);
4621 }
4622 sum = sum/fRebin;
4623 rehist->SetBinContent(k,sum);
4624 i += fRebin;
4625 }
4626
4627 return rehist;
4628
4629}
4630
4631//_____________________________________________________________________________
4632TH1F *AliTRDCalibraFit::CorrectTheError(TGraphErrors *hist)
4633{
4634 //
4635 // In the case of the vectors method the trees contains TGraphErrors for PH and PRF
4636 // to be able to add them after
4637 // We convert it to a TH1F to be able to applied the same fit function method
4638 // After having called this function you can not add the statistics anymore
4639 //
4640
4641 TH1F *rehist = 0x0;
4642
4643 Int_t nbins = hist->GetN();
4644 Double_t *x = hist->GetX();
4645 Double_t *entries = hist->GetEX();
4646 Double_t *mean = hist->GetY();
4647 Double_t *square = hist->GetEY();
4648 fEntriesCurrent = 0;
4649
4650 if (nbins < 2) {
4651 return rehist;
4652 }
4653
4654 Double_t step = x[1] - x[0];
4655 Double_t minvalue = x[0] - step/2;
4656 Double_t maxvalue = x[(nbins-1)] + step/2;
4657
4658 rehist = new TH1F("projcorrecterror","",nbins,minvalue,maxvalue);
4659
4660 for (Int_t k = 0; k < nbins; k++) {
4661 rehist->SetBinContent(k+1,mean[k]);
4662 if (entries[k] > 0.0) {
4663 fEntriesCurrent += (Int_t) entries[k];
4664 Double_t d = TMath::Abs(square[k] - (mean[k]*mean[k]));
4665 rehist->SetBinError(k+1,TMath::Sqrt(d/entries[k]));
4666 }
4667 else {
4668 rehist->SetBinError(k+1,0.0);
4669 }
4670 }
4671
4672 if(fEntriesCurrent > 0) fNumberEnt++;
4673
4674 return rehist;
4675
4676}
4677//
4678//____________Some basic geometry function_____________________________________
4679//
4680
4681//_____________________________________________________________________________
4682Int_t AliTRDCalibraFit::GetPlane(Int_t d) const
4683{
4684 //
4685 // Reconstruct the plane number from the detector number
4686 //
4687
4688 return ((Int_t) (d % 6));
4689
4690}
4691
4692//_____________________________________________________________________________
4693Int_t AliTRDCalibraFit::GetChamber(Int_t d) const
4694{
4695 //
4696 // Reconstruct the chamber number from the detector number
4697 //
4698 Int_t fgkNplan = 6;
4699
4700 return ((Int_t) (d % 30) / fgkNplan);
4701
4702}
4703
4704//_____________________________________________________________________________
4705Int_t AliTRDCalibraFit::GetSector(Int_t d) const
4706{
4707 //
4708 // Reconstruct the sector number from the detector number
4709 //
4710 Int_t fg = 30;
4711
4712 return ((Int_t) (d / fg));
4713
4714}
4715
4716//
4717//____________Fill and Init tree Gain, PRF, Vdrift and T0______________________
4718//
4719//_______________________________________________________________________________
4720void AliTRDCalibraFit::ResetVectorFit()
4721{
4722 //
4723 // Reset the VectorFits
4724 //
4725
4726 fVectorFit.SetOwner();
4727 fVectorFit.Clear();
4728 fVectorFit2.SetOwner();
4729 fVectorFit2.Clear();
4730
4731}
4732//
4733//____________Private Functions________________________________________________
4734//
4735
4736//_____________________________________________________________________________
4737Double_t AliTRDCalibraFit::PH(Double_t *x, Double_t *par)
4738{
4739 //
4740 // Function for the fit
4741 //
4742
4743 //TF1 *fAsymmGauss = new TF1("fAsymmGauss",AsymmGauss,0,4,6);
4744
4745 //PARAMETERS FOR FIT PH
4746 // PASAv.4
4747 //fAsymmGauss->SetParameter(0,0.113755);
4748 //fAsymmGauss->SetParameter(1,0.350706);
4749 //fAsymmGauss->SetParameter(2,0.0604244);
4750 //fAsymmGauss->SetParameter(3,7.65596);
4751 //fAsymmGauss->SetParameter(4,1.00124);
4752 //fAsymmGauss->SetParameter(5,0.870597); // No tail cancelation
4753
4754 Double_t xx = x[0];
4755
4756 if (xx < par[1]) {
4757 return par[5];
4758 }
4759
4760 Double_t dx = 0.005;
4761 Double_t xs = par[1];
4762 Double_t ss = 0.0;
4763 Double_t paras[2] = { 0.0, 0.0 };
4764
4765 while (xs < xx) {
4766 if ((xs >= par[1]) &&
4767 (xs < (par[1]+par[2]))) {
4768 //fAsymmGauss->SetParameter(0,par[0]);
4769 //fAsymmGauss->SetParameter(1,xs);
4770 //ss += fAsymmGauss->Eval(xx);
4771 paras[0] = par[0];
4772 paras[1] = xs;
4773 ss += AsymmGauss(&xx,paras);
4774 }
4775 if ((xs >= (par[1]+par[2])) &&
4776 (xs < (par[1]+par[2]+par[3]))) {
4777 //fAsymmGauss->SetParameter(0,par[0]*par[4]);
4778 //fAsymmGauss->SetParameter(1,xs);
4779 //ss += fAsymmGauss->Eval(xx);
4780 paras[0] = par[0]*par[4];
4781 paras[1] = xs;
4782 ss += AsymmGauss(&xx,paras);
4783 }
4784 xs += dx;
4785 }
4786
4787 return ss + par[5];
4788
4789}
4790
4791//_____________________________________________________________________________
4792Double_t AliTRDCalibraFit::AsymmGauss(Double_t *x, Double_t *par)
4793{
4794 //
4795 // Function for the fit
4796 //
4797
4798 //par[0] = normalization
4799 //par[1] = mean
4800 //par[2] = sigma
4801 //norm0 = 1
4802 //par[3] = lambda0
4803 //par[4] = norm1
4804 //par[5] = lambda1
4805
4806 Double_t par1save = par[1];
4807 //Double_t par2save = par[2];
4808 Double_t par2save = 0.0604244;
4809 //Double_t par3save = par[3];
4810 Double_t par3save = 7.65596;
4811 //Double_t par5save = par[5];
4812 Double_t par5save = 0.870597;
4813 Double_t dx = x[0] - par1save;
4814
4815 Double_t sigma2 = par2save*par2save;
4816 Double_t sqrt2 = TMath::Sqrt(2.0);
4817 Double_t exp1 = par3save * TMath::Exp(-par3save * (dx - 0.5 * par3save * sigma2))
4818 * (1.0 - TMath::Erf((par3save * sigma2 - dx) / (sqrt2 * par2save)));
4819 Double_t exp2 = par5save * TMath::Exp(-par5save * (dx - 0.5 * par5save * sigma2))
4820 * (1.0 - TMath::Erf((par5save * sigma2 - dx) / (sqrt2 * par2save)));
4821
4822 //return par[0]*(exp1+par[4]*exp2);
4823 return par[0] * (exp1 + 1.00124 * exp2);
4824
4825}
4826
4827//_____________________________________________________________________________
4828Double_t AliTRDCalibraFit::FuncLandauGaus(Double_t *x, Double_t *par)
4829{
4830 //
4831 // Sum Landau + Gaus with identical mean
4832 //
4833
4834 Double_t valLandau = par[0] * TMath::Landau(x[0],par[1],par[2]);
4835 //Double_t valGaus = par[3] * TMath::Gaus(x[0],par[4],par[5]);
4836 Double_t valGaus = par[3] * TMath::Gaus(x[0],par[1],par[4]);
4837 Double_t val = valLandau + valGaus;
4838
4839 return val;
4840
4841}
4842
4843//_____________________________________________________________________________
4844Double_t AliTRDCalibraFit::LanGauFun(Double_t *x, Double_t *par)
4845{
4846 //
4847 // Function for the fit
4848 //
4849 // Fit parameters:
4850 // par[0]=Width (scale) parameter of Landau density
4851 // par[1]=Most Probable (MP, location) parameter of Landau density
4852 // par[2]=Total area (integral -inf to inf, normalization constant)
4853 // par[3]=Width (sigma) of convoluted Gaussian function
4854 //
4855 // In the Landau distribution (represented by the CERNLIB approximation),
4856 // the maximum is located at x=-0.22278298 with the location parameter=0.
4857 // This shift is corrected within this function, so that the actual
4858 // maximum is identical to the MP parameter.
4859 //
4860
4861 // Numeric constants
4862 Double_t invsq2pi = 0.3989422804014; // (2 pi)^(-1/2)
4863 Double_t mpshift = -0.22278298; // Landau maximum location
4864
4865 // Control constants
4866 Double_t np = 100.0; // Number of convolution steps
4867 Double_t sc = 5.0; // Convolution extends to +-sc Gaussian sigmas
4868
4869 // Variables
4870 Double_t xx;
4871 Double_t mpc;
4872 Double_t fland;
4873 Double_t sum = 0.0;
4874 Double_t xlow;
4875 Double_t xupp;
4876 Double_t step;
4877 Double_t i;
4878
4879 // MP shift correction
4880 mpc = par[1] - mpshift * par[0];
4881
4882 // Range of convolution integral
4883 xlow = x[0] - sc * par[3];
4884 xupp = x[0] + sc * par[3];
4885
4886 step = (xupp - xlow) / np;
4887
4888 // Convolution integral of Landau and Gaussian by sum
4889 for (i = 1.0; i <= np/2; i++) {
4890
4891 xx = xlow + (i-.5) * step;
4892 fland = TMath::Landau(xx,mpc,par[0]) / par[0];
4893 sum += fland * TMath::Gaus(x[0],xx,par[3]);
4894
4895 xx = xupp - (i-.5) * step;
4896 fland = TMath::Landau(xx,mpc,par[0]) / par[0];
4897 sum += fland * TMath::Gaus(x[0],xx,par[3]);
4898
4899 }
4900
4901 return (par[2] * step * sum * invsq2pi / par[3]);
4902
4903}
4904//_____________________________________________________________________________
4905TF1 *AliTRDCalibraFit::LanGauFit(TH1 *his, Double_t *fitrange, Double_t *startvalues
4906 , Double_t *parlimitslo, Double_t *parlimitshi
4907 , Double_t *fitparams, Double_t *fiterrors
4908 , Double_t *chiSqr, Int_t *ndf) const
4909{
4910 //
4911 // Function for the fit
4912 //
4913
4914 Int_t i;
4915 Char_t funname[100];
4916
4917 TF1 *ffitold = (TF1 *) gROOT->GetListOfFunctions()->FindObject(funname);
4918 if (ffitold) {
4919 delete ffitold;
4920 }
4921
4922 TF1 *ffit = new TF1(funname,LanGauFun,fitrange[0],fitrange[1],4);
4923 ffit->SetParameters(startvalues);
4924 ffit->SetParNames("Width","MP","Area","GSigma");
4925
4926 for (i = 0; i < 4; i++) {
4927 ffit->SetParLimits(i,parlimitslo[i],parlimitshi[i]);
4928 }
4929
4930 his->Fit(funname,"RB0"); // Fit within specified range, use ParLimits, do not plot
4931
4932 ffit->GetParameters(fitparams); // Obtain fit parameters
4933 for (i = 0; i < 4; i++) {
4934 fiterrors[i] = ffit->GetParError(i); // Obtain fit parameter errors
4935 }
4936 chiSqr[0] = ffit->GetChisquare(); // Obtain chi^2
4937 ndf[0] = ffit->GetNDF(); // Obtain ndf
4938
4939 return (ffit); // Return fit function
4940
4941}
4942
4943//_____________________________________________________________________________
4944Int_t AliTRDCalibraFit::LanGauPro(Double_t *params, Double_t &maxx, Double_t &fwhm)
4945{
4946 //
4947 // Function for the fit
4948 //
4949
4950 Double_t p;
4951 Double_t x;
4952 Double_t fy;
4953 Double_t fxr;
4954 Double_t fxl;
4955 Double_t step;
4956 Double_t l;
4957 Double_t lold;
4958
4959 Int_t i = 0;
4960 Int_t maxcalls = 10000;
4961
4962 // Search for maximum
4963 p = params[1] - 0.1 * params[0];
4964 step = 0.05 * params[0];
4965 lold = -2.0;
4966 l = -1.0;
4967
4968 while ((l != lold) && (i < maxcalls)) {
4969 i++;
4970 lold = l;
4971 x = p + step;
4972 l = LanGauFun(&x,params);
4973 if (l < lold) {
4974 step = -step / 10.0;
4975 }
4976 p += step;
4977 }
4978
4979 if (i == maxcalls) {
4980 return (-1);
4981 }
4982 maxx = x;
4983 fy = l / 2.0;
4984
4985 // Search for right x location of fy
4986 p = maxx + params[0];
4987 step = params[0];
4988 lold = -2.0;
4989 l = -1e300;
4990 i = 0;
4991
4992 while ( (l != lold) && (i < maxcalls) ) {
4993 i++;
4994
4995 lold = l;
4996 x = p + step;
4997 l = TMath::Abs(LanGauFun(&x,params) - fy);
4998
4999 if (l > lold)
5000 step = -step/10;
5001
5002 p += step;
5003 }
5004
5005 if (i == maxcalls)
5006 return (-2);
5007
5008 fxr = x;
5009
5010
5011 // Search for left x location of fy
5012
5013 p = maxx - 0.5 * params[0];
5014 step = -params[0];
5015 lold = -2.0;
5016 l = -1.0e300;
5017 i = 0;
5018
5019 while ((l != lold) && (i < maxcalls)) {
5020 i++;
5021 lold = l;
5022 x = p + step;
5023 l = TMath::Abs(LanGauFun(&x,params) - fy);
5024 if (l > lold) {
5025 step = -step / 10.0;
5026 }
5027 p += step;
5028 }
5029
5030 if (i == maxcalls) {
5031 return (-3);
5032 }
5033
5034 fxl = x;
5035 fwhm = fxr - fxl;
5036
5037 return (0);
5038}
5039//_____________________________________________________________________________
5040Double_t AliTRDCalibraFit::GausConstant(Double_t *x, Double_t *par)
5041{
5042 //
5043 // Gaus with identical mean
5044 //
5045
5046 Double_t gauss = par[0] * TMath::Gaus(x[0],0.0,par[1])+par[2];
5047
5048 return gauss;
5049
5050}