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