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