1 /**************************************************************************
2 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
4 * Author: The ALICE Off-line Project. *
5 * Contributors are mentioned in the code where appropriate. *
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 **************************************************************************/
18 /////////////////////////////////////////////////////////////////////////////////
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
45 // R. Bailhache (R.Bailhache@gsi.de)
47 //////////////////////////////////////////////////////////////////////////////////////
52 #include <TProfile2D.h>
54 #include <TGraphErrors.h>
55 #include <TObjArray.h>
61 #include <TDirectory.h>
63 #include <TTreeStream.h>
64 #include <TLinearFitter.h>
69 #include "AliMathBase.h"
71 #include "AliTRDCalibraFit.h"
72 #include "AliTRDCalibraMode.h"
73 #include "AliTRDCalibraVector.h"
74 #include "AliTRDCalibraVdriftLinearFit.h"
75 #include "AliTRDcalibDB.h"
76 #include "AliTRDgeometry.h"
77 #include "AliTRDpadPlane.h"
78 #include "AliTRDgeometry.h"
79 #include "./Cal/AliTRDCalROC.h"
80 #include "./Cal/AliTRDCalPad.h"
81 #include "./Cal/AliTRDCalDet.h"
84 ClassImp(AliTRDCalibraFit)
86 AliTRDCalibraFit* AliTRDCalibraFit::fgInstance = 0;
87 Bool_t AliTRDCalibraFit::fgTerminated = kFALSE;
89 //_____________singleton implementation_________________________________________________
90 AliTRDCalibraFit *AliTRDCalibraFit::Instance()
93 // Singleton implementation
96 if (fgTerminated != kFALSE) {
100 if (fgInstance == 0) {
101 fgInstance = new AliTRDCalibraFit();
108 //______________________________________________________________________________________
109 void AliTRDCalibraFit::Terminate()
112 // Singleton implementation
113 // Deletes the instance of this class
116 fgTerminated = kTRUE;
118 if (fgInstance != 0) {
125 //______________________________________________________________________________________
126 AliTRDCalibraFit::AliTRDCalibraFit()
129 ,fNumberOfBinsExpected(0)
131 ,fBeginFitCharge(3.5)
133 ,fTakeTheMaxPH(kTRUE)
141 ,fNumberFitSuccess(0)
148 ,fCalibraMode(new AliTRDCalibraMode())
153 ,fScaleFitFactor(0.0)
161 ,fCurrentCoefDetector(0x0)
162 ,fCurrentCoefDetector2(0x0)
167 // Default constructor
170 fGeo = new AliTRDgeometry();
172 // Current variables initialised
173 for (Int_t k = 0; k < 2; k++) {
174 fCurrentCoef[k] = 0.0;
175 fCurrentCoef2[k] = 0.0;
177 for (Int_t i = 0; i < 3; i++) {
183 //______________________________________________________________________________________
184 AliTRDCalibraFit::AliTRDCalibraFit(const AliTRDCalibraFit &c)
187 ,fNumberOfBinsExpected(c.fNumberOfBinsExpected)
189 ,fBeginFitCharge(c.fBeginFitCharge)
190 ,fFitPHPeriode(c.fFitPHPeriode)
191 ,fTakeTheMaxPH(c.fTakeTheMaxPH)
192 ,fT0Shift0(c.fT0Shift0)
193 ,fT0Shift1(c.fT0Shift1)
194 ,fRangeFitPRF(c.fRangeFitPRF)
196 ,fMinEntries(c.fMinEntries)
198 ,fNumberFit(c.fNumberFit)
199 ,fNumberFitSuccess(c.fNumberFitSuccess)
200 ,fNumberEnt(c.fNumberEnt)
201 ,fStatisticMean(c.fStatisticMean)
203 ,fDebugLevel(c.fDebugLevel)
204 ,fFitVoir(c.fFitVoir)
205 ,fMagneticField(c.fMagneticField)
207 ,fCurrentCoefE(c.fCurrentCoefE)
208 ,fCurrentCoefE2(c.fCurrentCoefE2)
211 ,fScaleFitFactor(c.fScaleFitFactor)
212 ,fEntriesCurrent(c.fEntriesCurrent)
213 ,fCountDet(c.fCountDet)
219 ,fCurrentCoefDetector(0x0)
220 ,fCurrentCoefDetector2(0x0)
228 if(c.fCalibraMode) fCalibraMode = new AliTRDCalibraMode(*c.fCalibraMode);
230 //Current variables initialised
231 for (Int_t k = 0; k < 2; k++) {
232 fCurrentCoef[k] = 0.0;
233 fCurrentCoef2[k] = 0.0;
235 for (Int_t i = 0; i < 3; i++) {
239 if(c.fCalDet) fCalDet = new AliTRDCalDet(*c.fCalDet);
240 if(c.fCalDet2) fCalDet2 = new AliTRDCalDet(*c.fCalDet2);
242 if(c.fCalROC) fCalROC = new AliTRDCalROC(*c.fCalROC);
243 if(c.fCalROC2) fCalROC = new AliTRDCalROC(*c.fCalROC2);
245 fVectorFit.SetName(c.fVectorFit.GetName());
246 for(Int_t k = 0; k < c.fVectorFit.GetEntriesFast(); k++){
247 AliTRDFitInfo *fitInfo = new AliTRDFitInfo();
248 Int_t detector = ((AliTRDFitInfo *)c.fVectorFit.UncheckedAt(k))->GetDetector();
250 if (GetStack(detector) == 2) {
256 Float_t *coef = new Float_t[ntotal];
257 for (Int_t i = 0; i < ntotal; i++) {
258 coef[i] = ((AliTRDFitInfo *)c.fVectorFit.UncheckedAt(k))->GetCoef()[i];
260 fitInfo->SetCoef(coef);
261 fitInfo->SetDetector(detector);
262 fVectorFit.Add((TObject *) fitInfo);
264 fVectorFit.SetName(c.fVectorFit.GetName());
265 for(Int_t k = 0; k < c.fVectorFit2.GetEntriesFast(); k++){
266 AliTRDFitInfo *fitInfo = new AliTRDFitInfo();
267 Int_t detector = ((AliTRDFitInfo *)c.fVectorFit2.UncheckedAt(k))->GetDetector();
269 if (GetStack(detector) == 2) {
275 Float_t *coef = new Float_t[ntotal];
276 for (Int_t i = 0; i < ntotal; i++) {
277 coef[i] = ((AliTRDFitInfo *)c.fVectorFit2.UncheckedAt(k))->GetCoef()[i];
279 fitInfo->SetCoef(coef);
280 fitInfo->SetDetector(detector);
281 fVectorFit2.Add((TObject *) fitInfo);
286 fGeo = new AliTRDgeometry();
289 //____________________________________________________________________________________
290 AliTRDCalibraFit::~AliTRDCalibraFit()
293 // AliTRDCalibraFit destructor
295 if ( fDebugStreamer ) delete fDebugStreamer;
296 if ( fCalDet ) delete fCalDet;
297 if ( fCalDet2 ) delete fCalDet2;
298 if ( fCalROC ) delete fCalROC;
299 if ( fCalROC2 ) delete fCalROC2;
300 if( fCurrentCoefDetector ) delete [] fCurrentCoefDetector;
301 if( fCurrentCoefDetector2 ) delete [] fCurrentCoefDetector2;
303 fVectorFit2.Delete();
309 //_____________________________________________________________________________
310 void AliTRDCalibraFit::Destroy()
322 //__________________________________________________________________________________
323 void AliTRDCalibraFit::RangeChargeIntegration(Float_t vdrift, Float_t t0, Int_t &begin, Int_t &peak, Int_t &end)
326 // From the drift velocity and t0
327 // return the position of the peak and maximum negative slope
330 const Float_t kDrWidth = AliTRDgeometry::DrThick(); // drift region
331 Double_t widbins = 0.1; // 0.1 mus
333 //peak and maxnegslope in mus
334 Double_t begind = t0*widbins + fT0Shift0;
335 Double_t peakd = t0*widbins + fT0Shift1;
336 Double_t maxnegslope = (kDrWidth + vdrift*peakd)/vdrift;
338 // peak and maxnegslope in timebin
339 begin = TMath::Nint(begind*widbins);
340 peak = TMath::Nint(peakd*widbins);
341 end = TMath::Nint(maxnegslope*widbins);
344 //____________Functions fit Online CH2d________________________________________
345 Bool_t AliTRDCalibraFit::AnalyseCH(TH2I *ch)
348 // Fit the 1D histos, projections of the 2D ch on the Xaxis, for each
349 // calibration group normalized the resulted coefficients (to 1 normally)
352 // Set the calibration mode
353 const char *name = ch->GetTitle();
354 SetModeCalibration(name,0);
356 // Number of Ybins (detectors or groups of pads)
357 Int_t nbins = ch->GetNbinsX();// charge
358 Int_t nybins = ch->GetNbinsY();// groups number
359 if (!InitFit(nybins,0)) {
365 fStatisticMean = 0.0;
367 fNumberFitSuccess = 0;
369 // Init fCountDet and fCount
370 InitfCountDetAndfCount(0);
371 // Beginning of the loop betwwen dect1 and dect2
372 for (Int_t idect = fDect1; idect < fDect2; idect++) {
373 // Determination of fNnZ, fNnRphi, fNfragZ and fNfragRphi...
374 UpdatefCountDetAndfCount(idect,0);
375 ReconstructFitRowMinRowMax(idect, 0);
377 TH1I *projch = (TH1I *) ch->ProjectionX("projch",idect+1,idect+1,(Option_t *)"e");
378 projch->SetDirectory(0);
379 // Number of entries for this calibration group
380 Double_t nentries = 0.0;
382 for (Int_t k = 0; k < nbins; k++) {
383 Int_t binnb = (nbins+2)*(idect+1)+(k+1);
384 nentries += ch->GetBinContent(binnb);
385 mean += projch->GetBinCenter(k+1)*projch->GetBinContent(k+1);
386 projch->SetBinError(k+1,TMath::Sqrt(projch->GetBinContent(k+1)));
388 projch->SetEntries(nentries);
389 //printf("The number of entries for the group %d is %f\n",idect,nentries);
394 // Rebin and statistic stuff
396 projch = ReBin((TH1I *) projch);
398 // This detector has not enough statistics or was off
399 if (nentries <= fMinEntries) {
400 NotEnoughStatisticCH(idect);
401 if (fDebugLevel != 1) {
406 // Statistics of the group fitted
407 fStatisticMean += nentries;
412 case 0: FitMeanW((TH1 *) projch, nentries); break;
413 case 1: FitMean((TH1 *) projch, nentries, mean); break;
414 case 2: FitCH((TH1 *) projch, mean); break;
415 case 3: FitBisCH((TH1 *) projch, mean); break;
416 default: return kFALSE;
419 FillInfosFitCH(idect);
421 if (fDebugLevel != 1) {
426 if (fDebugLevel != 1) {
430 if (fNumberFit > 0) {
431 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));
432 fStatisticMean = fStatisticMean / fNumberFit;
435 AliInfo(Form("There are %d with at least one entries. There is no fit!",fNumberEnt));
437 delete fDebugStreamer;
438 fDebugStreamer = 0x0;
442 //____________Functions fit Online CH2d________________________________________
443 Bool_t AliTRDCalibraFit::AnalyseCH(AliTRDCalibraVector *calvect)
446 // Reconstruct a 1D histo from the vectorCH for each calibration group,
447 // fit the histo, normalized the resulted coefficients (to 1 normally)
450 // Set the calibraMode
451 const char *name = calvect->GetNameCH();
452 SetModeCalibration(name,0);
454 // Number of Xbins (detectors or groups of pads)
455 if (!InitFit((432*calvect->GetDetCha0(0)+108*calvect->GetDetCha2(0)),0)) {
461 fStatisticMean = 0.0;
463 fNumberFitSuccess = 0;
465 // Init fCountDet and fCount
466 InitfCountDetAndfCount(0);
467 // Beginning of the loop between dect1 and dect2
468 for (Int_t idect = fDect1; idect < fDect2; idect++) {
469 // Determination of fNnZ, fNnRphi, fNfragZ and fNfragRphi...........
470 UpdatefCountDetAndfCount(idect,0);
471 ReconstructFitRowMinRowMax(idect,0);
473 Double_t nentries = 0.0;
476 Bool_t something = kTRUE;
477 if(!calvect->GetCHEntries(fCountDet)) something = kFALSE;
481 projch = calvect->ConvertVectorCHHisto(fCountDet,(idect-(fCount-(fCalibraMode->GetNfragZ(0)*fCalibraMode->GetNfragRphi(0)))),(const char *) tname);
482 projch->SetDirectory(0);
483 for (Int_t k = 0; k < calvect->GetNumberBinCharge(); k++) {
484 nentries += projch->GetBinContent(k+1);
485 mean += projch->GetBinCenter(k+1)*projch->GetBinContent(k+1);
491 //printf("The number of entries for the group %d is %f\n",idect,nentries);
494 projch = ReBin((TH1F *) projch);
497 // This detector has not enough statistics or was not found in VectorCH
498 if (nentries <= fMinEntries) {
499 NotEnoughStatisticCH(idect);
500 if (fDebugLevel != 1) {
501 if(projch) delete projch;
505 // Statistic of the histos fitted
506 fStatisticMean += nentries;
511 case 0: FitMeanW((TH1 *) projch, nentries); break;
512 case 1: FitMean((TH1 *) projch, nentries, mean); break;
513 case 2: FitCH((TH1 *) projch, mean); break;
514 case 3: FitBisCH((TH1 *) projch, mean); break;
515 default: return kFALSE;
518 FillInfosFitCH(idect);
520 if (fDebugLevel != 1) {
525 if (fDebugLevel != 1) {
529 if (fNumberFit > 0) {
530 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));
531 fStatisticMean = fStatisticMean / fNumberFit;
534 AliInfo(Form("There are %d with at least one entries. There is no fit!",fNumberEnt));
536 delete fDebugStreamer;
537 fDebugStreamer = 0x0;
540 //________________functions fit Online PH2d____________________________________
541 Bool_t AliTRDCalibraFit::AnalysePH(TProfile2D *ph)
544 // Take the 1D profiles (average pulse height), projections of the 2D PH
545 // on the Xaxis, for each calibration group
546 // Reconstruct a drift velocity
547 // A first calibration of T0 is also made using the same method
550 // Set the calibration mode
551 const char *name = ph->GetTitle();
552 SetModeCalibration(name,1);
554 // Number of Xbins (detectors or groups of pads)
555 Int_t nbins = ph->GetNbinsX();// time
556 Int_t nybins = ph->GetNbinsY();// calibration group
557 if (!InitFit(nybins,1)) {
563 fStatisticMean = 0.0;
565 fNumberFitSuccess = 0;
567 // Init fCountDet and fCount
568 InitfCountDetAndfCount(1);
569 // Beginning of the loop
570 for (Int_t idect = fDect1; idect < fDect2; idect++) {
571 // Determination of fNnZ, fNnRphi, fNfragZ and fNfragRphi.......
572 UpdatefCountDetAndfCount(idect,1);
573 ReconstructFitRowMinRowMax(idect,1);
575 TH1D *projph = (TH1D *) ph->ProjectionX("projph",idect+1,idect+1,(Option_t *) "e");
576 projph->SetDirectory(0);
577 // Number of entries for this calibration group
578 Double_t nentries = 0;
579 for (Int_t k = 0; k < nbins; k++) {
580 Int_t binnb = (nbins+2)*(idect+1)+(k+1);
581 nentries += ph->GetBinEntries(binnb);
586 //printf("The number of entries for the group %d is %f\n",idect,nentries);
587 // This detector has not enough statistics or was off
588 if (nentries <= fMinEntries) {
589 //printf("Not enough statistic!\n");
590 NotEnoughStatisticPH(idect);
591 if (fDebugLevel != 1) {
596 // Statistics of the histos fitted
598 fStatisticMean += nentries;
599 // Calcul of "real" coef
600 CalculVdriftCoefMean();
605 case 0: FitLagrangePoly((TH1 *) projph); break;
606 case 1: FitPente((TH1 *) projph); break;
607 case 2: FitPH((TH1 *) projph,(Int_t) (idect - fDect1)); break;
608 default: return kFALSE;
610 // Fill the tree if end of a detector or only the pointer to the branch!!!
611 FillInfosFitPH(idect);
613 if (fDebugLevel != 1) {
618 if (fNumberFit > 0) {
619 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));
620 fStatisticMean = fStatisticMean / fNumberFit;
623 AliInfo(Form("There are %d with at least one entries. There is no fit!",fNumberEnt));
625 delete fDebugStreamer;
626 fDebugStreamer = 0x0;
629 //____________Functions fit Online PH2d________________________________________
630 Bool_t AliTRDCalibraFit::AnalysePH(AliTRDCalibraVector *calvect)
633 // Reconstruct the average pulse height from the vectorPH for each
635 // Reconstruct a drift velocity
636 // A first calibration of T0 is also made using the same method (slope method)
639 // Set the calibration mode
640 const char *name = calvect->GetNamePH();
641 SetModeCalibration(name,1);
643 // Number of Xbins (detectors or groups of pads)
644 if (!InitFit((432*calvect->GetDetCha0(1)+108*calvect->GetDetCha2(1)),1)) {
650 fStatisticMean = 0.0;
652 fNumberFitSuccess = 0;
654 // Init fCountDet and fCount
655 InitfCountDetAndfCount(1);
656 // Beginning of the loop
657 for (Int_t idect = fDect1; idect < fDect2; idect++) {
658 // Determination of fNnZ, fNnRphi, fNfragZ and fNfragRphi...........
659 UpdatefCountDetAndfCount(idect,1);
660 ReconstructFitRowMinRowMax(idect,1);
664 Bool_t something = kTRUE;
665 if(!calvect->GetPHEntries(fCountDet)) something = kFALSE;
669 projph = CorrectTheError((TGraphErrors *) (calvect->ConvertVectorPHTGraphErrors(fCountDet,(idect-(fCount-(fCalibraMode->GetNfragZ(1)*fCalibraMode->GetNfragRphi(1)))),(const char *) tname)));
670 projph->SetDirectory(0);
672 //printf("The number of entries for the group %d is %d\n",idect,fEntriesCurrent);
673 // This detector has not enough statistics or was off
674 if (fEntriesCurrent <= fMinEntries) {
675 //printf("Not enough stat!\n");
676 NotEnoughStatisticPH(idect);
677 if (fDebugLevel != 1) {
678 if(projph) delete projph;
682 // Statistic of the histos fitted
684 fStatisticMean += fEntriesCurrent;
685 // Calcul of "real" coef
686 CalculVdriftCoefMean();
691 case 0: FitLagrangePoly((TH1 *) projph); break;
692 case 1: FitPente((TH1 *) projph); break;
693 case 2: FitPH((TH1 *) projph,(Int_t) (idect - fDect1)); break;
694 default: return kFALSE;
696 // Fill the tree if end of a detector or only the pointer to the branch!!!
697 FillInfosFitPH(idect);
699 if (fDebugLevel != 1) {
705 if (fNumberFit > 0) {
706 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));
707 fStatisticMean = fStatisticMean / fNumberFit;
710 AliInfo(Form("There are %d with at least one entries. There is no fit!",fNumberEnt));
712 delete fDebugStreamer;
713 fDebugStreamer = 0x0;
716 //____________Functions fit Online PRF2d_______________________________________
717 Bool_t AliTRDCalibraFit::AnalysePRF(TProfile2D *prf)
720 // Take the 1D profiles (pad response function), projections of the 2D PRF
721 // on the Xaxis, for each calibration group
722 // Fit with a gaussian to reconstruct the sigma of the pad response function
725 // Set the calibration mode
726 const char *name = prf->GetTitle();
727 SetModeCalibration(name,2);
729 // Number of Ybins (detectors or groups of pads)
730 Int_t nybins = prf->GetNbinsY();// calibration groups
731 Int_t nbins = prf->GetNbinsX();// bins
732 Int_t nbg = GetNumberOfGroupsPRF((const char *)prf->GetTitle());
733 if((nbg > 0) || (nbg == -1)) return kFALSE;
734 if (!InitFit(nybins,2)) {
740 fStatisticMean = 0.0;
742 fNumberFitSuccess = 0;
744 // Init fCountDet and fCount
745 InitfCountDetAndfCount(2);
746 // Beginning of the loop
747 for (Int_t idect = fDect1; idect < fDect2; idect++) {
748 // Determination of fNnZ, fNnRphi, fNfragZ and fNfragRphi......
749 UpdatefCountDetAndfCount(idect,2);
750 ReconstructFitRowMinRowMax(idect,2);
752 TH1D *projprf = (TH1D *) prf->ProjectionX("projprf",idect+1,idect+1,(Option_t *) "e");
753 projprf->SetDirectory(0);
754 // Number of entries for this calibration group
755 Double_t nentries = 0;
756 for (Int_t k = 0; k < nbins; k++) {
757 Int_t binnb = (nbins+2)*(idect+1)+(k+1);
758 nentries += prf->GetBinEntries(binnb);
760 if(nentries > 0) fNumberEnt++;
761 // This detector has not enough statistics or was off
762 if (nentries <= fMinEntries) {
763 NotEnoughStatisticPRF(idect);
764 if (fDebugLevel != 1) {
769 // Statistics of the histos fitted
771 fStatisticMean += nentries;
772 // Calcul of "real" coef
777 case 0: FitPRF((TH1 *) projprf); break;
778 case 1: RmsPRF((TH1 *) projprf); break;
779 default: return kFALSE;
781 // Fill the tree if end of a detector or only the pointer to the branch!!!
782 FillInfosFitPRF(idect);
784 if (fDebugLevel != 1) {
789 if (fNumberFit > 0) {
790 AliInfo(Form("There are %d with at least one entries.",fNumberEnt));
791 AliInfo(Form("%d fits have been proceeded (sucessfully or not...).",fNumberFit));
792 AliInfo(Form("There is a mean statistic of: %d over these fitted histograms and %d successfulled fits"
793 ,(Int_t) fStatisticMean/fNumberFit,fNumberFitSuccess));
794 fStatisticMean = fStatisticMean / fNumberFit;
797 AliInfo(Form("There are %d with at least one entries. There is no fit!",fNumberEnt));
799 delete fDebugStreamer;
800 fDebugStreamer = 0x0;
803 //____________Functions fit Online PRF2d_______________________________________
804 Bool_t AliTRDCalibraFit::AnalysePRFMarianFit(TProfile2D *prf)
807 // Take the 1D profiles (pad response function), projections of the 2D PRF
808 // on the Xaxis, for each calibration group
809 // Fit with a gaussian to reconstruct the sigma of the pad response function
812 // Set the calibration mode
813 const char *name = prf->GetTitle();
814 SetModeCalibration(name,2);
816 // Number of Ybins (detectors or groups of pads)
817 TAxis *xprf = prf->GetXaxis();
818 TAxis *yprf = prf->GetYaxis();
819 Int_t nybins = yprf->GetNbins();// calibration groups
820 Int_t nbins = xprf->GetNbins();// bins
821 Float_t lowedge = (Float_t) xprf->GetBinLowEdge(1);//lowedge in bins
822 Float_t upedge = (Float_t) xprf->GetBinUpEdge(nbins);//upedge in bins
823 Int_t nbg = GetNumberOfGroupsPRF((const char *)name);
824 if(nbg == -1) return kFALSE;
825 if(nbg > 0) fMethod = 1;
827 if (!InitFit(nybins,2)) {
833 fStatisticMean = 0.0;
835 fNumberFitSuccess = 0;
837 // Init fCountDet and fCount
838 InitfCountDetAndfCount(2);
839 // Beginning of the loop
840 for (Int_t idect = fDect1; idect < fDect2; idect++) {
841 // Determination of fNnZ, fNnRphi, fNfragZ and fNfragRphi.......
842 UpdatefCountDetAndfCount(idect,2);
843 ReconstructFitRowMinRowMax(idect,2);
844 // Build the array of entries and sum
845 TArrayD arraye = TArrayD(nbins);
846 TArrayD arraym = TArrayD(nbins);
847 TArrayD arrayme = TArrayD(nbins);
848 Double_t nentries = 0;
849 //printf("nbins %d\n",nbins);
850 for (Int_t k = 0; k < nbins; k++) {
851 Int_t binnb = (nbins+2)*(idect+1)+(k+1);
852 Double_t entries = (Double_t)prf->GetBinEntries(binnb);
853 Double_t mean = (Double_t)prf->GetBinContent(binnb);
854 Double_t error = (Double_t)prf->GetBinError(binnb);
855 //printf("for %d we have %f\n",k,entries);
857 arraye.AddAt(entries,k);
858 arraym.AddAt(mean,k);
859 arrayme.AddAt(error,k);
861 if(nentries > 0) fNumberEnt++;
862 //printf("The number of entries for the group %d is %f\n",idect,nentries);
863 // This detector has not enough statistics or was off
864 if (nentries <= fMinEntries) {
865 NotEnoughStatisticPRF(idect);
868 // Statistics of the histos fitted
870 fStatisticMean += nentries;
871 // Calcul of "real" coef
876 case 0: FitPRFGausMI( arraye.GetArray(), arraym.GetArray(), arrayme.GetArray(), nbins, lowedge, upedge); break;
877 case 1: FitTnpRange( arraye.GetArray(), arraym.GetArray(), arrayme.GetArray(), nbg, nbins); break;
878 default: return kFALSE;
880 // Fill the tree if end of a detector or only the pointer to the branch!!!
881 FillInfosFitPRF(idect);
884 if (fNumberFit > 0) {
885 AliInfo(Form("There are %d with at least one entries.",fNumberEnt));
886 AliInfo(Form("%d fits have been proceeded (sucessfully or not...).",fNumberFit));
887 AliInfo(Form("There is a mean statistic of: %d over these fitted histograms and %d successfulled fits"
888 ,(Int_t) fStatisticMean/fNumberFit,fNumberFitSuccess));
889 fStatisticMean = fStatisticMean / fNumberFit;
892 AliInfo(Form("There are %d with at least one entries. There is no fit!",fNumberEnt));
894 delete fDebugStreamer;
895 fDebugStreamer = 0x0;
898 //____________Functions fit Online PRF2d_______________________________________
899 Bool_t AliTRDCalibraFit::AnalysePRF(AliTRDCalibraVector *calvect)
902 // Reconstruct the 1D histo (pad response function) from the vectorPRD for
903 // each calibration group
904 // Fit with a gaussian to reconstruct the sigma of the pad response function
907 // Set the calibra mode
908 const char *name = calvect->GetNamePRF();
909 SetModeCalibration(name,2);
910 //printf("test0 %s\n",name);
912 // Number of Xbins (detectors or groups of pads)
913 if (!InitFit((432*calvect->GetDetCha0(2)+108*calvect->GetDetCha2(2)),2)) {
918 ///printf("test2\n");
921 fStatisticMean = 0.0;
923 fNumberFitSuccess = 0;
925 // Init fCountDet and fCount
926 InitfCountDetAndfCount(2);
927 // Beginning of the loop
928 for (Int_t idect = fDect1; idect < fDect2; idect++) {
929 // Determination of fNnZ, fNnRphi, fNfragZ and fNfragRphi........
930 UpdatefCountDetAndfCount(idect,2);
931 ReconstructFitRowMinRowMax(idect,2);
935 Bool_t something = kTRUE;
936 if(!calvect->GetPRFEntries(fCountDet)) something = kFALSE;
938 TString tname("PRF");
940 projprf = CorrectTheError((TGraphErrors *) (calvect->ConvertVectorPRFTGraphErrors(fCountDet,(idect-(fCount-(fCalibraMode->GetNfragZ(1)*fCalibraMode->GetNfragRphi(1)))),(const char *) tname)));
941 projprf->SetDirectory(0);
943 // This detector has not enough statistics or was off
944 if (fEntriesCurrent <= fMinEntries) {
945 NotEnoughStatisticPRF(idect);
946 if (fDebugLevel != 1) {
947 if(projprf) delete projprf;
951 // Statistic of the histos fitted
953 fStatisticMean += fEntriesCurrent;
954 // Calcul of "real" coef
959 case 1: FitPRF((TH1 *) projprf); break;
960 case 2: RmsPRF((TH1 *) projprf); break;
961 default: return kFALSE;
963 // Fill the tree if end of a detector or only the pointer to the branch!!!
964 FillInfosFitPRF(idect);
966 if (fDebugLevel != 1) {
971 if (fNumberFit > 0) {
972 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));
975 AliInfo(Form("There are %d with at least one entries. There is no fit!",fNumberEnt));
977 delete fDebugStreamer;
978 fDebugStreamer = 0x0;
981 //____________Functions fit Online PRF2d_______________________________________
982 Bool_t AliTRDCalibraFit::AnalysePRFMarianFit(AliTRDCalibraVector *calvect)
985 // Reconstruct the 1D histo (pad response function) from the vectorPRD for
986 // each calibration group
987 // Fit with a gaussian to reconstruct the sigma of the pad response function
990 // Set the calibra mode
991 const char *name = calvect->GetNamePRF();
992 SetModeCalibration(name,2);
993 //printf("test0 %s\n",name);
994 Int_t nbg = GetNumberOfGroupsPRF((const char *)name);
995 printf("test1 %d\n",nbg);
996 if(nbg == -1) return kFALSE;
997 if(nbg > 0) fMethod = 1;
999 // Number of Xbins (detectors or groups of pads)
1000 if (!InitFit((432*calvect->GetDetCha0(2)+108*calvect->GetDetCha2(2)),2)) {
1001 //printf("test2\n");
1004 if (!InitFitPRF()) {
1005 //printf("test3\n");
1008 fStatisticMean = 0.0;
1010 fNumberFitSuccess = 0;
1014 Double_t *arrayx = 0;
1015 Double_t *arraye = 0;
1016 Double_t *arraym = 0;
1017 Double_t *arrayme = 0;
1018 Float_t lowedge = 0.0;
1019 Float_t upedge = 0.0;
1020 // Init fCountDet and fCount
1021 InitfCountDetAndfCount(2);
1022 // Beginning of the loop
1023 for (Int_t idect = fDect1; idect < fDect2; idect++) {
1024 // Determination of fNnZ, fNnRphi, fNfragZ and fNfragRphi......
1025 UpdatefCountDetAndfCount(idect,2);
1026 ReconstructFitRowMinRowMax(idect,2);
1028 TGraphErrors *projprftree = 0x0;
1029 fEntriesCurrent = 0;
1030 Bool_t something = kTRUE;
1031 if(!calvect->GetPRFEntries(fCountDet)) something = kFALSE;
1033 TString tname("PRF");
1035 projprftree = calvect->ConvertVectorPRFTGraphErrors(fCountDet,(idect-(fCount-(fCalibraMode->GetNfragZ(1)*fCalibraMode->GetNfragRphi(1)))),(const char *) tname);
1036 nbins = projprftree->GetN();
1037 arrayx = (Double_t *)projprftree->GetX();
1038 arraye = (Double_t *)projprftree->GetEX();
1039 arraym = (Double_t *)projprftree->GetY();
1040 arrayme = (Double_t *)projprftree->GetEY();
1041 Float_t step = arrayx[1]-arrayx[0];
1042 lowedge = arrayx[0] - step/2.0;
1043 upedge = arrayx[(nbins-1)] + step/2.0;
1044 //printf("nbins est %d\n",nbins);
1045 for(Int_t k = 0; k < nbins; k++){
1046 fEntriesCurrent += (Int_t)arraye[k];
1047 //printf("for %d we have %f, %f\n",k,arraye[k],((projprftree->GetEX())[k]));
1048 if(arraye[k]>0.0) arrayme[k] = TMath::Sqrt(TMath::Abs(arrayme[k]-arraym[k]*arraym[k])/arraye[k]);
1050 if(fEntriesCurrent > 0) fNumberEnt++;
1052 //printf("The number of entries for the group %d is %d\n",idect,fEntriesCurrent);
1053 // This detector has not enough statistics or was off
1054 if (fEntriesCurrent <= fMinEntries) {
1055 NotEnoughStatisticPRF(idect);
1056 if(projprftree) delete projprftree;
1059 // Statistic of the histos fitted
1061 fStatisticMean += fEntriesCurrent;
1062 // Calcul of "real" coef
1063 CalculPRFCoefMean();
1067 case 0: FitPRFGausMI(arraye,arraym,arrayme,nbins,lowedge,upedge); break;
1068 case 1: FitTnpRange(arraye,arraym,arrayme,nbg,nbins); break;
1069 default: return kFALSE;
1071 // Fill the tree if end of a detector or only the pointer to the branch!!!
1072 FillInfosFitPRF(idect);
1074 if (fDebugLevel != 1) {
1079 if (fNumberFit > 0) {
1080 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));
1083 AliInfo(Form("There are %d with at least one entries. There is no fit!",fNumberEnt));
1085 delete fDebugStreamer;
1086 fDebugStreamer = 0x0;
1089 //____________Functions fit Online CH2d________________________________________
1090 Bool_t AliTRDCalibraFit::AnalyseLinearFitters(AliTRDCalibraVdriftLinearFit *calivdli)
1093 // The linear method
1096 fStatisticMean = 0.0;
1098 fNumberFitSuccess = 0;
1100 if(!InitFitLinearFitter()) return kFALSE;
1103 for(Int_t idet = 0; idet < 540; idet++){
1106 //printf("detector number %d\n",idet);
1111 fEntriesCurrent = 0;
1113 Bool_t here = calivdli->GetParam(idet,¶m);
1114 Bool_t heree = calivdli->GetError(idet,&error);
1115 //printf("here %d and heree %d\n",here, heree);
1117 fEntriesCurrent = (Int_t) error[2];
1120 //printf("Number of entries %d\n",fEntriesCurrent);
1121 // Nothing found or not enough statistic
1122 if((!heree) || (!here) || (fEntriesCurrent <= fMinEntries)) {
1123 NotEnoughStatisticLinearFitter();
1130 fStatisticMean += fEntriesCurrent;
1133 if((-(param[1])) <= 0.0) {
1134 NotEnoughStatisticLinearFitter();
1138 // CalculDatabaseVdriftandTan
1139 CalculVdriftLorentzCoef();
1142 fNumberFitSuccess ++;
1144 // Put the fCurrentCoef
1145 fCurrentCoef[0] = -param[1];
1146 // here the database must be the one of the reconstruction for the lorentz angle....
1147 fCurrentCoef2[0] = (param[0]+fCurrentCoef[1]*fCurrentCoef2[1])/fCurrentCoef[0];
1148 fCurrentCoefE = error[1];
1149 fCurrentCoefE2 = error[0];
1150 if((fCurrentCoef2[0] != 0.0) && (param[0] != 0.0)){
1151 fCurrentCoefE2 = (fCurrentCoefE2/param[0]+fCurrentCoefE/fCurrentCoef[0])*fCurrentCoef2[0];
1155 FillInfosFitLinearFitter();
1160 if (fNumberFit > 0) {
1161 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));
1164 AliInfo(Form("There are %d with at least one entries. There is no fit!",fNumberEnt));
1166 delete fDebugStreamer;
1167 fDebugStreamer = 0x0;
1171 //____________Functions for seeing if the pad is really okey___________________
1172 //_____________________________________________________________________________
1173 Int_t AliTRDCalibraFit::GetNumberOfGroupsPRF(const char* nametitle)
1176 // Get numberofgroupsprf
1180 const Char_t *pattern0 = "Ngp0";
1181 const Char_t *pattern1 = "Ngp1";
1182 const Char_t *pattern2 = "Ngp2";
1183 const Char_t *pattern3 = "Ngp3";
1184 const Char_t *pattern4 = "Ngp4";
1185 const Char_t *pattern5 = "Ngp5";
1186 const Char_t *pattern6 = "Ngp6";
1189 if (strstr(nametitle,pattern0)) {
1192 if (strstr(nametitle,pattern1)) {
1195 if (strstr(nametitle,pattern2)) {
1198 if (strstr(nametitle,pattern3)) {
1201 if (strstr(nametitle,pattern4)) {
1204 if (strstr(nametitle,pattern5)) {
1207 if (strstr(nametitle,pattern6)){
1214 //_____________________________________________________________________________
1215 Bool_t AliTRDCalibraFit::SetModeCalibration(const char *name, Int_t i)
1218 // Set fNz[i] and fNrphi[i] of the AliTRDCalibraFit::Instance()
1219 // corresponding to the given name
1222 if(!SetNzFromTObject(name,i)) return kFALSE;
1223 if(!SetNrphiFromTObject(name,i)) return kFALSE;
1228 //_____________________________________________________________________________
1229 Bool_t AliTRDCalibraFit::SetNrphiFromTObject(const char *name, Int_t i)
1232 // Set fNrphi[i] of the AliTRDCalibraFit::Instance()
1233 // corresponding to the given TObject
1237 const Char_t *patternrphi0 = "Nrphi0";
1238 const Char_t *patternrphi1 = "Nrphi1";
1239 const Char_t *patternrphi2 = "Nrphi2";
1240 const Char_t *patternrphi3 = "Nrphi3";
1241 const Char_t *patternrphi4 = "Nrphi4";
1242 const Char_t *patternrphi5 = "Nrphi5";
1243 const Char_t *patternrphi6 = "Nrphi6";
1246 if (strstr(name,patternrphi0)) {
1247 fCalibraMode->SetNrphi(i ,0);
1250 if (strstr(name,patternrphi1)) {
1251 fCalibraMode->SetNrphi(i, 1);
1254 if (strstr(name,patternrphi2)) {
1255 fCalibraMode->SetNrphi(i, 2);
1258 if (strstr(name,patternrphi3)) {
1259 fCalibraMode->SetNrphi(i, 3);
1262 if (strstr(name,patternrphi4)) {
1263 fCalibraMode->SetNrphi(i, 4);
1266 if (strstr(name,patternrphi5)) {
1267 fCalibraMode->SetNrphi(i, 5);
1270 if (strstr(name,patternrphi6)) {
1271 fCalibraMode->SetNrphi(i, 6);
1275 fCalibraMode->SetNrphi(i ,0);
1279 //_____________________________________________________________________________
1280 Bool_t AliTRDCalibraFit::SetNzFromTObject(const char *name, Int_t i)
1283 // Set fNz[i] of the AliTRDCalibraFit::Instance()
1284 // corresponding to the given TObject
1288 const Char_t *patternz0 = "Nz0";
1289 const Char_t *patternz1 = "Nz1";
1290 const Char_t *patternz2 = "Nz2";
1291 const Char_t *patternz3 = "Nz3";
1292 const Char_t *patternz4 = "Nz4";
1294 if (strstr(name,patternz0)) {
1295 fCalibraMode->SetNz(i, 0);
1298 if (strstr(name,patternz1)) {
1299 fCalibraMode->SetNz(i ,1);
1302 if (strstr(name,patternz2)) {
1303 fCalibraMode->SetNz(i ,2);
1306 if (strstr(name,patternz3)) {
1307 fCalibraMode->SetNz(i ,3);
1310 if (strstr(name,patternz4)) {
1311 fCalibraMode->SetNz(i ,4);
1315 fCalibraMode->SetNz(i ,0);
1318 //_____________________________________________________________________________
1319 AliTRDCalDet *AliTRDCalibraFit::CreateDetObjectVdrift(TObjArray *vectorFit, Bool_t perdetector)
1322 // It creates the AliTRDCalDet object from the AliTRDFitInfo
1323 // It takes the mean value of the coefficients per detector
1324 // This object has to be written in the database
1327 // Create the DetObject
1328 AliTRDCalDet *object = new AliTRDCalDet("ChamberVdrift","TRD drift velocities (detector value)");
1330 Int_t loop = (Int_t) vectorFit->GetEntriesFast();
1331 if(loop != 540) AliInfo("The Vector Fit is not complete!");
1332 Int_t detector = -1;
1333 Float_t value = 0.0;
1335 for (Int_t k = 0; k < loop; k++) {
1336 detector = ((AliTRDFitInfo *) vectorFit->At(k))->GetDetector();
1339 mean = TMath::Abs(((AliTRDFitInfo *) vectorFit->At(k))->GetCoef()[0]);
1343 Int_t rowMax = fGeo->GetRowMax(GetLayer(detector),GetStack(detector),GetSector(detector));
1344 Int_t colMax = fGeo->GetColMax(GetLayer(detector));
1345 for (Int_t row = 0; row < rowMax; row++) {
1346 for (Int_t col = 0; col < colMax; col++) {
1347 value = ((AliTRDFitInfo *) vectorFit->At(k))->GetCoef()[(Int_t)(col*rowMax+row)];
1348 mean += TMath::Abs(value);
1352 if(count > 0) mean = mean/count;
1354 object->SetValue(detector,mean);
1359 //_____________________________________________________________________________
1360 AliTRDCalDet *AliTRDCalibraFit::CreateDetObjectGain(TObjArray *vectorFit, Double_t scaleFitFactor, Bool_t perdetector)
1363 // It creates the AliTRDCalDet object from the AliTRDFitInfo
1364 // It takes the mean value of the coefficients per detector
1365 // This object has to be written in the database
1368 // Create the DetObject
1369 AliTRDCalDet *object = new AliTRDCalDet("ChamberGainFactor","GainFactor (detector value)");
1372 Int_t loop = (Int_t) vectorFit->GetEntriesFast();
1373 if(loop != 540) AliInfo("The Vector Fit is not complete!");
1374 Int_t detector = -1;
1375 Float_t value = 0.0;
1377 for (Int_t k = 0; k < loop; k++) {
1378 detector = ((AliTRDFitInfo *) vectorFit->At(k))->GetDetector();
1381 value = ((AliTRDFitInfo *) vectorFit->At(k))->GetCoef()[0];
1382 if(value > 0) value = value*scaleFitFactor;
1383 mean = TMath::Abs(value);
1387 Int_t rowMax = fGeo->GetRowMax(GetLayer(detector),GetStack(detector),GetSector(detector));
1388 Int_t colMax = fGeo->GetColMax(GetLayer(detector));
1389 for (Int_t row = 0; row < rowMax; row++) {
1390 for (Int_t col = 0; col < colMax; col++) {
1391 value = ((AliTRDFitInfo *) vectorFit->At(k))->GetCoef()[(Int_t)(col*rowMax+row)];
1392 if(value > 0) value = value*scaleFitFactor;
1393 mean += TMath::Abs(value);
1397 if(count > 0) mean = mean/count;
1399 object->SetValue(detector,mean);
1404 //_____________________________________________________________________________
1405 AliTRDCalDet *AliTRDCalibraFit::CreateDetObjectT0(TObjArray *vectorFit, Bool_t perdetector)
1408 // It creates the AliTRDCalDet object from the AliTRDFitInfo2
1409 // It takes the min value of the coefficients per detector
1410 // This object has to be written in the database
1413 // Create the DetObject
1414 AliTRDCalDet *object = new AliTRDCalDet("ChamberT0","T0 (detector value)");
1416 Int_t loop = (Int_t) vectorFit->GetEntriesFast();
1417 if(loop != 540) AliInfo("The Vector Fit is not complete!");
1418 Int_t detector = -1;
1419 Float_t value = 0.0;
1421 for (Int_t k = 0; k < loop; k++) {
1422 detector = ((AliTRDFitInfo *) vectorFit->At(k))->GetDetector();
1423 Float_t min = 100.0;
1425 min = ((AliTRDFitInfo *) vectorFit->At(k))->GetCoef()[0];
1428 Int_t rowMax = fGeo->GetRowMax(GetLayer(detector),GetStack(detector),GetSector(detector));
1429 Int_t colMax = fGeo->GetColMax(GetLayer(detector));
1430 for (Int_t row = 0; row < rowMax; row++) {
1431 for (Int_t col = 0; col < colMax; col++) {
1432 value = ((AliTRDFitInfo *) vectorFit->At(k))->GetCoef()[(Int_t)(col*rowMax+row)];
1433 if(min > value) min = value;
1437 object->SetValue(detector,min);
1443 //_____________________________________________________________________________
1444 AliTRDCalDet *AliTRDCalibraFit::CreateDetObjectLorentzAngle(TObjArray *vectorFit)
1447 // It creates the AliTRDCalDet object from the AliTRDFitInfo2
1448 // It takes the min value of the coefficients per detector
1449 // This object has to be written in the database
1452 // Create the DetObject
1453 AliTRDCalDet *object = new AliTRDCalDet("tan(lorentzangle)","tan(lorentzangle) (detector value)");
1456 Int_t loop = (Int_t) vectorFit->GetEntriesFast();
1457 if(loop != 540) AliInfo("The Vector Fit is not complete!");
1458 Int_t detector = -1;
1459 Float_t value = 0.0;
1461 for (Int_t k = 0; k < loop; k++) {
1462 detector = ((AliTRDFitInfo *) vectorFit->At(k))->GetDetector();
1464 Int_t rowMax = fGeo->GetRowMax(GetLayer(detector),GetStack(detector),GetSector(detector));
1465 Int_t colMax = fGeo->GetColMax(GetLayer(detector));
1466 Float_t min = 100.0;
1467 for (Int_t row = 0; row < rowMax; row++) {
1468 for (Int_t col = 0; col < colMax; col++) {
1469 value = ((AliTRDFitInfo *) fVectorFit2.At(k))->GetCoef()[(Int_t)(col*rowMax+row)];
1470 mean += -TMath::Abs(value);
1474 if(count > 0) mean = mean/count;
1476 value = ((AliTRDFitInfo *) vectorFit->At(k))->GetCoef()[0];
1477 object->SetValue(detector,-TMath::Abs(value));
1483 //_____________________________________________________________________________
1484 TObject *AliTRDCalibraFit::CreatePadObjectGain(TObjArray *vectorFit, Double_t scaleFitFactor, AliTRDCalDet *detobject)
1487 // It Creates the AliTRDCalPad object from AliTRDFitInfo
1488 // You need first to create the object for the detectors,
1489 // where the mean value is put.
1490 // This object has to be written in the database
1493 // Create the DetObject
1494 AliTRDCalPad *object = new AliTRDCalPad("GainFactor","GainFactor (local variations)");
1497 for(Int_t k = 0; k < 540; k++){
1498 AliTRDCalROC *calROC = object->GetCalROC(k);
1499 Int_t nchannels = calROC->GetNchannels();
1500 for(Int_t ch = 0; ch < nchannels; ch++){
1501 calROC->SetValue(ch,1.0);
1507 Int_t loop = (Int_t) vectorFit->GetEntriesFast();
1508 if(loop != 540) AliInfo("The Vector Fit is not complete!");
1509 Int_t detector = -1;
1510 Float_t value = 0.0;
1512 for (Int_t k = 0; k < loop; k++) {
1513 detector = ((AliTRDFitInfo *) vectorFit->At(k))->GetDetector();
1514 AliTRDCalROC *calROC = object->GetCalROC(detector);
1515 Float_t mean = detobject->GetValue(detector);
1516 if(mean == 0) continue;
1517 Int_t rowMax = calROC->GetNrows();
1518 Int_t colMax = calROC->GetNcols();
1519 for (Int_t row = 0; row < rowMax; row++) {
1520 for (Int_t col = 0; col < colMax; col++) {
1521 value = ((AliTRDFitInfo *) vectorFit->At(k))->GetCoef()[(Int_t)(col*rowMax+row)];
1522 if(value > 0) value = value*scaleFitFactor;
1523 calROC->SetValue(col,row,TMath::Abs(value)/mean);
1531 //_____________________________________________________________________________
1532 TObject *AliTRDCalibraFit::CreatePadObjectVdrift(TObjArray *vectorFit, AliTRDCalDet *detobject)
1535 // It Creates the AliTRDCalPad object from AliTRDFitInfo
1536 // You need first to create the object for the detectors,
1537 // where the mean value is put.
1538 // This object has to be written in the database
1541 // Create the DetObject
1542 AliTRDCalPad *object = new AliTRDCalPad("LocalVdrift","TRD drift velocities (local variations)");
1545 for(Int_t k = 0; k < 540; k++){
1546 AliTRDCalROC *calROC = object->GetCalROC(k);
1547 Int_t nchannels = calROC->GetNchannels();
1548 for(Int_t ch = 0; ch < nchannels; ch++){
1549 calROC->SetValue(ch,1.0);
1555 Int_t loop = (Int_t) vectorFit->GetEntriesFast();
1556 if(loop != 540) AliInfo("The Vector Fit is not complete!");
1557 Int_t detector = -1;
1558 Float_t value = 0.0;
1560 for (Int_t k = 0; k < loop; k++) {
1561 detector = ((AliTRDFitInfo *) vectorFit->At(k))->GetDetector();
1562 AliTRDCalROC *calROC = object->GetCalROC(detector);
1563 Float_t mean = detobject->GetValue(detector);
1564 if(mean == 0) continue;
1565 Int_t rowMax = calROC->GetNrows();
1566 Int_t colMax = calROC->GetNcols();
1567 for (Int_t row = 0; row < rowMax; row++) {
1568 for (Int_t col = 0; col < colMax; col++) {
1569 value = ((AliTRDFitInfo *) vectorFit->At(k))->GetCoef()[(Int_t)(col*rowMax+row)];
1570 calROC->SetValue(col,row,TMath::Abs(value)/mean);
1578 //_____________________________________________________________________________
1579 TObject *AliTRDCalibraFit::CreatePadObjectT0(TObjArray *vectorFit, AliTRDCalDet *detobject)
1582 // It Creates the AliTRDCalPad object from AliTRDFitInfo2
1583 // You need first to create the object for the detectors,
1584 // where the mean value is put.
1585 // This object has to be written in the database
1588 // Create the DetObject
1589 AliTRDCalPad *object = new AliTRDCalPad("LocalT0","T0 (local variations)");
1592 for(Int_t k = 0; k < 540; k++){
1593 AliTRDCalROC *calROC = object->GetCalROC(k);
1594 Int_t nchannels = calROC->GetNchannels();
1595 for(Int_t ch = 0; ch < nchannels; ch++){
1596 calROC->SetValue(ch,0.0);
1602 Int_t loop = (Int_t) vectorFit->GetEntriesFast();
1603 if(loop != 540) AliInfo("The Vector Fit is not complete!");
1604 Int_t detector = -1;
1605 Float_t value = 0.0;
1607 for (Int_t k = 0; k < loop; k++) {
1608 detector = ((AliTRDFitInfo *) vectorFit->At(k))->GetDetector();
1609 AliTRDCalROC *calROC = object->GetCalROC(detector);
1610 Float_t min = detobject->GetValue(detector);
1611 Int_t rowMax = calROC->GetNrows();
1612 Int_t colMax = calROC->GetNcols();
1613 for (Int_t row = 0; row < rowMax; row++) {
1614 for (Int_t col = 0; col < colMax; col++) {
1615 value = ((AliTRDFitInfo *) vectorFit->At(k))->GetCoef()[(Int_t)(col*rowMax+row)];
1616 calROC->SetValue(col,row,value-min);
1624 //_____________________________________________________________________________
1625 TObject *AliTRDCalibraFit::CreatePadObjectPRF(TObjArray *vectorFit)
1628 // It Creates the AliTRDCalPad object from AliTRDFitInfo
1629 // This object has to be written in the database
1632 // Create the DetObject
1633 AliTRDCalPad *object = new AliTRDCalPad("PRFWidth","PRFWidth");
1635 Int_t loop = (Int_t) vectorFit->GetEntriesFast();
1636 if(loop != 540) AliInfo("The Vector Fit is not complete!");
1637 Int_t detector = -1;
1638 Float_t value = 0.0;
1640 for (Int_t k = 0; k < loop; k++) {
1641 detector = ((AliTRDFitInfo *) vectorFit->At(k))->GetDetector();
1642 AliTRDCalROC *calROC = object->GetCalROC(detector);
1643 Int_t rowMax = calROC->GetNrows();
1644 Int_t colMax = calROC->GetNcols();
1645 for (Int_t row = 0; row < rowMax; row++) {
1646 for (Int_t col = 0; col < colMax; col++) {
1647 value = ((AliTRDFitInfo *) vectorFit->At(k))->GetCoef()[(Int_t)(col*rowMax+row)];
1648 calROC->SetValue(col,row,TMath::Abs(value));
1656 //_____________________________________________________________________________
1657 AliTRDCalDet *AliTRDCalibraFit::MakeOutliersStatDet(TObjArray *vectorFit, const char *name, Double_t &mean)
1660 // It Creates the AliTRDCalDet object from AliTRDFitInfo
1661 // 0 successful fit 1 not successful fit
1662 // mean is the mean value over the successful fit
1663 // do not use it for t0: no meaning
1666 // Create the CalObject
1667 AliTRDCalDet *object = new AliTRDCalDet(name,name);
1671 Int_t loop = (Int_t) vectorFit->GetEntriesFast();
1673 AliInfo("The Vector Fit is not complete! We initialise all outliers");
1674 for(Int_t k = 0; k < 540; k++){
1675 object->SetValue(k,1.0);
1678 Int_t detector = -1;
1679 Float_t value = 0.0;
1681 for (Int_t k = 0; k < loop; k++) {
1682 detector = ((AliTRDFitInfo *) vectorFit->At(k))->GetDetector();
1683 value = ((AliTRDFitInfo *) vectorFit->At(k))->GetCoef()[0];
1684 if(value <= 0) object->SetValue(detector,1.0);
1686 object->SetValue(detector,0.0);
1691 if(count > 0) mean /= count;
1694 //_____________________________________________________________________________
1695 TObject *AliTRDCalibraFit::MakeOutliersStatPad(TObjArray *vectorFit, const char *name, Double_t &mean)
1698 // It Creates the AliTRDCalPad object from AliTRDFitInfo
1699 // 0 not successful fit 1 successful fit
1700 // mean mean value over the successful fit
1703 // Create the CalObject
1704 AliTRDCalPad *object = new AliTRDCalPad(name,name);
1708 Int_t loop = (Int_t) vectorFit->GetEntriesFast();
1710 AliInfo("The Vector Fit is not complete! We initialise all outliers");
1711 for(Int_t k = 0; k < 540; k++){
1712 AliTRDCalROC *calROC = object->GetCalROC(k);
1713 Int_t nchannels = calROC->GetNchannels();
1714 for(Int_t ch = 0; ch < nchannels; ch++){
1715 calROC->SetValue(ch,1.0);
1719 Int_t detector = -1;
1720 Float_t value = 0.0;
1722 for (Int_t k = 0; k < loop; k++) {
1723 detector = ((AliTRDFitInfo *) vectorFit->At(k))->GetDetector();
1724 AliTRDCalROC *calROC = object->GetCalROC(detector);
1725 Int_t nchannels = calROC->GetNchannels();
1726 for (Int_t ch = 0; ch < nchannels; ch++) {
1727 value = ((AliTRDFitInfo *) vectorFit->At(k))->GetCoef()[ch];
1728 if(value <= 0) calROC->SetValue(ch,1.0);
1730 calROC->SetValue(ch,0.0);
1736 if(count > 0) mean /= count;
1739 //_____________________________________________________________________________
1740 void AliTRDCalibraFit::SetPeriodeFitPH(Int_t periodeFitPH)
1743 // Set FitPH if 1 then each detector will be fitted
1746 if (periodeFitPH > 0) {
1747 fFitPHPeriode = periodeFitPH;
1750 AliInfo("periodeFitPH must be higher than 0!");
1754 //_____________________________________________________________________________
1755 void AliTRDCalibraFit::SetBeginFitCharge(Float_t beginFitCharge)
1758 // The fit of the deposited charge distribution begins at
1759 // histo->Mean()/beginFitCharge
1760 // You can here set beginFitCharge
1763 if (beginFitCharge > 0) {
1764 fBeginFitCharge = beginFitCharge;
1767 AliInfo("beginFitCharge must be strict positif!");
1772 //_____________________________________________________________________________
1773 void AliTRDCalibraFit::SetT0Shift0(Float_t t0Shift)
1776 // The t0 calculated with the maximum positif slope is shift from t0Shift0
1777 // You can here set t0Shift0
1781 fT0Shift0 = t0Shift;
1784 AliInfo("t0Shift0 must be strict positif!");
1789 //_____________________________________________________________________________
1790 void AliTRDCalibraFit::SetT0Shift1(Float_t t0Shift)
1793 // The t0 calculated with the maximum of the amplification region is shift from t0Shift1
1794 // You can here set t0Shift1
1798 fT0Shift1 = t0Shift;
1801 AliInfo("t0Shift must be strict positif!");
1806 //_____________________________________________________________________________
1807 void AliTRDCalibraFit::SetRangeFitPRF(Float_t rangeFitPRF)
1810 // The fit of the PRF is from -rangeFitPRF to rangeFitPRF
1811 // You can here set rangeFitPRF
1814 if ((rangeFitPRF > 0) &&
1815 (rangeFitPRF <= 1.5)) {
1816 fRangeFitPRF = rangeFitPRF;
1819 AliInfo("rangeFitPRF must be between 0 and 1.0");
1824 //_____________________________________________________________________________
1825 void AliTRDCalibraFit::SetMinEntries(Int_t minEntries)
1828 // Minimum entries for fitting
1831 if (minEntries > 0) {
1832 fMinEntries = minEntries;
1835 AliInfo("fMinEntries must be >= 0.");
1840 //_____________________________________________________________________________
1841 void AliTRDCalibraFit::SetRebin(Short_t rebin)
1844 // Rebin with rebin time less bins the Ch histo
1845 // You can set here rebin that should divide the number of bins of CH histo
1850 AliInfo("You have to be sure that fRebin divides fNumberBinCharge used!");
1853 AliInfo("You have to choose a positiv value!");
1857 //_____________________________________________________________________________
1858 Bool_t AliTRDCalibraFit::FillVectorFit()
1861 // For the Fit functions fill the vector Fit
1864 AliTRDFitInfo *fitInfo = new AliTRDFitInfo();
1867 if (GetStack(fCountDet) == 2) {
1874 //printf("For the detector %d , ntotal %d and fCoefCH[0] %f\n",countdet,ntotal,fCoefCH[0]);
1875 Float_t *coef = new Float_t[ntotal];
1876 for (Int_t i = 0; i < ntotal; i++) {
1877 coef[i] = fCurrentCoefDetector[i];
1880 Int_t detector = fCountDet;
1882 fitInfo->SetCoef(coef);
1883 fitInfo->SetDetector(detector);
1884 fVectorFit.Add((TObject *) fitInfo);
1889 //_____________________________________________________________________________
1890 Bool_t AliTRDCalibraFit::FillVectorFit2()
1893 // For the Fit functions fill the vector Fit
1896 AliTRDFitInfo *fitInfo = new AliTRDFitInfo();
1899 if (GetStack(fCountDet) == 2) {
1906 //printf("For the detector %d , ntotal %d and fCoefCH[0] %f\n",countdet,ntotal,fCoefCH[0]);
1907 Float_t *coef = new Float_t[ntotal];
1908 for (Int_t i = 0; i < ntotal; i++) {
1909 coef[i] = fCurrentCoefDetector2[i];
1912 Int_t detector = fCountDet;
1914 fitInfo->SetCoef(coef);
1915 fitInfo->SetDetector(detector);
1916 fVectorFit2.Add((TObject *) fitInfo);
1921 //____________Functions for initialising the AliTRDCalibraFit in the code_________
1922 Bool_t AliTRDCalibraFit::InitFit(Int_t nbins, Int_t i)
1925 // Init the number of expected bins and fDect1[i] fDect2[i]
1928 gStyle->SetPalette(1);
1929 gStyle->SetOptStat(1111);
1930 gStyle->SetPadBorderMode(0);
1931 gStyle->SetCanvasColor(10);
1932 gStyle->SetPadLeftMargin(0.13);
1933 gStyle->SetPadRightMargin(0.01);
1935 // Mode groups of pads: the total number of bins!
1936 CalculNumberOfBinsExpected(i);
1938 // Quick verification that we have the good pad calibration mode!
1939 if (fNumberOfBinsExpected != nbins) {
1940 AliInfo("It doesn't correspond to the mode of pad group calibration!");
1944 // Security for fDebug 3 and 4
1945 if ((fDebugLevel >= 3) &&
1949 AliInfo("This detector doesn't exit!");
1953 // Determine fDet1 and fDet2 and set the fNfragZ and fNfragRphi for debug 3 and 4
1954 CalculDect1Dect2(i);
1959 //____________Functions for initialising the AliTRDCalibraFit in the code_________
1960 Bool_t AliTRDCalibraFit::InitFitCH()
1963 // Init the fVectorFitCH for normalisation
1964 // Init the histo for debugging
1969 fScaleFitFactor = 0.0;
1970 fCurrentCoefDetector = new Float_t[2304];
1971 for (Int_t k = 0; k < 2304; k++) {
1972 fCurrentCoefDetector[k] = 0.0;
1974 fVectorFit.SetName("gainfactorscoefficients");
1976 // fDebug == 0 nothing
1977 // fDebug == 1 and fFitVoir no histo
1978 if (fDebugLevel == 1) {
1979 if(!CheckFitVoir()) return kFALSE;
1981 //Get the CalDet object
1983 AliTRDcalibDB *cal = AliTRDcalibDB::Instance();
1985 AliInfo("Could not get calibDB");
1988 if(fCalDet) delete fCalDet;
1989 fCalDet = new AliTRDCalDet(*(cal->GetGainFactorDet()));
1992 Float_t devalue = 1.0;
1993 if(fCalDet) delete fCalDet;
1994 fCalDet = new AliTRDCalDet("ChamberGainFactor","GainFactor (detector value)");
1995 for(Int_t k = 0; k < 540; k++){
1996 fCalDet->SetValue(k,devalue);
2002 //____________Functions for initialising the AliTRDCalibraFit in the code_________
2003 Bool_t AliTRDCalibraFit::InitFitPH()
2006 // Init the arrays of results
2007 // Init the histos for debugging
2011 fVectorFit.SetName("driftvelocitycoefficients");
2012 fVectorFit2.SetName("t0coefficients");
2014 fCurrentCoefDetector = new Float_t[2304];
2015 for (Int_t k = 0; k < 2304; k++) {
2016 fCurrentCoefDetector[k] = 0.0;
2019 fCurrentCoefDetector2 = new Float_t[2304];
2020 for (Int_t k = 0; k < 2304; k++) {
2021 fCurrentCoefDetector2[k] = 0.0;
2024 //fDebug == 0 nothing
2025 // fDebug == 1 and fFitVoir no histo
2026 if (fDebugLevel == 1) {
2027 if(!CheckFitVoir()) return kFALSE;
2029 //Get the CalDet object
2031 AliTRDcalibDB *cal = AliTRDcalibDB::Instance();
2033 AliInfo("Could not get calibDB");
2036 if(fCalDet) delete fCalDet;
2037 if(fCalDet2) delete fCalDet2;
2038 fCalDet = new AliTRDCalDet(*(cal->GetVdriftDet()));
2039 fCalDet2 = new AliTRDCalDet(*(cal->GetT0Det()));
2042 Float_t devalue = 1.5;
2043 Float_t devalue2 = 0.0;
2044 if(fCalDet) delete fCalDet;
2045 if(fCalDet2) delete fCalDet2;
2046 fCalDet = new AliTRDCalDet("ChamberVdrift","TRD drift velocities (detector value)");
2047 fCalDet2 = new AliTRDCalDet("ChamberT0","T0 (detector value)");
2048 for(Int_t k = 0; k < 540; k++){
2049 fCalDet->SetValue(k,devalue);
2050 fCalDet2->SetValue(k,devalue2);
2055 //____________Functions for initialising the AliTRDCalibraFit in the code_________
2056 Bool_t AliTRDCalibraFit::InitFitPRF()
2059 // Init the calibration mode (Nz, Nrphi), the histograms for
2060 // debugging the fit methods if fDebug > 0,
2064 fVectorFit.SetName("prfwidthcoefficients");
2066 fCurrentCoefDetector = new Float_t[2304];
2067 for (Int_t k = 0; k < 2304; k++) {
2068 fCurrentCoefDetector[k] = 0.0;
2071 // fDebug == 0 nothing
2072 // fDebug == 1 and fFitVoir no histo
2073 if (fDebugLevel == 1) {
2074 if(!CheckFitVoir()) return kFALSE;
2078 //____________Functions for initialising the AliTRDCalibraFit in the code_________
2079 Bool_t AliTRDCalibraFit::InitFitLinearFitter()
2082 // Init the fCalDet, fVectorFit fCurrentCoefDetector
2087 fCurrentCoefDetector = new Float_t[2304];
2088 fCurrentCoefDetector2 = new Float_t[2304];
2089 for (Int_t k = 0; k < 2304; k++) {
2090 fCurrentCoefDetector[k] = 0.0;
2091 fCurrentCoefDetector2[k] = 0.0;
2094 //printf("test0\n");
2096 AliTRDcalibDB *cal = AliTRDcalibDB::Instance();
2098 AliInfo("Could not get calibDB");
2102 //Get the CalDet object
2104 if(fCalDet) delete fCalDet;
2105 if(fCalDet2) delete fCalDet2;
2106 fCalDet = new AliTRDCalDet(*(cal->GetVdriftDet()));
2107 //printf("test1\n");
2108 fCalDet2 = new AliTRDCalDet("lorentz angle tan","lorentz angle tan (detector value)");
2109 //printf("test2\n");
2110 for(Int_t k = 0; k < 540; k++){
2111 fCalDet2->SetValue(k,cal->GetOmegaTau(fCalDet->GetValue(k),-fMagneticField));
2113 //printf("test3\n");
2116 Float_t devalue = 1.5;
2117 Float_t devalue2 = cal->GetOmegaTau(1.5,-fMagneticField);
2118 if(fCalDet) delete fCalDet;
2119 if(fCalDet2) delete fCalDet2;
2120 //printf("test1\n");
2121 fCalDet = new AliTRDCalDet("ChamberVdrift","TRD drift velocities (detector value)");
2122 fCalDet2 = new AliTRDCalDet("lorentz angle tan","lorentz angle tan (detector value)");
2123 //printf("test2\n");
2124 for(Int_t k = 0; k < 540; k++){
2125 fCalDet->SetValue(k,devalue);
2126 fCalDet2->SetValue(k,devalue2);
2128 //printf("test3\n");
2133 //____________Functions for initialising the AliTRDCalibraFit in the code_________
2134 void AliTRDCalibraFit::InitfCountDetAndfCount(Int_t i)
2137 // Init the current detector where we are fCountDet and the
2138 // next fCount for the functions Fit...
2141 // Loop on the Xbins of ch!!
2142 fCountDet = -1; // Current detector
2143 fCount = 0; // To find the next detector
2146 if (fDebugLevel >= 3) {
2147 // Set countdet to the detector
2148 fCountDet = AliTRDgeometry::GetDetector(fDet[0],fDet[1],fDet[2]);
2149 // Set counter to write at the end of the detector
2151 // Get the right calib objects
2154 if(fDebugLevel == 1) {
2156 fCalibraMode->CalculXBins(fCountDet,i);
2157 while(fCalibraMode->GetXbins(i) <=fFitVoir){
2159 fCalibraMode->CalculXBins(fCountDet,i);
2161 fCount = fCalibraMode->GetXbins(i);
2163 // Determination of fNnZ, fNnRphi, fNfragZ and fNfragRphi
2164 fCalibraMode->ModePadCalibration((Int_t) GetStack(fCountDet),i);
2165 fCalibraMode->ModePadFragmentation((Int_t) GetLayer(fCountDet)
2166 ,(Int_t) GetStack(fCountDet)
2167 ,(Int_t) GetSector(fCountDet),i);
2170 //_______________________________________________________________________________
2171 void AliTRDCalibraFit::CalculNumberOfBinsExpected(Int_t i)
2174 // Calculate the number of bins expected (calibration groups)
2177 fNumberOfBinsExpected = 0;
2178 fCalibraMode->ModePadCalibration(2,i);
2179 fCalibraMode->ModePadFragmentation(0,2,0,i);
2180 fCalibraMode->SetDetChamb2(i);
2181 if (fDebugLevel > 1) {
2182 AliInfo(Form("For the chamber 2: %d",fCalibraMode->GetDetChamb2(i)));
2184 fNumberOfBinsExpected += 6 * 18 * fCalibraMode->GetDetChamb2(i);
2185 fCalibraMode->ModePadCalibration(0,i);
2186 fCalibraMode->ModePadFragmentation(0,0,0,i);
2187 fCalibraMode->SetDetChamb0(i);
2188 if (fDebugLevel > 1) {
2189 AliInfo(Form("For the other chamber 0: %d",fCalibraMode->GetDetChamb0(i)));
2191 fNumberOfBinsExpected += 6 * 4 * 18 * fCalibraMode->GetDetChamb0(i);
2194 //_______________________________________________________________________________
2195 void AliTRDCalibraFit::CalculDect1Dect2(Int_t i)
2198 // Calculate the range of fits
2203 if (fDebugLevel == 1) {
2207 if ((fDebugLevel == 2) || (fDebugLevel == 0)) {
2209 fDect2 = fNumberOfBinsExpected;
2211 if (fDebugLevel >= 3) {
2212 fCountDet = AliTRDgeometry::GetDetector(fDet[0],fDet[1],fDet[2]);
2213 fCalibraMode->CalculXBins(fCountDet,i);
2214 fDect1 = fCalibraMode->GetXbins(i);
2215 // Determination of fNnZ, fNnRphi, fNfragZ and fNfragRphi
2216 fCalibraMode->ModePadCalibration((Int_t) GetStack(fCountDet),i);
2217 fCalibraMode->ModePadFragmentation((Int_t) GetLayer(fCountDet)
2218 ,(Int_t) GetStack(fCountDet)
2219 ,(Int_t) GetSector(fCountDet),i);
2220 // Set for the next detector
2221 fDect2 = fDect1 + fCalibraMode->GetNfragZ(i)*fCalibraMode->GetNfragRphi(i);
2224 //_______________________________________________________________________________
2225 Bool_t AliTRDCalibraFit::CheckFitVoir()
2228 // Check if fFitVoir is in the range
2231 if (fFitVoir < fNumberOfBinsExpected) {
2232 AliInfo(Form("We will see the fit of the object %d",fFitVoir));
2235 AliInfo("fFitVoir is out of range of the histo!");
2240 //____________Functions for initialising the AliTRDCalibraFit in the code_________
2241 void AliTRDCalibraFit::UpdatefCountDetAndfCount(Int_t idect, Int_t i)
2244 // See if we are in a new detector and update the
2245 // variables fNfragZ and fNfragRphi if yes
2246 // Will never happen for only one detector (3 and 4)
2247 // Doesn't matter for 2
2249 if (fCount == idect) {
2250 // On en est au detector
2252 // Determination of fNnZ, fNnRphi, fNfragZ and fNfragRphi
2253 fCalibraMode->ModePadCalibration((Int_t) GetStack(fCountDet),i);
2254 fCalibraMode->ModePadFragmentation((Int_t) GetLayer(fCountDet)
2255 ,(Int_t) GetStack(fCountDet)
2256 ,(Int_t) GetSector(fCountDet),i);
2257 // Set for the next detector
2258 fCount += fCalibraMode->GetNfragZ(i)*fCalibraMode->GetNfragRphi(i);
2263 //____________Functions for initialising the AliTRDCalibraFit in the code_________
2264 void AliTRDCalibraFit::ReconstructFitRowMinRowMax(Int_t idect, Int_t i)
2267 // Reconstruct the min pad row, max pad row, min pad col and
2268 // max pad col of the calibration group for the Fit functions
2270 if (fDebugLevel != 1) {
2271 fCalibraMode->ReconstructionRowPadGroup((Int_t) (idect-(fCount-(fCalibraMode->GetNfragZ(i)*fCalibraMode->GetNfragRphi(i)))),i);
2274 //____________Functions for initialising the AliTRDCalibraFit in the code_________
2275 Bool_t AliTRDCalibraFit::NotEnoughStatisticCH(Int_t idect)
2278 // For the case where there are not enough entries in the histograms
2279 // of the calibration group, the value present in the choosen database
2280 // will be put. A negativ sign enables to know that a fit was not possible.
2283 if (fDebugLevel == 1) {
2284 AliInfo("The element has not enough statistic to be fitted");
2289 AliInfo(Form("The element %d in this detector %d has not enough statistic to be fitted"
2290 ,idect-(fCount-(fCalibraMode->GetNfragZ(0)*fCalibraMode->GetNfragRphi(0))),fCountDet));
2292 // Calcul the coef from the database choosen
2293 CalculChargeCoefMean(kFALSE);
2295 //stack 2, not stack 2
2297 if(GetStack(fCountDet) == 2) factor = 12;
2300 // Fill the fCurrentCoefDetector with negative value to say: not fitted
2301 for (Int_t k = fCalibraMode->GetRowMin(0); k < fCalibraMode->GetRowMax(0); k++) {
2302 for (Int_t j = fCalibraMode->GetColMin(0); j < fCalibraMode->GetColMax(0); j++) {
2303 fCurrentCoefDetector[(Int_t)(j*factor+k)] = -TMath::Abs(fCurrentCoef[1]);
2307 //Put default value negative
2308 fCurrentCoef[0] = -TMath::Abs(fCurrentCoef[1]);
2309 fCurrentCoefE = 0.0;
2318 //____________Functions for initialising the AliTRDCalibraFit in the code_________
2319 Bool_t AliTRDCalibraFit::NotEnoughStatisticPH(Int_t idect)
2322 // For the case where there are not enough entries in the histograms
2323 // of the calibration group, the value present in the choosen database
2324 // will be put. A negativ sign enables to know that a fit was not possible.
2326 if (fDebugLevel == 1) {
2327 AliInfo("The element has not enough statistic to be fitted");
2331 AliInfo(Form("The element %d in this detector %d has not enough statistic to be fitted"
2332 ,idect-(fCount-(fCalibraMode->GetNfragZ(1)*fCalibraMode->GetNfragRphi(1))),fCountDet));
2334 CalculVdriftCoefMean();
2337 //stack 2 and not stack 2
2339 if(GetStack(fCountDet) == 2) factor = 12;
2343 // Fill the fCurrentCoefDetector 2
2344 for (Int_t k = fCalibraMode->GetRowMin(1); k < fCalibraMode->GetRowMax(1); k++) {
2345 for (Int_t j = fCalibraMode->GetColMin(1); j < fCalibraMode->GetColMax(1); j++) {
2346 fCurrentCoefDetector[(Int_t)(j*factor+k)] = -TMath::Abs(fCurrentCoef[1]);
2347 fCurrentCoefDetector2[(Int_t)(j*factor+k)] = fCurrentCoef2[1];
2351 // Put the default value
2352 fCurrentCoef[0] = -TMath::Abs(fCurrentCoef[1]);
2353 fCurrentCoefE = 0.0;
2354 fCurrentCoef2[0] = fCurrentCoef2[1];
2355 fCurrentCoefE2 = 0.0;
2366 //____________Functions for initialising the AliTRDCalibraFit in the code_________
2367 Bool_t AliTRDCalibraFit::NotEnoughStatisticPRF(Int_t idect)
2370 // For the case where there are not enough entries in the histograms
2371 // of the calibration group, the value present in the choosen database
2372 // will be put. A negativ sign enables to know that a fit was not possible.
2375 if (fDebugLevel == 1) {
2376 AliInfo("The element has not enough statistic to be fitted");
2380 AliInfo(Form("The element %d in this detector %d has not enough statistic to be fitted"
2381 ,idect-(fCount-(fCalibraMode->GetNfragZ(2)*fCalibraMode->GetNfragRphi(2))),fCountDet));
2383 CalculPRFCoefMean();
2385 // stack 2 and not stack 2
2387 if(GetStack(fCountDet) == 2) factor = 12;
2391 // Fill the fCurrentCoefDetector
2392 for (Int_t k = fCalibraMode->GetRowMin(2); k < fCalibraMode->GetRowMax(2); k++) {
2393 for (Int_t j = fCalibraMode->GetColMin(2); j < fCalibraMode->GetColMax(2); j++) {
2394 fCurrentCoefDetector[(Int_t)(j*factor+k)] = -fCurrentCoef[1];
2398 // Put the default value
2399 fCurrentCoef[0] = -fCurrentCoef[1];
2400 fCurrentCoefE = 0.0;
2408 //____________Functions for initialising the AliTRDCalibraFit in the code_________
2409 Bool_t AliTRDCalibraFit::NotEnoughStatisticLinearFitter()
2412 // For the case where there are not enough entries in the histograms
2413 // of the calibration group, the value present in the choosen database
2414 // will be put. A negativ sign enables to know that a fit was not possible.
2417 // Calcul the coef from the database choosen
2418 CalculVdriftLorentzCoef();
2421 if(GetStack(fCountDet) == 2) factor = 1728;
2425 // Fill the fCurrentCoefDetector
2426 for (Int_t k = 0; k < factor; k++) {
2427 fCurrentCoefDetector[k] = -TMath::Abs(fCurrentCoef[1]);
2428 // should be negative
2429 fCurrentCoefDetector2[k] = +TMath::Abs(fCurrentCoef2[1]);
2433 //Put default opposite sign
2434 fCurrentCoef[0] = -TMath::Abs(fCurrentCoef[1]);
2435 fCurrentCoefE = 0.0;
2436 fCurrentCoef2[0] = +TMath::Abs(fCurrentCoef2[1]);
2437 fCurrentCoefE2 = 0.0;
2439 FillFillLinearFitter();
2444 //____________Functions for initialising the AliTRDCalibraFit in the code_________
2445 Bool_t AliTRDCalibraFit::FillInfosFitCH(Int_t idect)
2448 // Fill the coefficients found with the fits or other
2449 // methods from the Fit functions
2452 if (fDebugLevel != 1) {
2455 if(GetStack(fCountDet) == 2) factor = 12;
2458 for (Int_t k = fCalibraMode->GetRowMin(0); k < fCalibraMode->GetRowMax(0); k++) {
2459 for (Int_t j = fCalibraMode->GetColMin(0); j < fCalibraMode->GetColMax(0); j++) {
2460 fCurrentCoefDetector[(Int_t)(j*factor+k)] = fCurrentCoef[0];
2471 //____________Functions for initialising the AliTRDCalibraFit in the code_________
2472 Bool_t AliTRDCalibraFit::FillInfosFitPH(Int_t idect)
2475 // Fill the coefficients found with the fits or other
2476 // methods from the Fit functions
2479 if (fDebugLevel != 1) {
2482 if(GetStack(fCountDet) == 2) factor = 12;
2485 for (Int_t k = fCalibraMode->GetRowMin(1); k < fCalibraMode->GetRowMax(1); k++) {
2486 for (Int_t j = fCalibraMode->GetColMin(1); j < fCalibraMode->GetColMax(1); j++) {
2487 fCurrentCoefDetector[(Int_t)(j*factor+k)] = fCurrentCoef[0];
2488 fCurrentCoefDetector2[(Int_t)(j*factor+k)] = fCurrentCoef2[0];
2495 //____________Functions for initialising the AliTRDCalibraFit in the code_________
2496 Bool_t AliTRDCalibraFit::FillInfosFitPRF(Int_t idect)
2499 // Fill the coefficients found with the fits or other
2500 // methods from the Fit functions
2503 if (fDebugLevel != 1) {
2506 if(GetStack(fCountDet) == 2) factor = 12;
2509 // Pointer to the branch
2510 for (Int_t k = fCalibraMode->GetRowMin(2); k < fCalibraMode->GetRowMax(2); k++) {
2511 for (Int_t j = fCalibraMode->GetColMin(2); j < fCalibraMode->GetColMax(2); j++) {
2512 fCurrentCoefDetector[(Int_t)(j*factor+k)] = fCurrentCoef[0];
2521 //____________Functions for initialising the AliTRDCalibraFit in the code_________
2522 Bool_t AliTRDCalibraFit::FillInfosFitLinearFitter()
2525 // Fill the coefficients found with the fits or other
2526 // methods from the Fit functions
2530 if(GetStack(fCountDet) == 2) factor = 1728;
2533 // Pointer to the branch
2534 for (Int_t k = 0; k < factor; k++) {
2535 fCurrentCoefDetector[k] = fCurrentCoef[0];
2536 fCurrentCoefDetector2[k] = fCurrentCoef2[0];
2539 FillFillLinearFitter();
2544 //________________________________________________________________________________
2545 void AliTRDCalibraFit::FillFillCH(Int_t idect)
2548 // DebugStream and fVectorFit
2551 // End of one detector
2552 if ((idect == (fCount-1))) {
2555 for (Int_t k = 0; k < 2304; k++) {
2556 fCurrentCoefDetector[k] = 0.0;
2560 if(fDebugLevel > 1){
2562 if ( !fDebugStreamer ) {
2564 TDirectory *backup = gDirectory;
2565 fDebugStreamer = new TTreeSRedirector("TRDDebugFitCH.root");
2566 if ( backup ) backup->cd(); //we don't want to be cd'd to the debug streamer
2569 Int_t detector = fCountDet;
2570 Int_t caligroup = idect;
2571 Short_t rowmin = fCalibraMode->GetRowMin(0);
2572 Short_t rowmax = fCalibraMode->GetRowMax(0);
2573 Short_t colmin = fCalibraMode->GetColMin(0);
2574 Short_t colmax = fCalibraMode->GetColMax(0);
2575 Float_t gf = fCurrentCoef[0];
2576 Float_t gfs = fCurrentCoef[1];
2577 Float_t gfE = fCurrentCoefE;
2579 (*fDebugStreamer) << "FillFillCH" <<
2580 "detector=" << detector <<
2581 "caligroup=" << caligroup <<
2582 "rowmin=" << rowmin <<
2583 "rowmax=" << rowmax <<
2584 "colmin=" << colmin <<
2585 "colmax=" << colmax <<
2593 //________________________________________________________________________________
2594 void AliTRDCalibraFit::FillFillPH(Int_t idect)
2597 // DebugStream and fVectorFit and fVectorFit2
2600 // End of one detector
2601 if ((idect == (fCount-1))) {
2605 for (Int_t k = 0; k < 2304; k++) {
2606 fCurrentCoefDetector[k] = 0.0;
2607 fCurrentCoefDetector2[k] = 0.0;
2611 if(fDebugLevel > 1){
2613 if ( !fDebugStreamer ) {
2615 TDirectory *backup = gDirectory;
2616 fDebugStreamer = new TTreeSRedirector("TRDDebugFitPH.root");
2617 if ( backup ) backup->cd(); //we don't want to be cd'd to the debug streamer
2621 Int_t detector = fCountDet;
2622 Int_t caligroup = idect;
2623 Short_t rowmin = fCalibraMode->GetRowMin(1);
2624 Short_t rowmax = fCalibraMode->GetRowMax(1);
2625 Short_t colmin = fCalibraMode->GetColMin(1);
2626 Short_t colmax = fCalibraMode->GetColMax(1);
2627 Float_t vf = fCurrentCoef[0];
2628 Float_t vs = fCurrentCoef[1];
2629 Float_t vfE = fCurrentCoefE;
2630 Float_t t0f = fCurrentCoef2[0];
2631 Float_t t0s = fCurrentCoef2[1];
2632 Float_t t0E = fCurrentCoefE2;
2636 (* fDebugStreamer) << "FillFillPH"<<
2637 "detector="<<detector<<
2638 "caligroup="<<caligroup<<
2653 //________________________________________________________________________________
2654 void AliTRDCalibraFit::FillFillPRF(Int_t idect)
2657 // DebugStream and fVectorFit
2660 // End of one detector
2661 if ((idect == (fCount-1))) {
2664 for (Int_t k = 0; k < 2304; k++) {
2665 fCurrentCoefDetector[k] = 0.0;
2670 if(fDebugLevel > 1){
2672 if ( !fDebugStreamer ) {
2674 TDirectory *backup = gDirectory;
2675 fDebugStreamer = new TTreeSRedirector("TRDDebugFitPRF.root");
2676 if ( backup ) backup->cd(); //we don't want to be cd'd to the debug streamer
2679 Int_t detector = fCountDet;
2680 Int_t layer = GetLayer(fCountDet);
2681 Int_t caligroup = idect;
2682 Short_t rowmin = fCalibraMode->GetRowMin(2);
2683 Short_t rowmax = fCalibraMode->GetRowMax(2);
2684 Short_t colmin = fCalibraMode->GetColMin(2);
2685 Short_t colmax = fCalibraMode->GetColMax(2);
2686 Float_t widf = fCurrentCoef[0];
2687 Float_t wids = fCurrentCoef[1];
2688 Float_t widfE = fCurrentCoefE;
2690 (* fDebugStreamer) << "FillFillPRF"<<
2691 "detector="<<detector<<
2693 "caligroup="<<caligroup<<
2705 //________________________________________________________________________________
2706 void AliTRDCalibraFit::FillFillLinearFitter()
2709 // DebugStream and fVectorFit
2712 // End of one detector
2718 for (Int_t k = 0; k < 2304; k++) {
2719 fCurrentCoefDetector[k] = 0.0;
2720 fCurrentCoefDetector2[k] = 0.0;
2724 if(fDebugLevel > 1){
2726 if ( !fDebugStreamer ) {
2728 TDirectory *backup = gDirectory;
2729 fDebugStreamer = new TTreeSRedirector("TRDDebugFitLinearFitter.root");
2730 if ( backup ) backup->cd(); //we don't want to be cd'd to the debug streamer
2733 //Debug: comparaison of the different methods (okey for first time but not for iterative procedure)
2734 AliTRDpadPlane *padplane = fGeo->GetPadPlane(GetLayer(fCountDet),GetStack(fCountDet));
2735 Float_t rowmd = (padplane->GetRow0()+padplane->GetRowEnd())/2.;
2736 Float_t r = AliTRDgeometry::GetTime0(GetLayer(fCountDet));
2737 Float_t tiltangle = padplane->GetTiltingAngle();
2738 Int_t detector = fCountDet;
2739 Int_t stack = GetStack(fCountDet);
2740 Int_t layer = GetLayer(fCountDet);
2741 Float_t vf = fCurrentCoef[0];
2742 Float_t vs = fCurrentCoef[1];
2743 Float_t vfE = fCurrentCoefE;
2744 Float_t lorentzangler = fCurrentCoef2[0];
2745 Float_t elorentzangler = fCurrentCoefE2;
2746 Float_t lorentzangles = fCurrentCoef2[1];
2748 (* fDebugStreamer) << "FillFillLinearFitter"<<
2749 "detector="<<detector<<
2754 "tiltangle="<<tiltangle<<
2758 "lorentzangler="<<lorentzangler<<
2759 "Elorentzangler="<<elorentzangler<<
2760 "lorentzangles="<<lorentzangles<<
2766 //____________Calcul Coef Mean_________________________________________________
2768 //_____________________________________________________________________________
2769 Bool_t AliTRDCalibraFit::CalculT0CoefMean()
2772 // For the detector Dect calcul the mean time 0
2773 // for the calibration group idect from the choosen database
2776 fCurrentCoef2[1] = 0.0;
2777 if(fDebugLevel != 1){
2778 if((fCalibraMode->GetNz(1) > 0) ||
2779 (fCalibraMode->GetNrphi(1) > 0)) {
2780 for (Int_t row = fCalibraMode->GetRowMin(1); row < fCalibraMode->GetRowMax(1); row++) {
2781 for (Int_t col = fCalibraMode->GetColMin(1); col < fCalibraMode->GetColMax(1); col++) {
2782 fCurrentCoef2[1] += (Float_t) (fCalROC2->GetValue(col,row)+fCalDet2->GetValue(fCountDet));
2785 fCurrentCoef2[1] = fCurrentCoef2[1] / ((fCalibraMode->GetColMax(1)-fCalibraMode->GetColMin(1))*(fCalibraMode->GetRowMax(1)-fCalibraMode->GetRowMin(1)));
2789 fCurrentCoef2[1] = fCalDet2->GetValue(fCountDet);
2792 for(Int_t row = 0; row < fGeo->GetRowMax(GetLayer(fCountDet),GetStack(fCountDet),GetSector(fCountDet)); row++){
2793 for(Int_t col = 0; col < fGeo->GetColMax(GetLayer(fCountDet)); col++){
2794 fCurrentCoef2[1] += (Float_t) (fCalROC2->GetValue(col,row)+fCalDet2->GetValue(fCountDet));
2797 fCurrentCoef2[1] = fCurrentCoef2[1] / ((fGeo->GetRowMax(GetLayer(fCountDet),GetStack(fCountDet),GetSector(fCountDet)))*(fGeo->GetColMax(GetLayer(fCountDet))));
2804 //_____________________________________________________________________________
2805 Bool_t AliTRDCalibraFit::CalculChargeCoefMean(Bool_t vrai)
2808 // For the detector Dect calcul the mean gain factor
2809 // for the calibration group idect from the choosen database
2812 fCurrentCoef[1] = 0.0;
2813 if(fDebugLevel != 1){
2814 if ((fCalibraMode->GetNz(0) > 0) ||
2815 (fCalibraMode->GetNrphi(0) > 0)) {
2816 for (Int_t row = fCalibraMode->GetRowMin(0); row < fCalibraMode->GetRowMax(0); row++) {
2817 for (Int_t col = fCalibraMode->GetColMin(0); col < fCalibraMode->GetColMax(0); col++) {
2818 fCurrentCoef[1] += (Float_t) (fCalROC->GetValue(col,row)*fCalDet->GetValue(fCountDet));
2819 if (vrai) fScaleFitFactor += (Float_t) (fCalROC->GetValue(col,row)*fCalDet->GetValue(fCountDet));
2822 fCurrentCoef[1] = fCurrentCoef[1] / ((fCalibraMode->GetColMax(0)-fCalibraMode->GetColMin(0))*(fCalibraMode->GetRowMax(0)-fCalibraMode->GetRowMin(0)));
2826 fCurrentCoef[1] = (Float_t) fCalDet->GetValue(fCountDet);
2827 if (vrai) fScaleFitFactor += ((Float_t) fCalDet->GetValue(fCountDet))*(fCalibraMode->GetColMax(0)-fCalibraMode->GetColMin(0))*(fCalibraMode->GetRowMax(0)-fCalibraMode->GetRowMin(0));
2832 //_____________________________________________________________________________
2833 Bool_t AliTRDCalibraFit::CalculPRFCoefMean()
2836 // For the detector Dect calcul the mean sigma of pad response
2837 // function for the calibration group idect from the choosen database
2840 fCurrentCoef[1] = 0.0;
2841 if(fDebugLevel != 1){
2842 for (Int_t row = fCalibraMode->GetRowMin(2); row < fCalibraMode->GetRowMax(2); row++) {
2843 for (Int_t col = fCalibraMode->GetColMin(2); col < fCalibraMode->GetColMax(2); col++) {
2844 fCurrentCoef[1] += (Float_t) fCalROC->GetValue(col,row);
2847 fCurrentCoef[1] = fCurrentCoef[1] / ((fCalibraMode->GetColMax(2)-fCalibraMode->GetColMin(2))*(fCalibraMode->GetRowMax(2)-fCalibraMode->GetRowMin(2)));
2851 //_____________________________________________________________________________
2852 Bool_t AliTRDCalibraFit::CalculVdriftCoefMean()
2855 // For the detector dect calcul the mean drift velocity for the
2856 // calibration group idect from the choosen database
2859 fCurrentCoef[1] = 0.0;
2860 if(fDebugLevel != 1){
2861 if ((fCalibraMode->GetNz(1) > 0) ||
2862 (fCalibraMode->GetNrphi(1) > 0)) {
2863 for (Int_t row = fCalibraMode->GetRowMin(1); row < fCalibraMode->GetRowMax(1); row++) {
2864 for (Int_t col = fCalibraMode->GetColMin(1); col < fCalibraMode->GetColMax(1); col++) {
2865 fCurrentCoef[1] += (Float_t) (fCalROC->GetValue(col,row)*fCalDet->GetValue(fCountDet));
2868 fCurrentCoef[1] = fCurrentCoef[1] / ((fCalibraMode->GetColMax(1)-fCalibraMode->GetColMin(1))*(fCalibraMode->GetRowMax(1)-fCalibraMode->GetRowMin(1)));
2872 fCurrentCoef[1] = (Float_t) fCalDet->GetValue(fCountDet);
2877 //_____________________________________________________________________________
2878 Bool_t AliTRDCalibraFit::CalculVdriftLorentzCoef()
2881 // For the detector fCountDet, mean drift velocity and tan lorentzangle
2884 fCurrentCoef[1] = fCalDet->GetValue(fCountDet);
2885 fCurrentCoef2[1] = fCalDet2->GetValue(fCountDet);
2889 //_____________________________________________________________________________
2890 Float_t AliTRDCalibraFit::GetPRFDefault(Int_t layer) const
2893 // Default width of the PRF if there is no database as reference
2898 //case 0: return 0.515;
2899 //case 1: return 0.502;
2900 //case 2: return 0.491;
2901 //case 3: return 0.481;
2902 //case 4: return 0.471;
2903 //case 5: return 0.463;
2904 //default: return 0.0;
2907 case 0: return 0.538429;
2908 case 1: return 0.524302;
2909 case 2: return 0.511591;
2910 case 3: return 0.500140;
2911 case 4: return 0.489821;
2912 case 5: return 0.480524;
2913 default: return 0.0;
2916 //________________________________________________________________________________
2917 void AliTRDCalibraFit::SetCalROC(Int_t i)
2920 // Set the calib object for fCountDet
2923 Float_t value = 0.0;
2925 //Get the CalDet object
2927 AliTRDcalibDB *cal = AliTRDcalibDB::Instance();
2929 AliInfo("Could not get calibDB");
2935 if(fCalROC) delete fCalROC;
2936 fCalROC = new AliTRDCalROC(*(cal->GetGainFactorROC(fCountDet)));
2939 if(fCalROC) delete fCalROC;
2940 if(fCalROC2) delete fCalROC2;
2941 fCalROC = new AliTRDCalROC(*(cal->GetVdriftROC(fCountDet)));
2942 fCalROC2 = new AliTRDCalROC(*(cal->GetT0ROC(fCountDet)));
2945 if(fCalROC) delete fCalROC;
2946 fCalROC = new AliTRDCalROC(*(cal->GetPRFROC(fCountDet))); break;
2954 if(fCalROC) delete fCalROC;
2955 fCalROC = new AliTRDCalROC(GetLayer(fCountDet),GetStack(fCountDet));
2956 for(Int_t k = 0; k < fCalROC->GetNchannels(); k++){
2957 fCalROC->SetValue(k,1.0);
2961 if(fCalROC) delete fCalROC;
2962 if(fCalROC2) delete fCalROC2;
2963 fCalROC = new AliTRDCalROC(GetLayer(fCountDet),GetStack(fCountDet));
2964 fCalROC2 = new AliTRDCalROC(GetLayer(fCountDet),GetStack(fCountDet));
2965 for(Int_t k = 0; k < fCalROC->GetNchannels(); k++){
2966 fCalROC->SetValue(k,1.0);
2967 fCalROC2->SetValue(k,0.0);
2971 if(fCalROC) delete fCalROC;
2972 value = GetPRFDefault(GetLayer(fCountDet));
2973 fCalROC = new AliTRDCalROC(GetLayer(fCountDet),GetStack(fCountDet));
2974 for(Int_t k = 0; k < fCalROC->GetNchannels(); k++){
2975 fCalROC->SetValue(k,value);
2983 //____________Fit Methods______________________________________________________
2985 //_____________________________________________________________________________
2986 void AliTRDCalibraFit::FitPente(TH1* projPH)
2989 // Slope methode for the drift velocity
2993 const Float_t kDrWidth = AliTRDgeometry::DrThick();
3000 fCurrentCoefE = 0.0;
3001 fCurrentCoefE2 = 0.0;
3002 fCurrentCoef[0] = 0.0;
3003 fCurrentCoef2[0] = 0.0;
3004 TLine *line = new TLine();
3007 TAxis *xpph = projPH->GetXaxis();
3008 Int_t nbins = xpph->GetNbins();
3009 Double_t lowedge = xpph->GetBinLowEdge(1);
3010 Double_t upedge = xpph->GetBinUpEdge(xpph->GetNbins());
3011 Double_t widbins = (upedge - lowedge) / nbins;
3012 Double_t limit = upedge + 0.5 * widbins;
3015 // Beginning of the signal
3016 TH1D *pentea = new TH1D("pentea","pentea",projPH->GetNbinsX(),0,(Float_t) limit);
3017 for (Int_t k = 1; k < projPH->GetNbinsX(); k++) {
3018 pentea->SetBinContent(k,(Double_t) (projPH->GetBinContent(k+1) - projPH->GetBinContent(k)));
3020 binmax = (Int_t) pentea->GetMaximumBin();
3023 AliInfo("Put the binmax from 1 to 2 to enable the fit");
3025 if (binmax >= nbins) {
3028 AliInfo("Put the binmax from nbins-1 to nbins-2 to enable the fit");
3030 pentea->Fit("pol2","0MR","",TMath::Max(pentea->GetBinCenter(binmax-1),0.0),pentea->GetBinCenter(binmax+1));
3031 Float_t l3P1am = pentea->GetFunction("pol2")->GetParameter(1);
3032 Float_t l3P2am = pentea->GetFunction("pol2")->GetParameter(2);
3033 Float_t l3P1amE = pentea->GetFunction("pol2")->GetParError(1);
3034 Float_t l3P2amE = pentea->GetFunction("pol2")->GetParError(2);
3036 fPhd[0] = -(l3P1am / (2 * l3P2am));
3039 if((l3P1am != 0.0) && (l3P2am != 0.0)){
3040 fCurrentCoefE2 = (l3P1amE/l3P1am + l3P2amE/l3P2am)*fPhd[0];
3043 // Amplification region
3046 for (Int_t kbin = 1; kbin < projPH->GetNbinsX(); kbin ++) {
3047 if (((projPH->GetBinContent(kbin+1) - projPH->GetBinContent(kbin)) <= 0.0) && (ju == 0) && (kbin > (fPhd[0]/widbins))) {
3054 AliInfo("Put the binmax from 1 to 2 to enable the fit");
3056 if (binmax >= nbins) {
3059 AliInfo("Put the binmax from nbins-1 to nbins-2 to enable the fit");
3061 projPH->Fit("pol2","0MR","",TMath::Max(projPH->GetBinCenter(binmax-1),0.0),projPH->GetBinCenter(binmax+1));
3062 Float_t l3P1amf = projPH->GetFunction("pol2")->GetParameter(1);
3063 Float_t l3P2amf = projPH->GetFunction("pol2")->GetParameter(2);
3064 Float_t l3P1amfE = projPH->GetFunction("pol2")->GetParError(1);
3065 Float_t l3P2amfE = projPH->GetFunction("pol2")->GetParError(2);
3067 fPhd[1] = -(l3P1amf / (2 * l3P2amf));
3069 if((l3P1amf != 0.0) && (l3P2amf != 0.0)){
3070 fCurrentCoefE = (l3P1amfE/l3P1amf + l3P2amfE/l3P2amf)*fPhd[1];
3073 fCurrentCoefE2 = fCurrentCoefE;
3076 TH1D *pente = new TH1D("pente","pente",projPH->GetNbinsX(),0,(Float_t) limit);
3077 for (Int_t k = TMath::Min(binmax+4,projPH->GetNbinsX()); k < projPH->GetNbinsX(); k++) {
3078 pente->SetBinContent(k,(Double_t) (projPH->GetBinContent(k+1) - projPH->GetBinContent(k)));
3081 if(pente->GetEntries() > 0) binmin = (Int_t) pente->GetMinimumBin();
3084 AliInfo("Put the binmax from 1 to 2 to enable the fit");
3086 if (binmin >= nbins) {
3089 AliInfo("Put the binmax from nbins-1 to nbins-2 to enable the fit");
3094 ,TMath::Max(pente->GetBinCenter(binmin-1), 0.0)
3095 ,TMath::Min(pente->GetBinCenter(binmin+1),(Double_t) limit));
3096 Float_t l3P1dr = pente->GetFunction("pol2")->GetParameter(1);
3097 Float_t l3P2dr = pente->GetFunction("pol2")->GetParameter(2);
3098 Float_t l3P1drE = pente->GetFunction("pol2")->GetParError(1);
3099 Float_t l3P2drE = pente->GetFunction("pol2")->GetParError(2);
3101 fPhd[2] = -(l3P1dr / (2 * l3P2dr));
3103 if((l3P1dr != 0.0) && (l3P2dr != 0.0)){
3104 fCurrentCoefE += (l3P1drE/l3P1dr + l3P2drE/l3P2dr)*fPhd[2];
3106 Float_t fPhdt0 = 0.0;
3107 Float_t t0Shift = 0.0;
3110 t0Shift = fT0Shift1;
3114 t0Shift = fT0Shift0;
3117 if ((fPhd[2] > fPhd[0]) &&
3118 (fPhd[2] > fPhd[1]) &&
3119 (fPhd[1] > fPhd[0]) &&
3121 fCurrentCoef[0] = (kDrWidth) / (fPhd[2]-fPhd[1]);
3122 fNumberFitSuccess++;
3124 if (fPhdt0 >= 0.0) {
3125 fCurrentCoef2[0] = (fPhdt0 - t0Shift) / widbins;
3126 if (fCurrentCoef2[0] < -1.0) {
3127 fCurrentCoef2[0] = fCurrentCoef2[1];
3131 fCurrentCoef2[0] = fCurrentCoef2[1];
3136 fCurrentCoef[0] = -TMath::Abs(fCurrentCoef[1]);
3137 fCurrentCoef2[0] = fCurrentCoef2[1];
3140 if (fDebugLevel == 1) {
3141 TCanvas *cpentei = new TCanvas("cpentei","cpentei",50,50,600,800);
3144 line->SetLineColor(2);
3145 line->DrawLine(fPhd[0],0,fPhd[0],projPH->GetMaximum());
3146 line->DrawLine(fPhd[1],0,fPhd[1],projPH->GetMaximum());
3147 line->DrawLine(fPhd[2],0,fPhd[2],projPH->GetMaximum());
3148 AliInfo(Form("fPhd[0] (beginning of the signal): %f" ,(Float_t) fPhd[0]));
3149 AliInfo(Form("fPhd[1] (end of the amplification region): %f" ,(Float_t) fPhd[1]));
3150 AliInfo(Form("fPhd[2] (end of the drift region): %f" ,(Float_t) fPhd[2]));
3151 AliInfo(Form("fVriftCoef[1] (with only the drift region(default)): %f",(Float_t) fCurrentCoef[0]));
3152 TCanvas *cpentei2 = new TCanvas("cpentei2","cpentei2",50,50,600,800);
3155 TCanvas *cpentei3 = new TCanvas("cpentei3","cpentei3",50,50,600,800);
3164 //_____________________________________________________________________________
3165 void AliTRDCalibraFit::FitLagrangePoly(TH1* projPH)
3168 // Slope methode but with polynomes de Lagrange
3172 const Float_t kDrWidth = AliTRDgeometry::DrThick();
3175 Double_t *x = new Double_t[5];
3176 Double_t *y = new Double_t[5];
3191 fCurrentCoefE = 0.0;
3192 fCurrentCoefE2 = 1.0;
3193 fCurrentCoef[0] = 0.0;
3194 fCurrentCoef2[0] = 0.0;
3195 TLine *line = new TLine();
3196 TF1 * polynome = 0x0;
3197 TF1 * polynomea = 0x0;
3198 TF1 * polynomeb = 0x0;
3202 TAxis *xpph = projPH->GetXaxis();
3203 Int_t nbins = xpph->GetNbins();
3204 Double_t lowedge = xpph->GetBinLowEdge(1);
3205 Double_t upedge = xpph->GetBinUpEdge(xpph->GetNbins());
3206 Double_t widbins = (upedge - lowedge) / nbins;
3207 Double_t limit = upedge + 0.5 * widbins;
3212 // Beginning of the signal
3213 TH1D *pentea = new TH1D("pentea","pentea",projPH->GetNbinsX(),0,(Float_t) limit);
3214 for (Int_t k = 1; k < projPH->GetNbinsX(); k++) {
3215 pentea->SetBinContent(k,(Double_t) (projPH->GetBinContent(k+1) - projPH->GetBinContent(k)));
3218 binmax = (Int_t) pentea->GetMaximumBin();
3220 Double_t minnn = 0.0;
3221 Double_t maxxx = 0.0;
3223 Int_t kase = nbins-binmax;
3231 minnn = pentea->GetBinCenter(binmax-2);
3232 maxxx = pentea->GetBinCenter(binmax);
3233 x[0] = pentea->GetBinCenter(binmax-2);
3234 x[1] = pentea->GetBinCenter(binmax-1);
3235 x[2] = pentea->GetBinCenter(binmax);
3236 y[0] = pentea->GetBinContent(binmax-2);
3237 y[1] = pentea->GetBinContent(binmax-1);
3238 y[2] = pentea->GetBinContent(binmax);
3239 c = CalculPolynomeLagrange2(x,y);
3240 AliInfo("At the limit for beginning!");
3243 minnn = pentea->GetBinCenter(binmax-2);
3244 maxxx = pentea->GetBinCenter(binmax+1);
3245 x[0] = pentea->GetBinCenter(binmax-2);
3246 x[1] = pentea->GetBinCenter(binmax-1);
3247 x[2] = pentea->GetBinCenter(binmax);
3248 x[3] = pentea->GetBinCenter(binmax+1);
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 c = CalculPolynomeLagrange3(x,y);
3261 minnn = pentea->GetBinCenter(binmax);
3262 maxxx = pentea->GetBinCenter(binmax+2);
3263 x[0] = pentea->GetBinCenter(binmax);
3264 x[1] = pentea->GetBinCenter(binmax+1);
3265 x[2] = pentea->GetBinCenter(binmax+2);
3266 y[0] = pentea->GetBinContent(binmax);
3267 y[1] = pentea->GetBinContent(binmax+1);
3268 y[2] = pentea->GetBinContent(binmax+2);
3269 c = CalculPolynomeLagrange2(x,y);
3272 minnn = pentea->GetBinCenter(binmax-1);
3273 maxxx = pentea->GetBinCenter(binmax+2);
3274 x[0] = pentea->GetBinCenter(binmax-1);
3275 x[1] = pentea->GetBinCenter(binmax);
3276 x[2] = pentea->GetBinCenter(binmax+1);
3277 x[3] = pentea->GetBinCenter(binmax+2);
3278 y[0] = pentea->GetBinContent(binmax-1);
3279 y[1] = pentea->GetBinContent(binmax);
3280 y[2] = pentea->GetBinContent(binmax+1);
3281 y[3] = pentea->GetBinContent(binmax+2);
3282 c = CalculPolynomeLagrange3(x,y);
3285 minnn = pentea->GetBinCenter(binmax-2);
3286 maxxx = pentea->GetBinCenter(binmax+2);
3287 x[0] = pentea->GetBinCenter(binmax-2);
3288 x[1] = pentea->GetBinCenter(binmax-1);
3289 x[2] = pentea->GetBinCenter(binmax);
3290 x[3] = pentea->GetBinCenter(binmax+1);
3291 x[4] = pentea->GetBinCenter(binmax+2);
3292 y[0] = pentea->GetBinContent(binmax-2);
3293 y[1] = pentea->GetBinContent(binmax-1);
3294 y[2] = pentea->GetBinContent(binmax);
3295 y[3] = pentea->GetBinContent(binmax+1);
3296 y[4] = pentea->GetBinContent(binmax+2);
3297 c = CalculPolynomeLagrange4(x,y);
3305 polynomeb = new TF1("polb","[0]+[1]*x+[2]*x*x+[3]*x*x*x+[4]*x*x*x*x",minnn,maxxx);
3306 polynomeb->SetParameters(c[0],c[1],c[2],c[3],c[4]);
3308 Double_t step = (maxxx-minnn)/10000;
3310 Double_t maxvalue = 0.0;
3311 Double_t placemaximum = minnn;
3312 for(Int_t o = 0; o < 10000; o++){
3313 if(o == 0) maxvalue = polynomeb->Eval(l);
3314 if(maxvalue < (polynomeb->Eval(l))){
3315 maxvalue = polynomeb->Eval(l);
3320 fPhd[0] = placemaximum;
3323 // Amplification region
3326 for (Int_t kbin = 1; kbin < projPH->GetNbinsX(); kbin ++) {
3327 if (((projPH->GetBinContent(kbin+1) - projPH->GetBinContent(kbin)) <= 0.0) && (ju == 0) && (kbin > (fPhd[0]/widbins))) {
3333 Double_t minn = 0.0;
3334 Double_t maxx = 0.0;
3346 Int_t kase1 = nbins - binmax;
3348 //Determination of minn and maxx
3349 //case binmax = nbins
3354 minn = projPH->GetBinCenter(binmax-2);
3355 maxx = projPH->GetBinCenter(binmax);
3356 x[0] = projPH->GetBinCenter(binmax-2);
3357 x[1] = projPH->GetBinCenter(binmax-1);
3358 x[2] = projPH->GetBinCenter(binmax);
3359 y[0] = projPH->GetBinContent(binmax-2);
3360 y[1] = projPH->GetBinContent(binmax-1);
3361 y[2] = projPH->GetBinContent(binmax);
3362 c = CalculPolynomeLagrange2(x,y);
3363 //AliInfo("At the limit for the drift!");
3366 minn = projPH->GetBinCenter(binmax-2);
3367 maxx = projPH->GetBinCenter(binmax+1);
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 y[0] = projPH->GetBinContent(binmax-2);
3373 y[1] = projPH->GetBinContent(binmax-1);
3374 y[2] = projPH->GetBinContent(binmax);
3375 y[3] = projPH->GetBinContent(binmax+1);
3376 c = CalculPolynomeLagrange3(x,y);
3385 minn = projPH->GetBinCenter(binmax);
3386 maxx = projPH->GetBinCenter(binmax+2);
3387 x[0] = projPH->GetBinCenter(binmax);
3388 x[1] = projPH->GetBinCenter(binmax+1);
3389 x[2] = projPH->GetBinCenter(binmax+2);
3390 y[0] = projPH->GetBinContent(binmax);
3391 y[1] = projPH->GetBinContent(binmax+1);
3392 y[2] = projPH->GetBinContent(binmax+2);
3393 c = CalculPolynomeLagrange2(x,y);
3396 minn = projPH->GetBinCenter(binmax-1);
3397 maxx = projPH->GetBinCenter(binmax+2);
3398 x[0] = projPH->GetBinCenter(binmax-1);
3399 x[1] = projPH->GetBinCenter(binmax);
3400 x[2] = projPH->GetBinCenter(binmax+1);
3401 x[3] = projPH->GetBinCenter(binmax+2);
3402 y[0] = projPH->GetBinContent(binmax-1);
3403 y[1] = projPH->GetBinContent(binmax);
3404 y[2] = projPH->GetBinContent(binmax+1);
3405 y[3] = projPH->GetBinContent(binmax+2);
3406 c = CalculPolynomeLagrange3(x,y);
3409 minn = projPH->GetBinCenter(binmax-2);
3410 maxx = projPH->GetBinCenter(binmax+2);
3411 x[0] = projPH->GetBinCenter(binmax-2);
3412 x[1] = projPH->GetBinCenter(binmax-1);
3413 x[2] = projPH->GetBinCenter(binmax);
3414 x[3] = projPH->GetBinCenter(binmax+1);
3415 x[4] = projPH->GetBinCenter(binmax+2);
3416 y[0] = projPH->GetBinContent(binmax-2);
3417 y[1] = projPH->GetBinContent(binmax-1);
3418 y[2] = projPH->GetBinContent(binmax);
3419 y[3] = projPH->GetBinContent(binmax+1);
3420 y[4] = projPH->GetBinContent(binmax+2);
3421 c = CalculPolynomeLagrange4(x,y);
3428 polynomea = new TF1("pola","[0]+[1]*x+[2]*x*x+[3]*x*x*x+[4]*x*x*x*x",minn,maxx);
3429 polynomea->SetParameters(c[0],c[1],c[2],c[3],c[4]);
3431 Double_t step = (maxx-minn)/1000;
3433 Double_t maxvalue = 0.0;
3434 Double_t placemaximum = minn;
3435 for(Int_t o = 0; o < 1000; o++){
3436 if(o == 0) maxvalue = polynomea->Eval(l);
3437 if(maxvalue < (polynomea->Eval(l))){
3438 maxvalue = polynomea->Eval(l);
3443 fPhd[1] = placemaximum;
3447 TH1D *pente = new TH1D("pente","pente", projPH->GetNbinsX(),0,(Float_t) limit);
3448 for (Int_t k = TMath::Min(binmax+4, projPH->GetNbinsX()); k < projPH->GetNbinsX(); k++) {
3449 pente->SetBinContent(k,(Double_t) (projPH->GetBinContent(k+1) - projPH->GetBinContent(k)));
3452 if(pente->GetEntries() > 0) binmin = (Int_t) pente->GetMinimumBin();
3458 AliInfo("Put the binmax from 1 to 2 to enable the fit");
3462 if((projPH->GetBinContent(binmin)-projPH->GetBinError(binmin)) < (projPH->GetBinContent(binmin+1))) put = kFALSE;
3463 if((projPH->GetBinContent(binmin)+projPH->GetBinError(binmin)) > (projPH->GetBinContent(binmin-1))) put = kFALSE;
3477 Bool_t case1 = kFALSE;
3478 Bool_t case2 = kFALSE;
3479 Bool_t case4 = kFALSE;
3481 //Determination of min and max
3482 //case binmin <= nbins-3
3484 if((binmin <= (nbins-3)) && ((binmin-2) >= TMath::Min(binmax+4, projPH->GetNbinsX()))){
3485 min = pente->GetBinCenter(binmin-2);
3486 max = pente->GetBinCenter(binmin+2);
3487 x[0] = pente->GetBinCenter(binmin-2);
3488 x[1] = pente->GetBinCenter(binmin-1);
3489 x[2] = pente->GetBinCenter(binmin);
3490 x[3] = pente->GetBinCenter(binmin+1);
3491 x[4] = pente->GetBinCenter(binmin+2);
3492 y[0] = pente->GetBinContent(binmin-2);
3493 y[1] = pente->GetBinContent(binmin-1);
3494 y[2] = pente->GetBinContent(binmin);
3495 y[3] = pente->GetBinContent(binmin+1);
3496 y[4] = pente->GetBinContent(binmin+2);
3497 //Calcul the polynome de Lagrange
3498 c = CalculPolynomeLagrange4(x,y);
3500 if((pente->GetBinContent(binmin+2) <= pente->GetBinContent(binmin+1)) &&
3501 (pente->GetBinContent(binmin-2) <= pente->GetBinContent(binmin-1))) put = kFALSE;
3502 if(((binmin+3) <= (nbins-1)) &&
3503 (pente->GetBinContent(binmin+3) <= pente->GetBinContent(binmin+2)) &&
3504 ((binmin-3) >= TMath::Min(binmax+4, projPH->GetNbinsX())) &&
3505 (pente->GetBinContent(binmin-3) <= pente->GetBinContent(binmin-2))) put = kFALSE;
3506 if((pente->GetBinContent(binmin+2) <= pente->GetBinContent(binmin+1)) &&
3507 (pente->GetBinContent(binmin-2) > pente->GetBinContent(binmin-1))) case1 = kTRUE;
3508 if((pente->GetBinContent(binmin+2) > pente->GetBinContent(binmin+1)) &&
3509 (pente->GetBinContent(binmin-2) <= pente->GetBinContent(binmin-1))) case4 = kTRUE;
3511 //case binmin = nbins-2
3513 if(((binmin == (nbins-2)) && ((binmin-2) >= TMath::Min(binmax+4, projPH->GetNbinsX()))) ||
3515 min = pente->GetBinCenter(binmin-2);
3516 max = pente->GetBinCenter(binmin+1);
3517 x[0] = pente->GetBinCenter(binmin-2);
3518 x[1] = pente->GetBinCenter(binmin-1);
3519 x[2] = pente->GetBinCenter(binmin);
3520 x[3] = pente->GetBinCenter(binmin+1);
3521 y[0] = pente->GetBinContent(binmin-2);
3522 y[1] = pente->GetBinContent(binmin-1);
3523 y[2] = pente->GetBinContent(binmin);
3524 y[3] = pente->GetBinContent(binmin+1);
3525 //Calcul the polynome de Lagrange
3526 c = CalculPolynomeLagrange3(x,y);
3527 //richtung +: nothing
3529 if((pente->GetBinContent(binmin-2) <= pente->GetBinContent(binmin-1))) case2 = kTRUE;
3532 if(((binmin <= (nbins-3)) && ((binmin-1) == TMath::Min(binmax+4, projPH->GetNbinsX()))) ||
3534 min = pente->GetBinCenter(binmin-1);
3535 max = pente->GetBinCenter(binmin+2);
3536 x[0] = pente->GetBinCenter(binmin-1);
3537 x[1] = pente->GetBinCenter(binmin);
3538 x[2] = pente->GetBinCenter(binmin+1);
3539 x[3] = pente->GetBinCenter(binmin+2);
3540 y[0] = pente->GetBinContent(binmin-1);
3541 y[1] = pente->GetBinContent(binmin);
3542 y[2] = pente->GetBinContent(binmin+1);
3543 y[3] = pente->GetBinContent(binmin+2);
3544 //Calcul the polynome de Lagrange
3545 c = CalculPolynomeLagrange3(x,y);
3547 if((pente->GetBinContent(binmin+2) <= pente->GetBinContent(binmin+1))) case2 = kTRUE;
3550 if((binmin <= (nbins-3)) && (binmin == TMath::Min(binmax+4, projPH->GetNbinsX()))){
3551 min = pente->GetBinCenter(binmin);
3552 max = pente->GetBinCenter(binmin+2);
3553 x[0] = pente->GetBinCenter(binmin);
3554 x[1] = pente->GetBinCenter(binmin+1);
3555 x[2] = pente->GetBinCenter(binmin+2);
3556 y[0] = pente->GetBinContent(binmin);
3557 y[1] = pente->GetBinContent(binmin+1);
3558 y[2] = pente->GetBinContent(binmin+2);
3559 //Calcul the polynome de Lagrange
3560 c = CalculPolynomeLagrange2(x,y);
3562 if((pente->GetBinContent(binmin+2) <= pente->GetBinContent(binmin+1))) put = kFALSE;
3565 if(((binmin == (nbins-2)) && ((binmin-1) == TMath::Min(binmax+4, projPH->GetNbinsX()))) ||
3567 min = pente->GetBinCenter(binmin-1);
3568 max = pente->GetBinCenter(binmin+1);
3569 x[0] = pente->GetBinCenter(binmin-1);
3570 x[1] = pente->GetBinCenter(binmin);
3571 x[2] = pente->GetBinCenter(binmin+1);
3572 y[0] = pente->GetBinContent(binmin-1);
3573 y[1] = pente->GetBinContent(binmin);
3574 y[2] = pente->GetBinContent(binmin+1);
3575 //Calcul the polynome de Lagrange
3576 c = CalculPolynomeLagrange2(x,y);
3577 //richtung +: nothing
3578 //richtung -: nothing
3580 //case binmin = nbins-1
3582 if((binmin == (nbins-1)) && ((binmin-2) >= TMath::Min(binmax+4, projPH->GetNbinsX()))){
3583 min = pente->GetBinCenter(binmin-2);
3584 max = pente->GetBinCenter(binmin);
3585 x[0] = pente->GetBinCenter(binmin-2);
3586 x[1] = pente->GetBinCenter(binmin-1);
3587 x[2] = pente->GetBinCenter(binmin);
3588 y[0] = pente->GetBinContent(binmin-2);
3589 y[1] = pente->GetBinContent(binmin-1);
3590 y[2] = pente->GetBinContent(binmin);
3591 //Calcul the polynome de Lagrange
3592 c = CalculPolynomeLagrange2(x,y);
3593 //AliInfo("At the limit for the drift!");
3594 //fluctuation too big!
3595 //richtung +: nothing
3597 if((pente->GetBinContent(binmin-2) <= pente->GetBinContent(binmin-1))) put = kFALSE;
3599 if((binmin == (nbins-1)) && ((binmin-2) < TMath::Min(binmax+4, projPH->GetNbinsX()))) {
3601 AliInfo("At the limit for the drift and not usable!");
3605 if((binmin == (nbins-2)) && ((binmin-1) < TMath::Min(binmax+4, projPH->GetNbinsX()))){
3607 AliInfo("For the drift...problem!");
3609 //pass but should not happen
3610 if((binmin <= (nbins-3)) && (binmin < TMath::Min(binmax+4, projPH->GetNbinsX()))){
3612 AliInfo("For the drift...problem!");
3616 polynome = new TF1("pol","[0]+[1]*x+[2]*x*x+[3]*x*x*x+[4]*x*x*x*x",min,max);
3617 polynome->SetParameters(c[0],c[1],c[2],c[3],c[4]);
3618 //AliInfo(Form("GetMinimum of the function %f",polynome->GetMinimumX()));
3619 Double_t step = (max-min)/1000;
3621 Double_t minvalue = 0.0;
3622 Double_t placeminimum = min;
3623 for(Int_t o = 0; o < 1000; o++){
3624 if(o == 0) minvalue = polynome->Eval(l);
3625 if(minvalue > (polynome->Eval(l))){
3626 minvalue = polynome->Eval(l);
3631 fPhd[2] = placeminimum;
3634 Float_t fPhdt0 = 0.0;
3635 Float_t t0Shift = 0.0;
3638 t0Shift = fT0Shift1;
3642 t0Shift = fT0Shift0;
3645 if ((fPhd[2] > fPhd[0]) &&
3646 (fPhd[2] > fPhd[1]) &&
3647 (fPhd[1] > fPhd[0]) &&
3649 fCurrentCoef[0] = (kDrWidth) / (fPhd[2]-fPhd[1]);
3650 fNumberFitSuccess++;
3651 if (fPhdt0 >= 0.0) {
3652 fCurrentCoef2[0] = (fPhdt0 - t0Shift) / widbins;
3653 if (fCurrentCoef2[0] < -1.0) {
3654 fCurrentCoef2[0] = fCurrentCoef2[1];
3658 fCurrentCoef2[0] = fCurrentCoef2[1];
3662 fCurrentCoef[0] = -TMath::Abs(fCurrentCoef[1]);
3663 fCurrentCoef2[0] = fCurrentCoef2[1];
3664 //printf("Fit failed!\n");
3667 if (fDebugLevel == 1) {
3668 TCanvas *cpentei = new TCanvas("cpentei","cpentei",50,50,600,800);
3671 line->SetLineColor(2);
3672 line->DrawLine(fPhd[0],0,fPhd[0],projPH->GetMaximum());
3673 line->DrawLine(fPhd[1],0,fPhd[1],projPH->GetMaximum());
3674 line->DrawLine(fPhd[2],0,fPhd[2],projPH->GetMaximum());
3675 AliInfo(Form("fPhd[0] (beginning of the signal): %f" ,(Float_t) fPhd[0]));
3676 AliInfo(Form("fPhd[1] (end of the amplification region): %f" ,(Float_t) fPhd[1]));
3677 AliInfo(Form("fPhd[2] (end of the drift region): %f" ,(Float_t) fPhd[2]));
3678 AliInfo(Form("fVriftCoef[3] (with only the drift region(default)): %f",(Float_t) fCurrentCoef[0]));
3679 TCanvas *cpentei2 = new TCanvas("cpentei2","cpentei2",50,50,600,800);
3682 TCanvas *cpentei3 = new TCanvas("cpentei3","cpentei3",50,50,600,800);
3687 if(pentea) delete pentea;
3688 if(pente) delete pente;
3689 if(polynome) delete polynome;
3690 if(polynomea) delete polynomea;
3691 if(polynomeb) delete polynomeb;
3695 if(line) delete line;
3699 projPH->SetDirectory(0);
3703 //_____________________________________________________________________________
3704 void AliTRDCalibraFit::FitPH(TH1* projPH, Int_t idect)
3707 // Fit methode for the drift velocity
3711 const Float_t kDrWidth = AliTRDgeometry::DrThick();
3714 TAxis *xpph = projPH->GetXaxis();
3715 Double_t upedge = xpph->GetBinUpEdge(xpph->GetNbins());
3717 TF1 *fPH = new TF1("fPH",AliTRDCalibraFit::PH,-0.05,3.2,6);
3718 fPH->SetParameter(0,0.469); // Scaling
3719 fPH->SetParameter(1,0.18); // Start
3720 fPH->SetParameter(2,0.0857325); // AR
3721 fPH->SetParameter(3,1.89); // DR
3722 fPH->SetParameter(4,0.08); // QA/QD
3723 fPH->SetParameter(5,0.0); // Baseline
3725 TLine *line = new TLine();
3727 fCurrentCoef[0] = 0.0;
3728 fCurrentCoef2[0] = 0.0;
3729 fCurrentCoefE = 0.0;
3730 fCurrentCoefE2 = 0.0;
3732 if (idect%fFitPHPeriode == 0) {
3734 AliInfo(Form("The detector %d will be fitted",idect));
3735 fPH->SetParameter(0,(projPH->Integral()-(projPH->GetBinContent(1)*projPH->GetNbinsX())) * 0.00028); // Scaling
3736 fPH->SetParameter(1,fPhd[0] - 0.1); // Start
3737 fPH->SetParameter(2,fPhd[1] - fPhd[0]); // AR
3738 fPH->SetParameter(3,fPhd[2] - fPhd[1]); // DR
3739 fPH->SetParameter(4,0.225); // QA/QD
3740 fPH->SetParameter(5,(Float_t) projPH->GetBinContent(1));
3742 if (fDebugLevel != 1) {
3743 projPH->Fit(fPH,"0M","",0.0,upedge);
3746 TCanvas *cpente = new TCanvas("cpente","cpente",50,50,600,800);
3748 projPH->Fit(fPH,"M+","",0.0,upedge);
3750 line->SetLineColor(4);
3751 line->DrawLine(fPH->GetParameter(1)
3753 ,fPH->GetParameter(1)
3754 ,projPH->GetMaximum());
3755 line->DrawLine(fPH->GetParameter(1)+fPH->GetParameter(2)
3757 ,fPH->GetParameter(1)+fPH->GetParameter(2)
3758 ,projPH->GetMaximum());
3759 line->DrawLine(fPH->GetParameter(1)+fPH->GetParameter(2)+fPH->GetParameter(3)
3761 ,fPH->GetParameter(1)+fPH->GetParameter(2)+fPH->GetParameter(3)
3762 ,projPH->GetMaximum());
3765 if (fPH->GetParameter(3) != 0) {
3766 fNumberFitSuccess++;
3767 fCurrentCoef[0] = kDrWidth / (fPH->GetParameter(3));
3768 fCurrentCoefE = (fPH->GetParError(3)/fPH->GetParameter(3))*fCurrentCoef[0];
3769 fCurrentCoef2[0] = fPH->GetParameter(1);
3770 fCurrentCoefE2 = fPH->GetParError(1);
3773 fCurrentCoef[0] = -TMath::Abs(fCurrentCoef[1]);
3774 fCurrentCoef2[0] = fCurrentCoef2[1];
3780 // Put the default value
3781 fCurrentCoef[0] = -TMath::Abs(fCurrentCoef[1]);
3782 fCurrentCoef2[0] = fCurrentCoef2[1];
3785 if (fDebugLevel != 1) {
3790 //_____________________________________________________________________________
3791 Bool_t AliTRDCalibraFit::FitPRFGausMI(Double_t *arraye, Double_t *arraym, Double_t *arrayme, Int_t nBins, Float_t xMin, Float_t xMax)
3794 // Fit methode for the sigma of the pad response function
3799 fCurrentCoef[0] = 0.0;
3800 fCurrentCoefE = 0.0;
3802 Double_t ret = FitGausMI(arraye, arraym, arrayme, nBins, xMin, xMax,¶m);
3805 fCurrentCoef[0] = -fCurrentCoef[1];
3809 fNumberFitSuccess++;
3810 fCurrentCoef[0] = param[2];
3811 fCurrentCoefE = ret;
3815 //_____________________________________________________________________________
3816 Double_t AliTRDCalibraFit::FitGausMI(Double_t *arraye, Double_t *arraym, Double_t *arrayme, Int_t nBins, Float_t xMin, Float_t xMax, TVectorD *param, Bool_t bError)
3819 // Fit methode for the sigma of the pad response function
3822 //We should have at least 3 points
3823 if(nBins <=3) return -4.0;
3825 TLinearFitter fitter(3,"pol2");
3826 fitter.StoreData(kFALSE);
3827 fitter.ClearPoints();
3829 Float_t binWidth = (xMax-xMin)/(Float_t)nBins;
3830 Float_t entries = 0;
3831 Int_t nbbinwithentries = 0;
3832 for (Int_t i=0; i<nBins; i++){
3834 if(arraye[i] > 15) nbbinwithentries++;
3835 //printf("entries for i %d: %f\n",i,arraye[i]);
3837 if ((entries<700) || (nbbinwithentries < ((Int_t)(nBins/2)))) return -4;
3838 //printf("entries %f\n",entries);
3839 //printf("nbbinwithentries %d\n",nbbinwithentries);
3842 Float_t errorm = 0.0;
3843 Float_t errorn = 0.0;
3844 Float_t error = 0.0;
3847 for (Int_t ibin=0;ibin<nBins; ibin++){
3848 Float_t entriesI = arraye[ibin];
3849 Float_t valueI = arraym[ibin];
3850 Double_t xcenter = 0.0;
3852 if ((entriesI>15) && (valueI>0.0)){
3853 xcenter = xMin+(ibin+0.5)*binWidth;
3858 if((valueI + 0.01) > 0.0) errorm = TMath::Log((valueI + 0.01)/valueI);
3859 if((valueI - 0.01) > 0.0) errorn = TMath::Log((valueI - 0.01)/valueI);
3860 error = TMath::Max(TMath::Abs(errorm),TMath::Abs(errorn));
3863 if((valueI + arrayme[ibin]) > 0.0) errorm = TMath::Log((valueI + arrayme[ibin])/valueI);
3864 if((valueI - arrayme[ibin]) > 0.0) errorn = TMath::Log((valueI - arrayme[ibin])/valueI);
3865 error = TMath::Max(TMath::Abs(errorm),TMath::Abs(errorn));
3867 if(error == 0.0) continue;
3868 val = TMath::Log(Float_t(valueI));
3869 fitter.AddPoint(&xcenter,val,error);
3873 if(fDebugLevel > 1){
3875 if ( !fDebugStreamer ) {
3877 TDirectory *backup = gDirectory;
3878 fDebugStreamer = new TTreeSRedirector("TRDDebugFitPRF.root");
3879 if ( backup ) backup->cd(); //we don't want to be cd'd to the debug streamer
3882 Int_t detector = fCountDet;
3883 Int_t layer = GetLayer(fCountDet);
3886 (* fDebugStreamer) << "FitGausMIFill"<<
3887 "detector="<<detector<<
3891 "entriesI="<<entriesI<<
3894 "xcenter="<<xcenter<<
3904 if(npoints <=3) return -4.0;
3909 fitter.GetParameters(par);
3910 chi2 = fitter.GetChisquare()/Float_t(npoints);
3913 if (!param) param = new TVectorD(3);
3914 if(par[2] == 0.0) return -4.0;
3915 Double_t x = TMath::Sqrt(TMath::Abs(-2*par[2]));
3916 Double_t deltax = (fitter.GetParError(2))/x;
3917 Double_t errorparam2 = TMath::Abs(deltax)/(x*x);
3920 (*param)[1] = par[1]/(-2.*par[2]);
3921 (*param)[2] = 1./TMath::Sqrt(TMath::Abs(-2.*par[2]));
3922 Double_t lnparam0 = par[0]+ par[1]* (*param)[1] + par[2]*(*param)[1]*(*param)[1];
3923 if ( lnparam0>307 ) return -4;
3924 (*param)[0] = TMath::Exp(lnparam0);
3926 if(fDebugLevel > 1){
3928 if ( !fDebugStreamer ) {
3930 TDirectory *backup = gDirectory;
3931 fDebugStreamer = new TTreeSRedirector("TRDDebugFitPRF.root");
3932 if ( backup ) backup->cd(); //we don't want to be cd'd to the debug streamer
3935 Int_t detector = fCountDet;
3936 Int_t layer = GetLayer(fCountDet);
3939 (* fDebugStreamer) << "FitGausMIFit"<<
3940 "detector="<<detector<<
3943 "errorsigma="<<chi2<<
3944 "mean="<<(*param)[1]<<
3945 "sigma="<<(*param)[2]<<
3946 "constant="<<(*param)[0]<<
3951 if((chi2/(*param)[2]) > 0.1){
3953 chi2 = FitGausMI(arraye,arraym,arrayme,nBins,xMin,xMax,param,kFALSE);
3958 if(fDebugLevel == 1){
3959 TString name("PRF");
3960 name += (Int_t)xMin;
3961 name += (Int_t)xMax;
3962 TCanvas *c1 = new TCanvas((const char *)name,(const char *)name,50,50,600,800);
3965 TH1F *histo = new TH1F((const char *)name,(const char *)name,nBins,xMin,xMax);
3966 for(Int_t k = 0; k < nBins; k++){
3967 histo->SetBinContent(k+1,arraym[k]);
3968 histo->SetBinError(k+1,arrayme[k]);
3971 name += "functionf";
3972 TF1 *f1= new TF1((const char*)name,"[0]*exp(-(x-[1])^2/(2*[2]*[2]))",xMin,xMax);
3973 f1->SetParameter(0, (*param)[0]);
3974 f1->SetParameter(1, (*param)[1]);
3975 f1->SetParameter(2, (*param)[2]);
3983 //_____________________________________________________________________________
3984 void AliTRDCalibraFit::FitPRF(TH1 *projPRF)
3987 // Fit methode for the sigma of the pad response function
3990 fCurrentCoef[0] = 0.0;
3991 fCurrentCoefE = 0.0;
3993 if (fDebugLevel != 1) {
3994 projPRF->Fit("gaus","0M","",-fRangeFitPRF,fRangeFitPRF);
3997 TCanvas *cfit = new TCanvas("cfit","cfit",50,50,600,800);
3999 projPRF->Fit("gaus","M+","",-fRangeFitPRF,fRangeFitPRF);
4002 fCurrentCoef[0] = projPRF->GetFunction("gaus")->GetParameter(2);
4003 fCurrentCoefE = projPRF->GetFunction("gaus")->GetParError(2);
4004 if(fCurrentCoef[0] <= 0.0) fCurrentCoef[0] = -fCurrentCoef[1];
4006 fNumberFitSuccess++;
4009 //_____________________________________________________________________________
4010 void AliTRDCalibraFit::RmsPRF(TH1 *projPRF)
4013 // Fit methode for the sigma of the pad response function
4015 fCurrentCoef[0] = 0.0;
4016 fCurrentCoefE = 0.0;
4017 if (fDebugLevel == 1) {
4018 TCanvas *cfit = new TCanvas("cfit","cfit",50,50,600,800);
4022 fCurrentCoef[0] = projPRF->GetRMS();
4023 if(fCurrentCoef[0] <= 0.0) fCurrentCoef[0] = -fCurrentCoef[1];
4025 fNumberFitSuccess++;
4028 //_____________________________________________________________________________
4029 void AliTRDCalibraFit::FitTnpRange(Double_t *arraye, Double_t *arraym, Double_t *arrayme, Int_t nbg, Int_t nybins)
4032 // Fit methode for the sigma of the pad response function with 2*nbg tan bins
4035 TLinearFitter linearfitter = TLinearFitter(3,"pol2");
4038 Int_t nbins = (Int_t)(nybins/(2*nbg));
4039 Float_t lowedge = -3.0*nbg;
4040 Float_t upedge = lowedge + 3.0;
4043 Double_t xvalues = -0.2*nbg+0.1;
4045 Int_t total = 2*nbg;
4048 for(Int_t k = 0; k < total; k++){
4049 if(FitPRFGausMI(arraye+offset, arraym+offset, arrayme+offset, nbins, lowedge, upedge)){
4051 y = fCurrentCoef[0]*fCurrentCoef[0];
4052 linearfitter.AddPoint(&xvalues,y,2*fCurrentCoefE*fCurrentCoef[0]);
4055 if(fDebugLevel > 1){
4057 if ( !fDebugStreamer ) {
4059 TDirectory *backup = gDirectory;
4060 fDebugStreamer = new TTreeSRedirector("TRDDebugFitPRF.root");
4061 if ( backup ) backup->cd(); //we don't want to be cd'd to the debug streamer
4064 Int_t detector = fCountDet;
4065 Int_t layer = GetLayer(fCountDet);
4066 Int_t nbtotal = total;
4068 Float_t low = lowedge;
4069 Float_t up = upedge;
4070 Float_t tnp = xvalues;
4071 Float_t wid = fCurrentCoef[0];
4072 Float_t widfE = fCurrentCoefE;
4074 (* fDebugStreamer) << "FitTnpRange0"<<
4075 "detector="<<detector<<
4077 "nbtotal="<<nbtotal<<
4095 fCurrentCoefE = 0.0;
4096 fCurrentCoef[0] = 0.0;
4098 //printf("npoints\n",npoints);
4101 fCurrentCoef[0] = -TMath::Abs(fCurrentCoef[1]);
4106 linearfitter.Eval();
4107 linearfitter.GetParameters(pars0);
4108 Double_t pointError0 = TMath::Sqrt(linearfitter.GetChisquare()/npoints);
4109 Double_t errorsx0 = linearfitter.GetParError(2)*pointError0;
4110 Double_t min0 = 0.0;
4111 Double_t ermin0 = 0.0;
4112 //Double_t prfe0 = 0.0;
4113 Double_t prf0 = 0.0;
4114 if((pars0[2] > 0.0) && (pars0[1] != 0.0)) {
4115 min0 = -pars0[1]/(2*pars0[2]);
4116 ermin0 = TMath::Abs(min0*(errorsx0/pars0[2]+linearfitter.GetParError(1)*pointError0/pars0[1]));
4117 prf0 = pars0[0]+pars0[1]*min0+pars0[2]*min0*min0;
4120 prfe0 = linearfitter->GetParError(0)*pointError0
4121 +(linearfitter->GetParError(1)*pointError0/pars0[1]+ermin0/min0)*pars0[1]*min0
4122 +(linearfitter->GetParError(2)*pointError0/pars0[2]+2*ermin0/min0)*pars0[2]*min0*min0;
4123 prfe0 = prfe0/(2*TMath::Sqrt(prf0));
4124 fCurrentCoefE = (Float_t) prfe0;
4126 fCurrentCoef[0] = (Float_t) TMath::Sqrt(TMath::Abs(prf0));
4129 fCurrentCoef[0] = -TMath::Abs(fCurrentCoef[1]);
4133 fCurrentCoef[0] = -TMath::Abs(fCurrentCoef[1]);
4136 if(fDebugLevel > 1){
4138 if ( !fDebugStreamer ) {
4140 TDirectory *backup = gDirectory;
4141 fDebugStreamer = new TTreeSRedirector("TRDDebugFitPRF.root");
4142 if ( backup ) backup->cd(); //we don't want to be cd'd to the debug streamer
4145 Int_t detector = fCountDet;
4146 Int_t layer = GetLayer(fCountDet);
4147 Int_t nbtotal = total;
4148 Double_t colsize[6] = {0.635,0.665,0.695,0.725,0.755,0.785};
4149 Double_t sigmax = TMath::Sqrt(TMath::Abs(pars0[2]))*10000*colsize[layer];
4151 (* fDebugStreamer) << "FitTnpRange1"<<
4152 "detector="<<detector<<
4154 "nbtotal="<<nbtotal<<
4158 "npoints="<<npoints<<
4161 "sigmaprf="<<fCurrentCoef[0]<<
4162 "sigprf="<<fCurrentCoef[1]<<
4169 //_____________________________________________________________________________
4170 void AliTRDCalibraFit::FitMean(TH1 *projch, Double_t nentries, Double_t mean)
4173 // Only mean methode for the gain factor
4176 fCurrentCoef[0] = mean;
4177 fCurrentCoefE = 0.0;
4178 if(nentries > 0) fCurrentCoefE = projch->GetRMS()/TMath::Sqrt(nentries);
4179 if (fDebugLevel == 1) {
4180 TCanvas *cpmean = new TCanvas("cpmean","cpmean",50,50,600,800);
4184 CalculChargeCoefMean(kTRUE);
4185 fNumberFitSuccess++;
4187 //_____________________________________________________________________________
4188 void AliTRDCalibraFit::FitMeanW(TH1 *projch, Double_t nentries)
4191 // mean w methode for the gain factor
4195 Int_t nybins = projch->GetNbinsX();
4197 //The weight function
4198 Double_t a = 0.00228515;
4199 Double_t b = -0.00231487;
4200 Double_t c = 0.00044298;
4201 Double_t d = -0.00379239;
4202 Double_t e = 0.00338349;
4212 //A arbitrary error for the moment
4213 fCurrentCoefE = 0.0;
4214 fCurrentCoef[0] = 0.0;
4217 Double_t sumw = 0.0;
4219 Float_t sumAll = (Float_t) nentries;
4220 Int_t sumCurrent = 0;
4221 for(Int_t k = 0; k <nybins; k++){
4222 Double_t fraction = Float_t(sumCurrent)/Float_t(sumAll);
4223 if (fraction>0.95) break;
4224 Double_t weight = a + b*fraction + c*fraction*fraction + d *fraction*fraction*fraction+
4225 e*fraction*fraction*fraction*fraction;
4226 sumw += weight*projch->GetBinContent(k+1)*projch->GetBinCenter(k+1);
4227 sum += weight*projch->GetBinContent(k+1);
4228 sumCurrent += (Int_t) projch->GetBinContent(k+1);
4229 //printf("fraction %f, weight %f, bincontent %f\n",fraction,weight,projch->GetBinContent(k+1));
4231 if(sum > 0.0) fCurrentCoef[0] = (sumw/sum);
4233 if (fDebugLevel == 1) {
4234 TCanvas *cpmeanw = new TCanvas("cpmeanw","cpmeanw",50,50,600,800);
4238 fNumberFitSuccess++;
4239 CalculChargeCoefMean(kTRUE);
4241 //_____________________________________________________________________________
4242 void AliTRDCalibraFit::FitMeanWSm(TH1 *projch, Float_t sumAll)
4245 // mean w methode for the gain factor
4249 Int_t nybins = projch->GetNbinsX();
4251 //The weight function
4252 Double_t a = 0.00228515;
4253 Double_t b = -0.00231487;
4254 Double_t c = 0.00044298;
4255 Double_t d = -0.00379239;
4256 Double_t e = 0.00338349;
4266 //A arbitrary error for the moment
4267 fCurrentCoefE = 0.0;
4268 fCurrentCoef[0] = 0.0;
4271 Double_t sumw = 0.0;
4273 Int_t sumCurrent = 0;
4274 for(Int_t k = 0; k <nybins; k++){
4275 Double_t fraction = Float_t(sumCurrent)/Float_t(sumAll);
4276 if (fraction>0.95) break;
4277 Double_t weight = a + b*fraction + c*fraction*fraction + d *fraction*fraction*fraction+
4278 e*fraction*fraction*fraction*fraction;
4279 sumw += weight*projch->GetBinContent(k+1)*projch->GetBinCenter(k+1);
4280 sum += weight*projch->GetBinContent(k+1);
4281 sumCurrent += (Int_t) projch->GetBinContent(k+1);
4282 //printf("fraction %f, weight %f, bincontent %f\n",fraction,weight,projch->GetBinContent(k+1));
4284 if(sum > 0.0) fCurrentCoef[0] = (sumw/sum);
4286 if (fDebugLevel == 1) {
4287 TCanvas *cpmeanw = new TCanvas("cpmeanw","cpmeanw",50,50,600,800);
4291 fNumberFitSuccess++;
4293 //_____________________________________________________________________________
4294 void AliTRDCalibraFit::FitCH(TH1 *projch, Double_t mean)
4297 // Fit methode for the gain factor
4300 fCurrentCoef[0] = 0.0;
4301 fCurrentCoefE = 0.0;
4302 Double_t chisqrl = 0.0;
4303 Double_t chisqrg = 0.0;
4304 Double_t chisqr = 0.0;
4305 TF1 *fLandauGaus = new TF1("fLandauGaus",FuncLandauGaus,0,300,5);
4307 projch->Fit("landau","0",""
4308 ,(Double_t) mean/fBeginFitCharge
4309 ,projch->GetBinCenter(projch->GetNbinsX()));
4310 Double_t l3P0 = projch->GetFunction("landau")->GetParameter(0);
4311 Double_t l3P1 = projch->GetFunction("landau")->GetParameter(1);
4312 Double_t l3P2 = projch->GetFunction("landau")->GetParameter(2);
4313 chisqrl = projch->GetFunction("landau")->GetChisquare();
4315 projch->Fit("gaus","0",""
4316 ,(Double_t) mean/fBeginFitCharge
4317 ,projch->GetBinCenter(projch->GetNbinsX()));
4318 Double_t g3P0 = projch->GetFunction("gaus")->GetParameter(0);
4319 Double_t g3P2 = projch->GetFunction("gaus")->GetParameter(2);
4320 chisqrg = projch->GetFunction("gaus")->GetChisquare();
4322 fLandauGaus->SetParameters(l3P0,l3P1,l3P2,g3P0,g3P2);
4323 if (fDebugLevel != 1) {
4324 projch->Fit("fLandauGaus","0",""
4325 ,(Double_t) mean/fBeginFitCharge
4326 ,projch->GetBinCenter(projch->GetNbinsX()));
4327 chisqr = projch->GetFunction("fLandauGaus")->GetChisquare();
4330 TCanvas *cp = new TCanvas("cp","cp",50,50,600,800);
4332 projch->Fit("fLandauGaus","+",""
4333 ,(Double_t) mean/fBeginFitCharge
4334 ,projch->GetBinCenter(projch->GetNbinsX()));
4335 chisqr = projch->GetFunction("fLandauGaus")->GetChisquare();
4337 fLandauGaus->Draw("same");
4340 if ((projch->GetFunction("fLandauGaus")->GetParameter(1) > 0) && (projch->GetFunction("fLandauGaus")->GetParError(1) < (0.05*projch->GetFunction("fLandauGaus")->GetParameter(1))) && (chisqr < chisqrl) && (chisqr < chisqrg)) {
4341 //if ((projch->GetFunction("fLandauGaus")->GetParameter(1) > 0) && (chisqr < chisqrl) && (chisqr < chisqrg)) {
4342 fNumberFitSuccess++;
4343 CalculChargeCoefMean(kTRUE);
4344 fCurrentCoef[0] = projch->GetFunction("fLandauGaus")->GetParameter(1);
4345 fCurrentCoefE = projch->GetFunction("fLandauGaus")->GetParError(1);
4348 CalculChargeCoefMean(kFALSE);
4349 fCurrentCoef[0] = -TMath::Abs(fCurrentCoef[1]);
4352 if (fDebugLevel != 1) {
4357 //_____________________________________________________________________________
4358 void AliTRDCalibraFit::FitBisCH(TH1* projch, Double_t mean)
4361 // Fit methode for the gain factor more time consuming
4365 //Some parameters to initialise
4366 Double_t widthLandau, widthGaus, mPV, integral;
4367 Double_t chisquarel = 0.0;
4368 Double_t chisquareg = 0.0;
4369 projch->Fit("landau","0M+",""
4371 ,projch->GetBinCenter(projch->GetNbinsX()));
4372 widthLandau = projch->GetFunction("landau")->GetParameter(2);
4373 chisquarel = projch->GetFunction("landau")->GetChisquare();
4374 projch->Fit("gaus","0M+",""
4376 ,projch->GetBinCenter(projch->GetNbinsX()));
4377 widthGaus = projch->GetFunction("gaus")->GetParameter(2);
4378 chisquareg = projch->GetFunction("gaus")->GetChisquare();
4380 mPV = (projch->GetFunction("landau")->GetParameter(1))/2;
4381 integral = (projch->GetFunction("gaus")->Integral(0.3*mean,3*mean)+projch->GetFunction("landau")->Integral(0.3*mean,3*mean))/2;
4383 // Setting fit range and start values
4385 //Double_t sv[4] = { l3P2, fChargeCoef[1], projch->Integral("width"), fG3P2 };
4386 //Double_t sv[4] = { fL3P2, fChargeCoef[1], fL3P0, fG3P2 };
4387 Double_t sv[4] = { widthLandau, mPV, integral, widthGaus};
4388 Double_t pllo[4] = { 0.001, 0.001, projch->Integral()/3, 0.001};
4389 Double_t plhi[4] = { 300.0, 300.0, 30*projch->Integral(), 300.0};
4390 Double_t fp[4] = { 1.0, 1.0, 1.0, 1.0 };
4391 Double_t fpe[4] = { 1.0, 1.0, 1.0, 1.0 };
4394 fCurrentCoef[0] = 0.0;
4395 fCurrentCoefE = 0.0;
4399 TF1 *fitsnr = LanGauFit(projch,&fr[0],&sv[0]
4404 Double_t projchPeak;
4405 Double_t projchFWHM;
4406 LanGauPro(fp,projchPeak,projchFWHM);
4408 if ((fp[1] > 0) && ((fpe[1] < (0.05*fp[1])) && (chisqr < chisquarel) && (chisqr < chisquareg))) {
4409 //if ((fp[1] > 0) && ((chisqr < chisquarel) && (chisqr < chisquareg))) {
4410 fNumberFitSuccess++;
4411 CalculChargeCoefMean(kTRUE);
4412 fCurrentCoef[0] = fp[1];
4413 fCurrentCoefE = fpe[1];
4414 //chargeCoefE2 = chisqr;
4417 CalculChargeCoefMean(kFALSE);
4418 fCurrentCoef[0] = -TMath::Abs(fCurrentCoef[1]);
4420 if (fDebugLevel == 1) {
4421 AliInfo(Form("fChargeCoef[0]: %f",(Float_t) fCurrentCoef[0]));
4422 TCanvas *cpy = new TCanvas("cpy","cpy",50,50,600,800);
4425 fitsnr->Draw("same");
4431 //_____________________________________________________________________________
4432 Double_t *AliTRDCalibraFit::CalculPolynomeLagrange2(Double_t *x, Double_t *y) const
4435 // Calcul the coefficients of the polynome passant par ces trois points de degre 2
4437 Double_t *c = new Double_t[5];
4438 Double_t x0 = y[0]/((x[0]-x[1])*(x[0]-x[2]));
4439 Double_t x1 = y[1]/((x[1]-x[0])*(x[1]-x[2]));
4440 Double_t x2 = y[2]/((x[2]-x[0])*(x[2]-x[1]));
4445 c[1] = -(x0*(x[1]+x[2])+x1*(x[0]+x[2])+x2*(x[0]+x[1]));
4446 c[0] = x0*x[1]*x[2]+x1*x[0]*x[2]+x2*x[0]*x[1];
4453 //_____________________________________________________________________________
4454 Double_t *AliTRDCalibraFit::CalculPolynomeLagrange3(Double_t *x, Double_t *y) const
4457 // Calcul the coefficients of the polynome passant par ces quatre points de degre 3
4459 Double_t *c = new Double_t[5];
4460 Double_t x0 = y[0]/((x[0]-x[1])*(x[0]-x[2])*(x[0]-x[3]));
4461 Double_t x1 = y[1]/((x[1]-x[0])*(x[1]-x[2])*(x[1]-x[3]));
4462 Double_t x2 = y[2]/((x[2]-x[0])*(x[2]-x[1])*(x[2]-x[3]));
4463 Double_t x3 = y[3]/((x[3]-x[0])*(x[3]-x[1])*(x[3]-x[2]));
4467 c[2] = -(x0*(x[1]+x[2]+x[3])
4468 +x1*(x[0]+x[2]+x[3])
4469 +x2*(x[0]+x[1]+x[3])
4470 +x3*(x[0]+x[1]+x[2]));
4471 c[1] = (x0*(x[1]*x[2]+x[1]*x[3]+x[2]*x[3])
4472 +x1*(x[0]*x[2]+x[0]*x[3]+x[2]*x[3])
4473 +x2*(x[0]*x[1]+x[0]*x[3]+x[1]*x[3])
4474 +x3*(x[0]*x[1]+x[0]*x[2]+x[1]*x[2]));
4476 c[0] = -(x0*x[1]*x[2]*x[3]
4479 +x3*x[0]*x[1]*x[2]);
4487 //_____________________________________________________________________________
4488 Double_t *AliTRDCalibraFit::CalculPolynomeLagrange4(Double_t *x, Double_t *y) const
4491 // Calcul the coefficients of the polynome passant par ces cinqs points de degre 4
4493 Double_t *c = new Double_t[5];
4494 Double_t x0 = y[0]/((x[0]-x[1])*(x[0]-x[2])*(x[0]-x[3])*(x[0]-x[4]));
4495 Double_t x1 = y[1]/((x[1]-x[0])*(x[1]-x[2])*(x[1]-x[3])*(x[1]-x[4]));
4496 Double_t x2 = y[2]/((x[2]-x[0])*(x[2]-x[1])*(x[2]-x[3])*(x[2]-x[4]));
4497 Double_t x3 = y[3]/((x[3]-x[0])*(x[3]-x[1])*(x[3]-x[2])*(x[3]-x[4]));
4498 Double_t x4 = y[4]/((x[4]-x[0])*(x[4]-x[1])*(x[4]-x[2])*(x[4]-x[3]));
4501 c[4] = x0+x1+x2+x3+x4;
4502 c[3] = -(x0*(x[1]+x[2]+x[3]+x[4])
4503 +x1*(x[0]+x[2]+x[3]+x[4])
4504 +x2*(x[0]+x[1]+x[3]+x[4])
4505 +x3*(x[0]+x[1]+x[2]+x[4])
4506 +x4*(x[0]+x[1]+x[2]+x[3]));
4507 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])
4508 +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])
4509 +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])
4510 +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])
4511 +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]));
4513 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])
4514 +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])
4515 +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])
4516 +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])
4517 +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]));
4519 c[0] = (x0*x[1]*x[2]*x[3]*x[4]
4520 +x1*x[0]*x[2]*x[3]*x[4]
4521 +x2*x[0]*x[1]*x[3]*x[4]
4522 +x3*x[0]*x[1]*x[2]*x[4]
4523 +x4*x[0]*x[1]*x[2]*x[3]);
4529 //_____________________________________________________________________________
4530 void AliTRDCalibraFit::NormierungCharge()
4533 // Normalisation of the gain factor resulting for the fits
4536 // Calcul of the mean of choosen method by fFitChargeNDB
4538 //printf("total number of entries %d\n",fVectorFitCH->GetEntriesFast());
4539 for (Int_t k = 0; k < (Int_t) fVectorFit.GetEntriesFast(); k++) {
4541 Int_t detector = ((AliTRDFitInfo *) fVectorFit.At(k))->GetDetector();
4542 Float_t *coef = ((AliTRDFitInfo *) fVectorFit.At(k))->GetCoef();
4543 //printf("detector %d coef[0] %f\n",detector,coef[0]);
4544 if (GetStack(detector) == 2) {
4547 if (GetStack(detector) != 2) {
4550 for (Int_t j = 0; j < total; j++) {
4558 fScaleFitFactor = fScaleFitFactor / sum;
4561 fScaleFitFactor = 1.0;
4564 //methode de boeuf mais bon...
4565 Double_t scalefactor = fScaleFitFactor;
4567 if(fDebugLevel > 1){
4569 if ( !fDebugStreamer ) {
4571 TDirectory *backup = gDirectory;
4572 fDebugStreamer = new TTreeSRedirector("TRDDebugFitCH.root");
4573 if ( backup ) backup->cd(); //we don't want to be cd'd to the debug streamer
4575 (* fDebugStreamer) << "NormierungCharge"<<
4576 "scalefactor="<<scalefactor<<
4580 //_____________________________________________________________________________
4581 TH1I *AliTRDCalibraFit::ReBin(TH1I *hist) const
4584 // Rebin of the 1D histo for the gain calibration if needed.
4585 // you have to choose fRebin, divider of fNumberBinCharge
4588 TAxis *xhist = hist->GetXaxis();
4589 TH1I *rehist = new TH1I("projrebin","",(Int_t) xhist->GetNbins()/fRebin
4590 ,xhist->GetBinLowEdge(1)
4591 ,xhist->GetBinUpEdge(xhist->GetNbins()));
4593 AliInfo(Form("fRebin: %d",fRebin));
4595 for (Int_t k = 1; k <= (Int_t) xhist->GetNbins()/fRebin; k++) {
4597 for (Int_t ji = i; ji < i+fRebin; ji++) {
4598 sum += hist->GetBinContent(ji);
4601 rehist->SetBinContent(k,sum);
4609 //_____________________________________________________________________________
4610 TH1F *AliTRDCalibraFit::ReBin(TH1F *hist) const
4613 // Rebin of the 1D histo for the gain calibration if needed
4614 // you have to choose fRebin divider of fNumberBinCharge
4617 TAxis *xhist = hist->GetXaxis();
4618 TH1F *rehist = new TH1F("projrebin","",(Int_t) xhist->GetNbins()/fRebin
4619 ,xhist->GetBinLowEdge(1)
4620 ,xhist->GetBinUpEdge(xhist->GetNbins()));
4622 AliInfo(Form("fRebin: %d",fRebin));
4624 for (Int_t k = 1; k <= (Int_t) xhist->GetNbins()/fRebin; k++) {
4626 for (Int_t ji = i; ji < i+fRebin; ji++) {
4627 sum += hist->GetBinContent(ji);
4630 rehist->SetBinContent(k,sum);
4638 //_____________________________________________________________________________
4639 TH1F *AliTRDCalibraFit::CorrectTheError(TGraphErrors *hist)
4642 // In the case of the vectors method the trees contains TGraphErrors for PH and PRF
4643 // to be able to add them after
4644 // We convert it to a TH1F to be able to applied the same fit function method
4645 // After having called this function you can not add the statistics anymore
4650 Int_t nbins = hist->GetN();
4651 Double_t *x = hist->GetX();
4652 Double_t *entries = hist->GetEX();
4653 Double_t *mean = hist->GetY();
4654 Double_t *square = hist->GetEY();
4655 fEntriesCurrent = 0;
4661 Double_t step = x[1] - x[0];
4662 Double_t minvalue = x[0] - step/2;
4663 Double_t maxvalue = x[(nbins-1)] + step/2;
4665 rehist = new TH1F("projcorrecterror","",nbins,minvalue,maxvalue);
4667 for (Int_t k = 0; k < nbins; k++) {
4668 rehist->SetBinContent(k+1,mean[k]);
4669 if (entries[k] > 0.0) {
4670 fEntriesCurrent += (Int_t) entries[k];
4671 Double_t d = TMath::Abs(square[k] - (mean[k]*mean[k]));
4672 rehist->SetBinError(k+1,TMath::Sqrt(d/entries[k]));
4675 rehist->SetBinError(k+1,0.0);
4679 if(fEntriesCurrent > 0) fNumberEnt++;
4685 //____________Some basic geometry function_____________________________________
4688 //_____________________________________________________________________________
4689 Int_t AliTRDCalibraFit::GetLayer(Int_t d) const
4692 // Reconstruct the plane number from the detector number
4695 return ((Int_t) (d % 6));
4699 //_____________________________________________________________________________
4700 Int_t AliTRDCalibraFit::GetStack(Int_t d) const
4703 // Reconstruct the stack number from the detector number
4705 const Int_t kNlayer = 6;
4707 return ((Int_t) (d % 30) / kNlayer);
4711 //_____________________________________________________________________________
4712 Int_t AliTRDCalibraFit::GetSector(Int_t d) const
4715 // Reconstruct the sector number from the detector number
4719 return ((Int_t) (d / fg));
4724 //____________Fill and Init tree Gain, PRF, Vdrift and T0______________________
4726 //_______________________________________________________________________________
4727 void AliTRDCalibraFit::ResetVectorFit()
4730 // Reset the VectorFits
4733 fVectorFit.SetOwner();
4735 fVectorFit2.SetOwner();
4736 fVectorFit2.Clear();
4740 //____________Private Functions________________________________________________
4743 //_____________________________________________________________________________
4744 Double_t AliTRDCalibraFit::PH(Double_t *x, Double_t *par)
4747 // Function for the fit
4750 //TF1 *fAsymmGauss = new TF1("fAsymmGauss",AsymmGauss,0,4,6);
4752 //PARAMETERS FOR FIT PH
4754 //fAsymmGauss->SetParameter(0,0.113755);
4755 //fAsymmGauss->SetParameter(1,0.350706);
4756 //fAsymmGauss->SetParameter(2,0.0604244);
4757 //fAsymmGauss->SetParameter(3,7.65596);
4758 //fAsymmGauss->SetParameter(4,1.00124);
4759 //fAsymmGauss->SetParameter(5,0.870597); // No tail cancelation
4767 Double_t dx = 0.005;
4768 Double_t xs = par[1];
4770 Double_t paras[2] = { 0.0, 0.0 };
4773 if ((xs >= par[1]) &&
4774 (xs < (par[1]+par[2]))) {
4775 //fAsymmGauss->SetParameter(0,par[0]);
4776 //fAsymmGauss->SetParameter(1,xs);
4777 //ss += fAsymmGauss->Eval(xx);
4780 ss += AsymmGauss(&xx,paras);
4782 if ((xs >= (par[1]+par[2])) &&
4783 (xs < (par[1]+par[2]+par[3]))) {
4784 //fAsymmGauss->SetParameter(0,par[0]*par[4]);
4785 //fAsymmGauss->SetParameter(1,xs);
4786 //ss += fAsymmGauss->Eval(xx);
4787 paras[0] = par[0]*par[4];
4789 ss += AsymmGauss(&xx,paras);
4798 //_____________________________________________________________________________
4799 Double_t AliTRDCalibraFit::AsymmGauss(Double_t *x, Double_t *par)
4802 // Function for the fit
4805 //par[0] = normalization
4813 Double_t par1save = par[1];
4814 //Double_t par2save = par[2];
4815 Double_t par2save = 0.0604244;
4816 //Double_t par3save = par[3];
4817 Double_t par3save = 7.65596;
4818 //Double_t par5save = par[5];
4819 Double_t par5save = 0.870597;
4820 Double_t dx = x[0] - par1save;
4822 Double_t sigma2 = par2save*par2save;
4823 Double_t sqrt2 = TMath::Sqrt(2.0);
4824 Double_t exp1 = par3save * TMath::Exp(-par3save * (dx - 0.5 * par3save * sigma2))
4825 * (1.0 - TMath::Erf((par3save * sigma2 - dx) / (sqrt2 * par2save)));
4826 Double_t exp2 = par5save * TMath::Exp(-par5save * (dx - 0.5 * par5save * sigma2))
4827 * (1.0 - TMath::Erf((par5save * sigma2 - dx) / (sqrt2 * par2save)));
4829 //return par[0]*(exp1+par[4]*exp2);
4830 return par[0] * (exp1 + 1.00124 * exp2);
4834 //_____________________________________________________________________________
4835 Double_t AliTRDCalibraFit::FuncLandauGaus(Double_t *x, Double_t *par)
4838 // Sum Landau + Gaus with identical mean
4841 Double_t valLandau = par[0] * TMath::Landau(x[0],par[1],par[2]);
4842 //Double_t valGaus = par[3] * TMath::Gaus(x[0],par[4],par[5]);
4843 Double_t valGaus = par[3] * TMath::Gaus(x[0],par[1],par[4]);
4844 Double_t val = valLandau + valGaus;
4850 //_____________________________________________________________________________
4851 Double_t AliTRDCalibraFit::LanGauFun(Double_t *x, Double_t *par)
4854 // Function for the fit
4857 // par[0]=Width (scale) parameter of Landau density
4858 // par[1]=Most Probable (MP, location) parameter of Landau density
4859 // par[2]=Total area (integral -inf to inf, normalization constant)
4860 // par[3]=Width (sigma) of convoluted Gaussian function
4862 // In the Landau distribution (represented by the CERNLIB approximation),
4863 // the maximum is located at x=-0.22278298 with the location parameter=0.
4864 // This shift is corrected within this function, so that the actual
4865 // maximum is identical to the MP parameter.
4868 // Numeric constants
4869 Double_t invsq2pi = 0.3989422804014; // (2 pi)^(-1/2)
4870 Double_t mpshift = -0.22278298; // Landau maximum location
4872 // Control constants
4873 Double_t np = 100.0; // Number of convolution steps
4874 Double_t sc = 5.0; // Convolution extends to +-sc Gaussian sigmas
4886 // MP shift correction
4887 mpc = par[1] - mpshift * par[0];
4889 // Range of convolution integral
4890 xlow = x[0] - sc * par[3];
4891 xupp = x[0] + sc * par[3];
4893 step = (xupp - xlow) / np;
4895 // Convolution integral of Landau and Gaussian by sum
4896 for (i = 1.0; i <= np/2; i++) {
4898 xx = xlow + (i-.5) * step;
4899 fland = TMath::Landau(xx,mpc,par[0]) / par[0];
4900 sum += fland * TMath::Gaus(x[0],xx,par[3]);
4902 xx = xupp - (i-.5) * step;
4903 fland = TMath::Landau(xx,mpc,par[0]) / par[0];
4904 sum += fland * TMath::Gaus(x[0],xx,par[3]);
4908 return (par[2] * step * sum * invsq2pi / par[3]);
4911 //_____________________________________________________________________________
4912 TF1 *AliTRDCalibraFit::LanGauFit(TH1 *his, Double_t *fitrange, Double_t *startvalues
4913 , Double_t *parlimitslo, Double_t *parlimitshi
4914 , Double_t *fitparams, Double_t *fiterrors
4915 , Double_t *chiSqr, Int_t *ndf) const
4918 // Function for the fit
4922 Char_t funname[100];
4924 TF1 *ffitold = (TF1 *) gROOT->GetListOfFunctions()->FindObject(funname);
4929 TF1 *ffit = new TF1(funname,LanGauFun,fitrange[0],fitrange[1],4);
4930 ffit->SetParameters(startvalues);
4931 ffit->SetParNames("Width","MP","Area","GSigma");
4933 for (i = 0; i < 4; i++) {
4934 ffit->SetParLimits(i,parlimitslo[i],parlimitshi[i]);
4937 his->Fit(funname,"RB0"); // Fit within specified range, use ParLimits, do not plot
4939 ffit->GetParameters(fitparams); // Obtain fit parameters
4940 for (i = 0; i < 4; i++) {
4941 fiterrors[i] = ffit->GetParError(i); // Obtain fit parameter errors
4943 chiSqr[0] = ffit->GetChisquare(); // Obtain chi^2
4944 ndf[0] = ffit->GetNDF(); // Obtain ndf
4946 return (ffit); // Return fit function
4950 //_____________________________________________________________________________
4951 Int_t AliTRDCalibraFit::LanGauPro(Double_t *params, Double_t &maxx, Double_t &fwhm)
4954 // Function for the fit
4967 Int_t maxcalls = 10000;
4969 // Search for maximum
4970 p = params[1] - 0.1 * params[0];
4971 step = 0.05 * params[0];
4975 while ((l != lold) && (i < maxcalls)) {
4979 l = LanGauFun(&x,params);
4981 step = -step / 10.0;
4986 if (i == maxcalls) {
4992 // Search for right x location of fy
4993 p = maxx + params[0];
4999 while ( (l != lold) && (i < maxcalls) ) {
5004 l = TMath::Abs(LanGauFun(&x,params) - fy);
5018 // Search for left x location of fy
5020 p = maxx - 0.5 * params[0];
5026 while ((l != lold) && (i < maxcalls)) {
5030 l = TMath::Abs(LanGauFun(&x,params) - fy);
5032 step = -step / 10.0;
5037 if (i == maxcalls) {
5046 //_____________________________________________________________________________
5047 Double_t AliTRDCalibraFit::GausConstant(Double_t *x, Double_t *par)
5050 // Gaus with identical mean
5053 Double_t gauss = par[0] * TMath::Gaus(x[0],0.0,par[1])+par[2];