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