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