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: only the fit of the choosen calibration group fFitVoir (SetFitVoir)
32 // _fDebug = 2: a comparaison of the coefficients found and the default values
33 // in the choosen database.
34 // fCoef , histogram of the coefs as function of the calibration group number
35 // fDelta , histogram of the relative difference of the coef with the default
36 // value in the database as function of the calibration group number
37 // fError , dirstribution of this relative difference
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>
62 #include <TTreeStream.h>
63 #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 "AliTRDCommonParam.h"
80 #include "./Cal/AliTRDCalROC.h"
81 #include "./Cal/AliTRDCalPad.h"
82 #include "./Cal/AliTRDCalDet.h"
85 ClassImp(AliTRDCalibraFit)
87 AliTRDCalibraFit* AliTRDCalibraFit::fgInstance = 0;
88 Bool_t AliTRDCalibraFit::fgTerminated = kFALSE;
90 //_____________singleton implementation_________________________________________________
91 AliTRDCalibraFit *AliTRDCalibraFit::Instance()
94 // Singleton implementation
97 if (fgTerminated != kFALSE) {
101 if (fgInstance == 0) {
102 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) {
124 //______________________________________________________________________________________
125 AliTRDCalibraFit::AliTRDCalibraFit()
128 ,fNumberOfBinsExpected(0)
130 ,fBeginFitCharge(3.5)
132 ,fTakeTheMaxPH(kTRUE)
141 ,fNumberFitSuccess(0)
148 ,fCalibraMode(new AliTRDCalibraMode())
153 ,fScaleFitFactor(0.0)
162 ,fCurrentCoefDetector(0x0)
163 ,fCurrentCoefDetector2(0x0)
168 // Default constructor
171 fGeo = new AliTRDgeometry();
173 // Current variables initialised
174 for (Int_t k = 0; k < 2; k++) {
175 fCurrentCoef[k] = 0.0;
176 fCurrentCoef2[k] = 0.0;
178 for (Int_t i = 0; i < 3; i++) {
184 //______________________________________________________________________________________
185 AliTRDCalibraFit::AliTRDCalibraFit(const AliTRDCalibraFit &c)
188 ,fNumberOfBinsExpected(c.fNumberOfBinsExpected)
190 ,fBeginFitCharge(c.fBeginFitCharge)
191 ,fFitPHPeriode(c.fFitPHPeriode)
192 ,fTakeTheMaxPH(c.fTakeTheMaxPH)
193 ,fT0Shift0(c.fT0Shift0)
194 ,fT0Shift1(c.fT0Shift1)
195 ,fRangeFitPRF(c.fRangeFitPRF)
197 ,fMinEntries(c.fMinEntries)
199 ,fScaleGain(c.fScaleGain)
200 ,fNumberFit(c.fNumberFit)
201 ,fNumberFitSuccess(c.fNumberFitSuccess)
202 ,fNumberEnt(c.fNumberEnt)
203 ,fStatisticMean(c.fStatisticMean)
205 ,fDebugLevel(c.fDebugLevel)
206 ,fFitVoir(c.fFitVoir)
207 ,fMagneticField(c.fMagneticField)
209 ,fCurrentCoefE(c.fCurrentCoefE)
210 ,fCurrentCoefE2(c.fCurrentCoefE2)
213 ,fScaleFitFactor(c.fScaleFitFactor)
214 ,fEntriesCurrent(c.fEntriesCurrent)
215 ,fCountDet(c.fCountDet)
222 ,fCurrentCoefDetector(0x0)
223 ,fCurrentCoefDetector2(0x0)
231 if(c.fCalibraMode) fCalibraMode = new AliTRDCalibraMode(*c.fCalibraMode);
233 //Current variables initialised
234 for (Int_t k = 0; k < 2; k++) {
235 fCurrentCoef[k] = 0.0;
236 fCurrentCoef2[k] = 0.0;
238 for (Int_t i = 0; i < 3; i++) {
242 if(c.fCalDet) fCalDet = new AliTRDCalDet(*c.fCalDet);
243 if(c.fCalDet2) fCalDet2 = new AliTRDCalDet(*c.fCalDet2);
245 if(c.fCalROC) fCalROC = new AliTRDCalROC(*c.fCalROC);
246 if(c.fCalROC2) fCalROC = new AliTRDCalROC(*c.fCalROC2);
248 fVectorFit.SetName(c.fVectorFit.GetName());
249 for(Int_t k = 0; k < c.fVectorFit.GetEntriesFast(); k++){
250 AliTRDFitInfo *fitInfo = new AliTRDFitInfo();
251 Int_t detector = ((AliTRDFitInfo *)c.fVectorFit.UncheckedAt(k))->GetDetector();
253 if (GetStack(detector) == 2) {
259 Float_t *coef = new Float_t[ntotal];
260 for (Int_t i = 0; i < ntotal; i++) {
261 coef[i] = ((AliTRDFitInfo *)c.fVectorFit.UncheckedAt(k))->GetCoef()[i];
263 fitInfo->SetCoef(coef);
264 fitInfo->SetDetector(detector);
265 fVectorFit.Add((TObject *) fitInfo);
267 fVectorFit.SetName(c.fVectorFit.GetName());
268 for(Int_t k = 0; k < c.fVectorFit2.GetEntriesFast(); k++){
269 AliTRDFitInfo *fitInfo = new AliTRDFitInfo();
270 Int_t detector = ((AliTRDFitInfo *)c.fVectorFit2.UncheckedAt(k))->GetDetector();
272 if (GetStack(detector) == 2) {
278 Float_t *coef = new Float_t[ntotal];
279 for (Int_t i = 0; i < ntotal; i++) {
280 coef[i] = ((AliTRDFitInfo *)c.fVectorFit2.UncheckedAt(k))->GetCoef()[i];
282 fitInfo->SetCoef(coef);
283 fitInfo->SetDetector(detector);
284 fVectorFit2.Add((TObject *) fitInfo);
289 fGeo = new AliTRDgeometry();
292 //____________________________________________________________________________________
293 AliTRDCalibraFit::~AliTRDCalibraFit()
296 // AliTRDCalibraFit destructor
298 if ( fDebugStreamer ) delete fDebugStreamer;
299 if ( fCalDet ) delete fCalDet;
300 if ( fCalDet2 ) delete fCalDet2;
301 if ( fCalROC ) delete fCalROC;
302 if ( fCalROC2 ) delete fCalROC2;
303 if( fCurrentCoefDetector ) delete [] fCurrentCoefDetector;
304 if( fCurrentCoefDetector2 ) delete [] fCurrentCoefDetector2;
306 fVectorFit2.Delete();
312 //_____________________________________________________________________________
313 void AliTRDCalibraFit::Destroy()
325 //_____________________________________________________________________________
326 void AliTRDCalibraFit::DestroyDebugStreamer()
329 // Delete DebugStreamer
332 if ( fDebugStreamer ) delete fDebugStreamer;
333 fDebugStreamer = 0x0;
336 //__________________________________________________________________________________
337 void AliTRDCalibraFit::RangeChargeIntegration(Float_t vdrift, Float_t t0, Int_t &begin, Int_t &peak, Int_t &end) const
340 // From the drift velocity and t0
341 // return the position of the peak and maximum negative slope
344 const Float_t kDrWidth = AliTRDgeometry::DrThick(); // drift region
345 Double_t widbins = 0.1; // 0.1 mus
347 //peak and maxnegslope in mus
348 Double_t begind = t0*widbins + fT0Shift0;
349 Double_t peakd = t0*widbins + fT0Shift1;
350 Double_t maxnegslope = (kDrWidth + vdrift*peakd)/vdrift;
352 // peak and maxnegslope in timebin
353 begin = TMath::Nint(begind*widbins);
354 peak = TMath::Nint(peakd*widbins);
355 end = TMath::Nint(maxnegslope*widbins);
358 //____________Functions fit Online CH2d________________________________________
359 Bool_t AliTRDCalibraFit::AnalyseCH(const TH2I *ch)
362 // Fit the 1D histos, projections of the 2D ch on the Xaxis, for each
363 // calibration group normalized the resulted coefficients (to 1 normally)
366 // Set the calibration mode
367 //const char *name = ch->GetTitle();
368 TString name = ch->GetTitle();
369 if(!SetModeCalibration(name,0)) return kFALSE;
371 // Number of Ybins (detectors or groups of pads)
372 Int_t nbins = ch->GetNbinsX();// charge
373 Int_t nybins = ch->GetNbinsY();// groups number
374 if (!InitFit(nybins,0)) {
380 fStatisticMean = 0.0;
382 fNumberFitSuccess = 0;
384 // Init fCountDet and fCount
385 InitfCountDetAndfCount(0);
386 // Beginning of the loop betwwen dect1 and dect2
387 for (Int_t idect = fDect1; idect < fDect2; idect++) {
388 // Determination of fNnZ, fNnRphi, fNfragZ and fNfragRphi...
389 UpdatefCountDetAndfCount(idect,0);
390 ReconstructFitRowMinRowMax(idect, 0);
392 TH1I *projch = (TH1I *) ch->ProjectionX("projch",idect+1,idect+1,(Option_t *)"e");
393 projch->SetDirectory(0);
394 // Number of entries for this calibration group
395 Double_t nentries = 0.0;
397 for (Int_t k = 0; k < nbins; k++) {
398 Int_t binnb = (nbins+2)*(idect+1)+(k+1);
399 nentries += ch->GetBinContent(binnb);
400 mean += projch->GetBinCenter(k+1)*projch->GetBinContent(k+1);
401 projch->SetBinError(k+1,TMath::Sqrt(projch->GetBinContent(k+1)));
403 projch->SetEntries(nentries);
404 //printf("The number of entries for the group %d is %f\n",idect,nentries);
409 // Rebin and statistic stuff
411 projch = ReBin((TH1I *) projch);
413 // This detector has not enough statistics or was off
414 if (nentries <= fMinEntries) {
415 NotEnoughStatisticCH(idect);
416 if (fDebugLevel != 1) {
421 // Statistics of the group fitted
422 fStatisticMean += nentries;
427 case 0: FitMeanW((TH1 *) projch, nentries); break;
428 case 1: FitMean((TH1 *) projch, nentries, mean); break;
429 case 2: FitCH((TH1 *) projch, mean); break;
430 case 3: FitBisCH((TH1 *) projch, mean); break;
431 default: return kFALSE;
434 FillInfosFitCH(idect);
436 if (fDebugLevel != 1) {
441 if (fDebugLevel != 1) {
445 if (fNumberFit > 0) {
446 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));
447 fStatisticMean = fStatisticMean / fNumberFit;
450 AliInfo(Form("There are %d with at least one entries. There is no fit!",fNumberEnt));
452 delete fDebugStreamer;
453 fDebugStreamer = 0x0;
457 //____________Functions fit Online CH2d________________________________________
458 Bool_t AliTRDCalibraFit::AnalyseCH(AliTRDCalibraVector *calvect)
461 // Reconstruct a 1D histo from the vectorCH for each calibration group,
462 // fit the histo, normalized the resulted coefficients (to 1 normally)
465 // Set the calibraMode
466 //const char *name = calvect->GetNameCH();
467 TString name = calvect->GetNameCH();
468 if(!SetModeCalibration(name,0)) return kFALSE;
470 // Number of Xbins (detectors or groups of pads)
471 if (!InitFit((432*calvect->GetDetCha0(0)+108*calvect->GetDetCha2(0)),0)) {
477 fStatisticMean = 0.0;
479 fNumberFitSuccess = 0;
481 // Init fCountDet and fCount
482 InitfCountDetAndfCount(0);
483 // Beginning of the loop between dect1 and dect2
484 for (Int_t idect = fDect1; idect < fDect2; idect++) {
485 // Determination of fNnZ, fNnRphi, fNfragZ and fNfragRphi...........
486 UpdatefCountDetAndfCount(idect,0);
487 ReconstructFitRowMinRowMax(idect,0);
489 Double_t nentries = 0.0;
491 if(!calvect->GetCHEntries(fCountDet)) {
492 NotEnoughStatisticCH(idect);
498 TH1F *projch = calvect->ConvertVectorCHHisto(fCountDet,(idect-(fCount-(fCalibraMode->GetNfragZ(0)*fCalibraMode->GetNfragRphi(0)))),(const char *) tname);
499 projch->SetDirectory(0);
500 for (Int_t k = 0; k < calvect->GetNumberBinCharge(); k++) {
501 nentries += projch->GetBinContent(k+1);
502 mean += projch->GetBinCenter(k+1)*projch->GetBinContent(k+1);
508 //printf("The number of entries for the group %d is %f\n",idect,nentries);
511 projch = ReBin((TH1F *) projch);
513 // This detector has not enough statistics or was not found in VectorCH
514 if (nentries <= fMinEntries) {
515 NotEnoughStatisticCH(idect);
518 // Statistic of the histos fitted
519 fStatisticMean += nentries;
524 case 0: FitMeanW((TH1 *) projch, nentries); break;
525 case 1: FitMean((TH1 *) projch, nentries, mean); break;
526 case 2: FitCH((TH1 *) projch, mean); break;
527 case 3: FitBisCH((TH1 *) projch, mean); break;
528 default: return kFALSE;
531 FillInfosFitCH(idect);
534 if (fDebugLevel != 1) {
538 if (fNumberFit > 0) {
539 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));
540 fStatisticMean = fStatisticMean / fNumberFit;
543 AliInfo(Form("There are %d with at least one entries. There is no fit!",fNumberEnt));
545 delete fDebugStreamer;
546 fDebugStreamer = 0x0;
549 //____________Functions fit Online CH2d________________________________________
550 Double_t AliTRDCalibraFit::AnalyseCHAllTogether(const TH2I *ch)
553 // Fit the 1D histos, projections of the 2D ch on the Xaxis, for each
554 // calibration group normalized the resulted coefficients (to 1 normally)
556 Int_t nbins = ch->GetNbinsX();// charge
557 Int_t nybins = ch->GetNbinsY();// groups number
559 TH1I *projch = (TH1I *) ch->ProjectionX("projch",1,nybins+1,(Option_t *)"e");
560 projch->SetDirectory(0);
561 // Number of entries for this calibration group
562 Double_t nentries = 0.0;
564 for (Int_t k = 0; k < nbins; k++) {
565 nentries += projch->GetBinContent(k+1);
566 mean += projch->GetBinCenter(k+1)*projch->GetBinContent(k+1);
567 projch->SetBinError(k+1,TMath::Sqrt(projch->GetBinContent(k+1)));
569 projch->SetEntries(nentries);
570 //printf("The number of entries for the group %d is %f\n",idect,nentries);
575 // This detector has not enough statistics or was off
576 if (nentries <= fMinEntries) {
578 AliInfo(Form("There are %d entries for all together, it is not enough (%d)",(Int_t) nentries, fMinEntries));
584 case 0: FitMeanW((TH1 *) projch, nentries); break;
585 case 1: FitMean((TH1 *) projch, nentries, mean); break;
586 case 2: FitCH((TH1 *) projch, mean); break;
587 case 3: FitBisCH((TH1 *) projch, mean); break;
588 default: return kFALSE;
590 delete fDebugStreamer;
591 fDebugStreamer = 0x0;
593 if(fCurrentCoef[0] > 0.0) return fCurrentCoef[0];
597 //________________functions fit Online PH2d____________________________________
598 Bool_t AliTRDCalibraFit::AnalysePH(const TProfile2D *ph)
601 // Take the 1D profiles (average pulse height), projections of the 2D PH
602 // on the Xaxis, for each calibration group
603 // Reconstruct a drift velocity
604 // A first calibration of T0 is also made using the same method
607 // Set the calibration mode
608 //const char *name = ph->GetTitle();
609 TString name = ph->GetTitle();
610 if(!SetModeCalibration(name,1)) return kFALSE;
612 //printf("Mode calibration set\n");
614 // Number of Xbins (detectors or groups of pads)
615 Int_t nbins = ph->GetNbinsX();// time
616 Int_t nybins = ph->GetNbinsY();// calibration group
617 if (!InitFit(nybins,1)) {
621 //printf("Init fit\n");
627 //printf("Init fit PH\n");
629 fStatisticMean = 0.0;
631 fNumberFitSuccess = 0;
633 // Init fCountDet and fCount
634 InitfCountDetAndfCount(1);
635 //printf("Init Count Det and fCount %d, %d\n",fDect1,fDect2);
637 // Beginning of the loop
638 for (Int_t idect = fDect1; idect < fDect2; idect++) {
639 //printf("idect = %d\n",idect);
640 // Determination of fNnZ, fNnRphi, fNfragZ and fNfragRphi.......
641 UpdatefCountDetAndfCount(idect,1);
642 ReconstructFitRowMinRowMax(idect,1);
644 TH1D *projph = (TH1D *) ph->ProjectionX("projph",idect+1,idect+1,(Option_t *) "e");
645 projph->SetDirectory(0);
646 // Number of entries for this calibration group
647 Double_t nentries = 0;
648 for (Int_t k = 0; k < nbins; k++) {
649 Int_t binnb = (nbins+2)*(idect+1)+(k+1);
650 nentries += ph->GetBinEntries(binnb);
655 //printf("The number of entries for the group %d is %f\n",idect,nentries);
656 // This detector has not enough statistics or was off
657 if (nentries <= fMinEntries) {
658 //printf("Not enough statistic!\n");
659 NotEnoughStatisticPH(idect,nentries);
660 if (fDebugLevel != 1) {
665 // Statistics of the histos fitted
667 fStatisticMean += nentries;
668 // Calcul of "real" coef
669 CalculVdriftCoefMean();
672 //printf("Method\n");
675 case 0: FitLagrangePoly((TH1 *) projph); break;
676 case 1: FitPente((TH1 *) projph); break;
677 case 2: FitPH((TH1 *) projph,(Int_t) (idect - fDect1)); break;
678 default: return kFALSE;
680 // Fill the tree if end of a detector or only the pointer to the branch!!!
681 FillInfosFitPH(idect,nentries);
683 if (fDebugLevel != 1) {
688 if (fNumberFit > 0) {
689 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));
690 fStatisticMean = fStatisticMean / fNumberFit;
693 AliInfo(Form("There are %d with at least one entries. There is no fit!",fNumberEnt));
695 delete fDebugStreamer;
696 fDebugStreamer = 0x0;
699 //____________Functions fit Online PH2d________________________________________
700 Bool_t AliTRDCalibraFit::AnalysePH(AliTRDCalibraVector *calvect)
703 // Reconstruct the average pulse height from the vectorPH for each
705 // Reconstruct a drift velocity
706 // A first calibration of T0 is also made using the same method (slope method)
709 // Set the calibration mode
710 //const char *name = calvect->GetNamePH();
711 TString name = calvect->GetNamePH();
712 if(!SetModeCalibration(name,1)) return kFALSE;
714 // Number of Xbins (detectors or groups of pads)
715 if (!InitFit((432*calvect->GetDetCha0(1)+108*calvect->GetDetCha2(1)),1)) {
721 fStatisticMean = 0.0;
723 fNumberFitSuccess = 0;
725 // Init fCountDet and fCount
726 InitfCountDetAndfCount(1);
727 // Beginning of the loop
728 for (Int_t idect = fDect1; idect < fDect2; idect++) {
729 // Determination of fNnZ, fNnRphi, fNfragZ and fNfragRphi...........
730 UpdatefCountDetAndfCount(idect,1);
731 ReconstructFitRowMinRowMax(idect,1);
734 if(!calvect->GetPHEntries(fCountDet)) {
735 NotEnoughStatisticPH(idect,fEntriesCurrent);
740 TH1F *projph = calvect->CorrectTheError((TGraphErrors *) (calvect->ConvertVectorPHTGraphErrors(fCountDet,(idect-(fCount-(fCalibraMode->GetNfragZ(1)*fCalibraMode->GetNfragRphi(1)))),(const char *) tname)),fEntriesCurrent);
741 projph->SetDirectory(0);
742 if(fEntriesCurrent > 0) fNumberEnt++;
743 //printf("The number of entries for the group %d is %d\n",idect,fEntriesCurrent);
744 // This detector has not enough statistics or was off
745 if (fEntriesCurrent <= fMinEntries) {
746 //printf("Not enough stat!\n");
747 NotEnoughStatisticPH(idect,fEntriesCurrent);
750 // Statistic of the histos fitted
752 fStatisticMean += fEntriesCurrent;
753 // Calcul of "real" coef
754 CalculVdriftCoefMean();
759 case 0: FitLagrangePoly((TH1 *) projph); break;
760 case 1: FitPente((TH1 *) projph); break;
761 case 2: FitPH((TH1 *) projph,(Int_t) (idect - fDect1)); break;
762 default: return kFALSE;
764 // Fill the tree if end of a detector or only the pointer to the branch!!!
765 FillInfosFitPH(idect,fEntriesCurrent);
769 if (fNumberFit > 0) {
770 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));
771 fStatisticMean = fStatisticMean / fNumberFit;
774 AliInfo(Form("There are %d with at least one entries. There is no fit!",fNumberEnt));
776 delete fDebugStreamer;
777 fDebugStreamer = 0x0;
780 //____________Functions fit Online PRF2d_______________________________________
781 Bool_t AliTRDCalibraFit::AnalysePRF(const TProfile2D *prf)
784 // Take the 1D profiles (pad response function), projections of the 2D PRF
785 // on the Xaxis, for each calibration group
786 // Fit with a gaussian to reconstruct the sigma of the pad response function
789 // Set the calibration mode
790 //const char *name = prf->GetTitle();
791 TString name = prf->GetTitle();
792 if(!SetModeCalibration(name,2)) return kFALSE;
794 // Number of Ybins (detectors or groups of pads)
795 Int_t nybins = prf->GetNbinsY();// calibration groups
796 Int_t nbins = prf->GetNbinsX();// bins
797 Int_t nbg = GetNumberOfGroupsPRF((const char *)prf->GetTitle());
798 if((nbg > 0) || (nbg == -1)) return kFALSE;
799 if (!InitFit(nybins,2)) {
805 fStatisticMean = 0.0;
807 fNumberFitSuccess = 0;
809 // Init fCountDet and fCount
810 InitfCountDetAndfCount(2);
811 // Beginning of the loop
812 for (Int_t idect = fDect1; idect < fDect2; idect++) {
813 // Determination of fNnZ, fNnRphi, fNfragZ and fNfragRphi......
814 UpdatefCountDetAndfCount(idect,2);
815 ReconstructFitRowMinRowMax(idect,2);
817 TH1D *projprf = (TH1D *) prf->ProjectionX("projprf",idect+1,idect+1,(Option_t *) "e");
818 projprf->SetDirectory(0);
819 // Number of entries for this calibration group
820 Double_t nentries = 0;
821 for (Int_t k = 0; k < nbins; k++) {
822 Int_t binnb = (nbins+2)*(idect+1)+(k+1);
823 nentries += prf->GetBinEntries(binnb);
825 if(nentries > 0) fNumberEnt++;
826 // This detector has not enough statistics or was off
827 if (nentries <= fMinEntries) {
828 NotEnoughStatisticPRF(idect);
829 if (fDebugLevel != 1) {
834 // Statistics of the histos fitted
836 fStatisticMean += nentries;
837 // Calcul of "real" coef
842 case 0: FitPRF((TH1 *) projprf); break;
843 case 1: RmsPRF((TH1 *) projprf); break;
844 default: return kFALSE;
846 // Fill the tree if end of a detector or only the pointer to the branch!!!
847 FillInfosFitPRF(idect);
849 if (fDebugLevel != 1) {
854 if (fNumberFit > 0) {
855 AliInfo(Form("There are %d with at least one entries.",fNumberEnt));
856 AliInfo(Form("%d fits have been proceeded (sucessfully or not...).",fNumberFit));
857 AliInfo(Form("There is a mean statistic of: %d over these fitted histograms and %d successfulled fits"
858 ,(Int_t) fStatisticMean/fNumberFit,fNumberFitSuccess));
859 fStatisticMean = fStatisticMean / fNumberFit;
862 AliInfo(Form("There are %d with at least one entries. There is no fit!",fNumberEnt));
864 delete fDebugStreamer;
865 fDebugStreamer = 0x0;
868 //____________Functions fit Online PRF2d_______________________________________
869 Bool_t AliTRDCalibraFit::AnalysePRFMarianFit(const TProfile2D *prf)
872 // Take the 1D profiles (pad response function), projections of the 2D PRF
873 // on the Xaxis, for each calibration group
874 // Fit with a gaussian to reconstruct the sigma of the pad response function
877 // Set the calibration mode
878 //const char *name = prf->GetTitle();
879 TString name = prf->GetTitle();
880 if(!SetModeCalibration(name,2)) return kFALSE;
882 // Number of Ybins (detectors or groups of pads)
883 TAxis *xprf = prf->GetXaxis();
884 TAxis *yprf = prf->GetYaxis();
885 Int_t nybins = yprf->GetNbins();// calibration groups
886 Int_t nbins = xprf->GetNbins();// bins
887 Float_t lowedge = (Float_t) xprf->GetBinLowEdge(1);//lowedge in bins
888 Float_t upedge = (Float_t) xprf->GetBinUpEdge(nbins);//upedge in bins
889 Int_t nbg = GetNumberOfGroupsPRF((const char *)name);
890 if(nbg == -1) return kFALSE;
891 if(nbg > 0) fMethod = 1;
893 if (!InitFit(nybins,2)) {
899 fStatisticMean = 0.0;
901 fNumberFitSuccess = 0;
903 // Init fCountDet and fCount
904 InitfCountDetAndfCount(2);
905 // Beginning of the loop
906 for (Int_t idect = fDect1; idect < fDect2; idect++) {
907 // Determination of fNnZ, fNnRphi, fNfragZ and fNfragRphi.......
908 UpdatefCountDetAndfCount(idect,2);
909 ReconstructFitRowMinRowMax(idect,2);
910 // Build the array of entries and sum
911 TArrayD arraye = TArrayD(nbins);
912 TArrayD arraym = TArrayD(nbins);
913 TArrayD arrayme = TArrayD(nbins);
914 Double_t nentries = 0;
915 //printf("nbins %d\n",nbins);
916 for (Int_t k = 0; k < nbins; k++) {
917 Int_t binnb = (nbins+2)*(idect+1)+(k+1);
918 Double_t entries = (Double_t)prf->GetBinEntries(binnb);
919 Double_t mean = (Double_t)prf->GetBinContent(binnb);
920 Double_t error = (Double_t)prf->GetBinError(binnb);
921 //printf("for %d we have %f\n",k,entries);
923 arraye.AddAt(entries,k);
924 arraym.AddAt(mean,k);
925 arrayme.AddAt(error,k);
927 if(nentries > 0) fNumberEnt++;
928 //printf("The number of entries for the group %d is %f\n",idect,nentries);
929 // This detector has not enough statistics or was off
930 if (nentries <= fMinEntries) {
931 NotEnoughStatisticPRF(idect);
934 // Statistics of the histos fitted
936 fStatisticMean += nentries;
937 // Calcul of "real" coef
942 case 0: FitPRFGausMI( arraye.GetArray(), arraym.GetArray(), arrayme.GetArray(), nbins, lowedge, upedge); break;
943 case 1: FitTnpRange( arraye.GetArray(), arraym.GetArray(), arrayme.GetArray(), nbg, nbins); break;
944 default: return kFALSE;
946 // Fill the tree if end of a detector or only the pointer to the branch!!!
947 FillInfosFitPRF(idect);
950 if (fNumberFit > 0) {
951 AliInfo(Form("There are %d with at least one entries.",fNumberEnt));
952 AliInfo(Form("%d fits have been proceeded (sucessfully or not...).",fNumberFit));
953 AliInfo(Form("There is a mean statistic of: %d over these fitted histograms and %d successfulled fits"
954 ,(Int_t) fStatisticMean/fNumberFit,fNumberFitSuccess));
955 fStatisticMean = fStatisticMean / fNumberFit;
958 AliInfo(Form("There are %d with at least one entries. There is no fit!",fNumberEnt));
960 delete fDebugStreamer;
961 fDebugStreamer = 0x0;
964 //____________Functions fit Online PRF2d_______________________________________
965 Bool_t AliTRDCalibraFit::AnalysePRF(AliTRDCalibraVector *calvect)
968 // Reconstruct the 1D histo (pad response function) from the vectorPRD for
969 // each calibration group
970 // Fit with a gaussian to reconstruct the sigma of the pad response function
973 // Set the calibra mode
974 //const char *name = calvect->GetNamePRF();
975 TString name = calvect->GetNamePRF();
976 if(!SetModeCalibration(name,2)) return kFALSE;
977 //printf("test0 %s\n",name);
979 // Number of Xbins (detectors or groups of pads)
980 if (!InitFit((432*calvect->GetDetCha0(2)+108*calvect->GetDetCha2(2)),2)) {
985 ///printf("test2\n");
988 fStatisticMean = 0.0;
990 fNumberFitSuccess = 0;
992 // Init fCountDet and fCount
993 InitfCountDetAndfCount(2);
994 // Beginning of the loop
995 for (Int_t idect = fDect1; idect < fDect2; idect++) {
996 // Determination of fNnZ, fNnRphi, fNfragZ and fNfragRphi........
997 UpdatefCountDetAndfCount(idect,2);
998 ReconstructFitRowMinRowMax(idect,2);
1000 fEntriesCurrent = 0;
1001 if(!calvect->GetPRFEntries(fCountDet)) {
1002 NotEnoughStatisticPRF(idect);
1005 TString tname("PRF");
1007 TH1F *projprf = calvect->CorrectTheError((TGraphErrors *) (calvect->ConvertVectorPRFTGraphErrors(fCountDet,(idect-(fCount-(fCalibraMode->GetNfragZ(1)*fCalibraMode->GetNfragRphi(1)))),(const char *) tname)),fEntriesCurrent);
1008 projprf->SetDirectory(0);
1009 if(fEntriesCurrent > 0) fNumberEnt++;
1010 // This detector has not enough statistics or was off
1011 if (fEntriesCurrent <= fMinEntries) {
1012 NotEnoughStatisticPRF(idect);
1015 // Statistic of the histos fitted
1017 fStatisticMean += fEntriesCurrent;
1018 // Calcul of "real" coef
1019 CalculPRFCoefMean();
1023 case 1: FitPRF((TH1 *) projprf); break;
1024 case 2: RmsPRF((TH1 *) projprf); break;
1025 default: return kFALSE;
1027 // Fill the tree if end of a detector or only the pointer to the branch!!!
1028 FillInfosFitPRF(idect);
1031 if (fNumberFit > 0) {
1032 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));
1035 AliInfo(Form("There are %d with at least one entries. There is no fit!",fNumberEnt));
1037 delete fDebugStreamer;
1038 fDebugStreamer = 0x0;
1041 //____________Functions fit Online PRF2d_______________________________________
1042 Bool_t AliTRDCalibraFit::AnalysePRFMarianFit(AliTRDCalibraVector *calvect)
1045 // Reconstruct the 1D histo (pad response function) from the vectorPRD for
1046 // each calibration group
1047 // Fit with a gaussian to reconstruct the sigma of the pad response function
1050 // Set the calibra mode
1051 //const char *name = calvect->GetNamePRF();
1052 TString name = calvect->GetNamePRF();
1053 if(!SetModeCalibration(name,2)) return kFALSE;
1054 //printf("test0 %s\n",name);
1055 Int_t nbg = GetNumberOfGroupsPRF((const char *)name);
1056 //printf("test1 %d\n",nbg);
1057 if(nbg == -1) return kFALSE;
1058 if(nbg > 0) fMethod = 1;
1060 // Number of Xbins (detectors or groups of pads)
1061 if (!InitFit((432*calvect->GetDetCha0(2)+108*calvect->GetDetCha2(2)),2)) {
1062 //printf("test2\n");
1065 if (!InitFitPRF()) {
1066 //printf("test3\n");
1069 fStatisticMean = 0.0;
1071 fNumberFitSuccess = 0;
1075 Double_t *arrayx = 0;
1076 Double_t *arraye = 0;
1077 Double_t *arraym = 0;
1078 Double_t *arrayme = 0;
1079 Float_t lowedge = 0.0;
1080 Float_t upedge = 0.0;
1081 // Init fCountDet and fCount
1082 InitfCountDetAndfCount(2);
1083 // Beginning of the loop
1084 for (Int_t idect = fDect1; idect < fDect2; idect++) {
1085 // Determination of fNnZ, fNnRphi, fNfragZ and fNfragRphi......
1086 UpdatefCountDetAndfCount(idect,2);
1087 ReconstructFitRowMinRowMax(idect,2);
1089 fEntriesCurrent = 0;
1090 if(!calvect->GetPRFEntries(fCountDet)) {
1091 NotEnoughStatisticPRF(idect);
1094 TString tname("PRF");
1096 TGraphErrors *projprftree = calvect->ConvertVectorPRFTGraphErrors(fCountDet,(idect-(fCount-(fCalibraMode->GetNfragZ(1)*fCalibraMode->GetNfragRphi(1)))),(const char *) tname);
1097 nbins = projprftree->GetN();
1098 arrayx = (Double_t *)projprftree->GetX();
1099 arraye = (Double_t *)projprftree->GetEX();
1100 arraym = (Double_t *)projprftree->GetY();
1101 arrayme = (Double_t *)projprftree->GetEY();
1102 Float_t step = arrayx[1]-arrayx[0];
1103 lowedge = arrayx[0] - step/2.0;
1104 upedge = arrayx[(nbins-1)] + step/2.0;
1105 //printf("nbins est %d\n",nbins);
1106 for(Int_t k = 0; k < nbins; k++){
1107 fEntriesCurrent += (Int_t)arraye[k];
1108 //printf("for %d we have %f, %f\n",k,arraye[k],((projprftree->GetEX())[k]));
1109 if(arraye[k]>0.0) arrayme[k] = TMath::Sqrt(TMath::Abs(arrayme[k]-arraym[k]*arraym[k])/arraye[k]);
1111 if(fEntriesCurrent > 0) fNumberEnt++;
1112 //printf("The number of entries for the group %d is %d\n",idect,fEntriesCurrent);
1113 // This detector has not enough statistics or was off
1114 if (fEntriesCurrent <= fMinEntries) {
1115 NotEnoughStatisticPRF(idect);
1118 // Statistic of the histos fitted
1120 fStatisticMean += fEntriesCurrent;
1121 // Calcul of "real" coef
1122 CalculPRFCoefMean();
1126 case 0: FitPRFGausMI(arraye,arraym,arrayme,nbins,lowedge,upedge); break;
1127 case 1: FitTnpRange(arraye,arraym,arrayme,nbg,nbins); break;
1128 default: return kFALSE;
1130 // Fill the tree if end of a detector or only the pointer to the branch!!!
1131 FillInfosFitPRF(idect);
1134 if (fNumberFit > 0) {
1135 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));
1138 AliInfo(Form("There are %d with at least one entries. There is no fit!",fNumberEnt));
1140 delete fDebugStreamer;
1141 fDebugStreamer = 0x0;
1144 //____________Functions fit Online CH2d________________________________________
1145 Bool_t AliTRDCalibraFit::AnalyseLinearFitters(AliTRDCalibraVdriftLinearFit *calivdli)
1148 // The linear method
1151 fStatisticMean = 0.0;
1153 fNumberFitSuccess = 0;
1155 if(!InitFitLinearFitter()) return kFALSE;
1158 for(Int_t idet = 0; idet < 540; idet++){
1161 //printf("detector number %d\n",idet);
1166 fEntriesCurrent = 0;
1168 Bool_t here = calivdli->GetParam(idet,¶m);
1169 Bool_t heree = calivdli->GetError(idet,&error);
1170 //printf("here %d and heree %d\n",here, heree);
1172 fEntriesCurrent = (Int_t) error[2];
1175 //printf("Number of entries %d\n",fEntriesCurrent);
1176 // Nothing found or not enough statistic
1177 if((!heree) || (!here) || (fEntriesCurrent <= fMinEntries)) {
1178 NotEnoughStatisticLinearFitter();
1185 fStatisticMean += fEntriesCurrent;
1188 if((-(param[1])) <= 0.0) {
1189 NotEnoughStatisticLinearFitter();
1193 // CalculDatabaseVdriftandTan
1194 CalculVdriftLorentzCoef();
1197 fNumberFitSuccess ++;
1199 // Put the fCurrentCoef
1200 fCurrentCoef[0] = -param[1];
1201 // here the database must be the one of the reconstruction for the lorentz angle....
1202 fCurrentCoef2[0] = (param[0]+fCurrentCoef[1]*fCurrentCoef2[1])/fCurrentCoef[0];
1203 fCurrentCoefE = error[1];
1204 fCurrentCoefE2 = error[0];
1205 if((TMath::Abs(fCurrentCoef2[0]) > 0.0000001) && (TMath::Abs(param[0]) > 0.0000001)){
1206 fCurrentCoefE2 = (fCurrentCoefE2/param[0]+fCurrentCoefE/fCurrentCoef[0])*fCurrentCoef2[0];
1210 FillInfosFitLinearFitter();
1215 if (fNumberFit > 0) {
1216 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));
1219 AliInfo(Form("There are %d with at least one entries. There is no fit!",fNumberEnt));
1221 delete fDebugStreamer;
1222 fDebugStreamer = 0x0;
1226 //____________Functions fit Online CH2d________________________________________
1227 Double_t AliTRDCalibraFit::AnalyseLinearFittersAllTogether(AliTRDCalibraVdriftLinearFit *calivdli)
1230 // The linear method
1235 TH2S *linearfitterhisto = 0x0;
1237 for(Int_t idet = 0; idet < 540; idet++){
1239 TH2S * u = calivdli->GetLinearFitterHistoForce(idet);
1240 if(idet == 0) linearfitterhisto = u;
1241 else linearfitterhisto->Add(u);
1248 TAxis *xaxis = linearfitterhisto->GetXaxis();
1249 TAxis *yaxis = linearfitterhisto->GetYaxis();
1250 TLinearFitter linearfitter = TLinearFitter(2,"pol1");
1252 Double_t integral = linearfitterhisto->Integral();
1253 //printf("Integral is %f\n",integral);
1254 Bool_t securitybreaking = kFALSE;
1255 if(TMath::Abs(integral-1199) < 0.00001) securitybreaking = kTRUE;
1256 for(Int_t ibinx = 0; ibinx < linearfitterhisto->GetNbinsX(); ibinx++){
1257 for(Int_t ibiny = 0; ibiny < linearfitterhisto->GetNbinsY(); ibiny++){
1258 if(linearfitterhisto->GetBinContent(ibinx+1,ibiny+1)>0){
1259 Double_t x = xaxis->GetBinCenter(ibinx+1);
1260 Double_t y = yaxis->GetBinCenter(ibiny+1);
1262 for(Int_t k = 0; k < (Int_t)linearfitterhisto->GetBinContent(ibinx+1,ibiny+1); k++){
1263 if(!securitybreaking){
1264 linearfitter.AddPoint(&x,y);
1269 linearfitter.AddPoint(&x,y);
1279 //printf("Find %d entries for the detector %d\n",arrayI[cb],cb);
1281 // Eval the linear fitter
1282 if(entries > fMinEntries){
1283 TVectorD par = TVectorD(2);
1285 if((linearfitter.EvalRobust(0.8)==0)) {
1286 //printf("Take the param\n");
1287 linearfitter.GetParameters(par);
1290 //printf("Finish\n");
1291 // Put the fCurrentCoef
1292 fCurrentCoef[0] = -par[1];
1293 // here the database must be the one of the reconstruction for the lorentz angle....
1294 fCurrentCoef2[0] = (par[0]+fCurrentCoef[1]*fCurrentCoef2[1])/fCurrentCoef[0];
1296 return fCurrentCoef[0];
1306 delete linearfitterhisto;
1307 delete fDebugStreamer;
1308 fDebugStreamer = 0x0;
1311 //____________Functions for seeing if the pad is really okey___________________
1312 //_____________________________________________________________________________
1313 Int_t AliTRDCalibraFit::GetNumberOfGroupsPRF(TString nametitle)
1316 // Get numberofgroupsprf
1320 const Char_t *pattern0 = "Ngp0";
1321 const Char_t *pattern1 = "Ngp1";
1322 const Char_t *pattern2 = "Ngp2";
1323 const Char_t *pattern3 = "Ngp3";
1324 const Char_t *pattern4 = "Ngp4";
1325 const Char_t *pattern5 = "Ngp5";
1326 const Char_t *pattern6 = "Ngp6";
1329 if (strstr(nametitle.Data(),pattern0)) {
1332 if (strstr(nametitle.Data(),pattern1)) {
1335 if (strstr(nametitle.Data(),pattern2)) {
1338 if (strstr(nametitle.Data(),pattern3)) {
1341 if (strstr(nametitle.Data(),pattern4)) {
1344 if (strstr(nametitle.Data(),pattern5)) {
1347 if (strstr(nametitle.Data(),pattern6)){
1354 //_____________________________________________________________________________
1355 Bool_t AliTRDCalibraFit::SetModeCalibration(TString name, Int_t i)
1358 // Set fNz[i] and fNrphi[i] of the AliTRDCalibraFit::Instance()
1359 // corresponding to the given name
1362 if(!SetNzFromTObject(name,i)) return kFALSE;
1363 if(!SetNrphiFromTObject(name,i)) return kFALSE;
1368 //_____________________________________________________________________________
1369 Bool_t AliTRDCalibraFit::SetNrphiFromTObject(TString name, Int_t i)
1372 // Set fNrphi[i] of the AliTRDCalibraFit::Instance()
1373 // corresponding to the given TObject
1377 const Char_t *patternrphi0 = "Nrphi0";
1378 const Char_t *patternrphi1 = "Nrphi1";
1379 const Char_t *patternrphi2 = "Nrphi2";
1380 const Char_t *patternrphi3 = "Nrphi3";
1381 const Char_t *patternrphi4 = "Nrphi4";
1382 const Char_t *patternrphi5 = "Nrphi5";
1383 const Char_t *patternrphi6 = "Nrphi6";
1386 const Char_t *patternrphi10 = "Nrphi10";
1387 const Char_t *patternrphi100 = "Nrphi100";
1388 const Char_t *patternz10 = "Nz10";
1389 const Char_t *patternz100 = "Nz100";
1392 if ((strstr(name.Data(),patternrphi100)) && (strstr(name.Data(),patternz100))) {
1393 fCalibraMode->SetAllTogether(i);
1395 if (fDebugLevel > 1) {
1396 AliInfo(Form("fNbDet %d and 100",fNbDet));
1400 if ((strstr(name.Data(),patternrphi10)) && (strstr(name.Data(),patternz10))) {
1401 fCalibraMode->SetPerSuperModule(i);
1403 if (fDebugLevel > 1) {
1404 AliInfo(Form("fNDet %d and 100",fNbDet));
1409 if (strstr(name.Data(),patternrphi0)) {
1410 fCalibraMode->SetNrphi(i ,0);
1411 if (fDebugLevel > 1) {
1412 AliInfo(Form("fNbDet %d and 0",fNbDet));
1416 if (strstr(name.Data(),patternrphi1)) {
1417 fCalibraMode->SetNrphi(i, 1);
1418 if (fDebugLevel > 1) {
1419 AliInfo(Form("fNbDet %d and 1",fNbDet));
1423 if (strstr(name.Data(),patternrphi2)) {
1424 fCalibraMode->SetNrphi(i, 2);
1425 if (fDebugLevel > 1) {
1426 AliInfo(Form("fNbDet %d and 2",fNbDet));
1430 if (strstr(name.Data(),patternrphi3)) {
1431 fCalibraMode->SetNrphi(i, 3);
1432 if (fDebugLevel > 1) {
1433 AliInfo(Form("fNbDet %d and 3",fNbDet));
1437 if (strstr(name.Data(),patternrphi4)) {
1438 fCalibraMode->SetNrphi(i, 4);
1439 if (fDebugLevel > 1) {
1440 AliInfo(Form("fNbDet %d and 4",fNbDet));
1444 if (strstr(name.Data(),patternrphi5)) {
1445 fCalibraMode->SetNrphi(i, 5);
1446 if (fDebugLevel > 1) {
1447 AliInfo(Form("fNbDet %d and 5",fNbDet));
1451 if (strstr(name.Data(),patternrphi6)) {
1452 fCalibraMode->SetNrphi(i, 6);
1453 if (fDebugLevel > 1) {
1454 AliInfo(Form("fNbDet %d and 6",fNbDet));
1459 if (fDebugLevel > 1) {
1460 AliInfo(Form("fNbDet %d and rest",fNbDet));
1462 fCalibraMode->SetNrphi(i ,0);
1466 //_____________________________________________________________________________
1467 Bool_t AliTRDCalibraFit::SetNzFromTObject(TString name, Int_t i)
1470 // Set fNz[i] of the AliTRDCalibraFit::Instance()
1471 // corresponding to the given TObject
1475 const Char_t *patternz0 = "Nz0";
1476 const Char_t *patternz1 = "Nz1";
1477 const Char_t *patternz2 = "Nz2";
1478 const Char_t *patternz3 = "Nz3";
1479 const Char_t *patternz4 = "Nz4";
1481 const Char_t *patternrphi10 = "Nrphi10";
1482 const Char_t *patternrphi100 = "Nrphi100";
1483 const Char_t *patternz10 = "Nz10";
1484 const Char_t *patternz100 = "Nz100";
1486 if ((strstr(name.Data(),patternrphi100)) && (strstr(name.Data(),patternz100))) {
1487 fCalibraMode->SetAllTogether(i);
1489 if (fDebugLevel > 1) {
1490 AliInfo(Form("fNbDet %d and 100",fNbDet));
1494 if ((strstr(name.Data(),patternrphi10)) && (strstr(name.Data(),patternz10))) {
1495 fCalibraMode->SetPerSuperModule(i);
1497 if (fDebugLevel > 1) {
1498 AliInfo(Form("fNbDet %d and 10",fNbDet));
1502 if (strstr(name.Data(),patternz0)) {
1503 fCalibraMode->SetNz(i, 0);
1504 if (fDebugLevel > 1) {
1505 AliInfo(Form("fNbDet %d and 0",fNbDet));
1509 if (strstr(name.Data(),patternz1)) {
1510 fCalibraMode->SetNz(i ,1);
1511 if (fDebugLevel > 1) {
1512 AliInfo(Form("fNbDet %d and 1",fNbDet));
1516 if (strstr(name.Data(),patternz2)) {
1517 fCalibraMode->SetNz(i ,2);
1518 if (fDebugLevel > 1) {
1519 AliInfo(Form("fNbDet %d and 2",fNbDet));
1523 if (strstr(name.Data(),patternz3)) {
1524 fCalibraMode->SetNz(i ,3);
1525 if (fDebugLevel > 1) {
1526 AliInfo(Form("fNbDet %d and 3",fNbDet));
1530 if (strstr(name.Data(),patternz4)) {
1531 fCalibraMode->SetNz(i ,4);
1532 if (fDebugLevel > 1) {
1533 AliInfo(Form("fNbDet %d and 4",fNbDet));
1538 if (fDebugLevel > 1) {
1539 AliInfo(Form("fNbDet %d and rest",fNbDet));
1541 fCalibraMode->SetNz(i ,0);
1544 //______________________________________________________________________
1545 void AliTRDCalibraFit::RemoveOutliers(Int_t type, Bool_t perdetector){
1547 // Remove the results too far from the mean value and rms
1548 // type: 0 gain, 1 vdrift
1552 Int_t loop = (Int_t) fVectorFit.GetEntriesFast();
1554 AliInfo("The Vector Fit is not complete!");
1557 Int_t detector = -1;
1559 Float_t value = 0.0;
1561 /////////////////////////////////
1562 // Calculate the mean values
1563 ////////////////////////////////
1565 ////////////////////////
1566 Double_t meanAll = 0.0;
1567 Double_t rmsAll = 0.0;
1572 for (Int_t k = 0; k < loop; k++) {
1573 detector = ((AliTRDFitInfo *) fVectorFit.At(k))->GetDetector();
1574 sector = GetSector(detector);
1576 value = ((AliTRDFitInfo *) fVectorFit.At(k))->GetCoef()[0];
1578 rmsAll += value*value;
1584 Int_t rowMax = fGeo->GetRowMax(GetLayer(detector),GetStack(detector),GetSector(detector));
1585 Int_t colMax = fGeo->GetColMax(GetLayer(detector));
1586 for (Int_t row = 0; row < rowMax; row++) {
1587 for (Int_t col = 0; col < colMax; col++) {
1588 value = ((AliTRDFitInfo *) fVectorFit.At(k))->GetCoef()[(Int_t)(col*rowMax+row)];
1590 rmsAll += value*value;
1600 meanAll = meanAll/countAll;
1601 rmsAll = TMath::Sqrt(TMath::Abs(rmsAll/countAll - (meanAll*meanAll)));
1603 //printf("RemoveOutliers: meanAll %f and rmsAll %f\n",meanAll,rmsAll);
1604 /////////////////////////////////////////////////
1606 ////////////////////////////////////////////////
1607 Double_t defaultvalue = -1.0;
1608 if(type==1) defaultvalue = -1.5;
1609 for (Int_t k = 0; k < loop; k++) {
1610 detector = ((AliTRDFitInfo *) fVectorFit.At(k))->GetDetector();
1611 sector = GetSector(detector);
1612 Int_t rowMax = fGeo->GetRowMax(GetLayer(detector),GetStack(detector),GetSector(detector));
1613 Int_t colMax = fGeo->GetColMax(GetLayer(detector));
1614 Float_t *coef = ((AliTRDFitInfo *) fVectorFit.At(k))->GetCoef();
1616 // remove the results too far away
1617 for (Int_t row = 0; row < rowMax; row++) {
1618 for (Int_t col = 0; col < colMax; col++) {
1619 value = coef[(Int_t)(col*rowMax+row)];
1620 if((value > 0.0) && (rmsAll > 0.0) && (TMath::Abs(value-meanAll) > (2*rmsAll))) {
1621 coef[(Int_t)(col*rowMax+row)] = defaultvalue;
1627 //______________________________________________________________________
1628 void AliTRDCalibraFit::RemoveOutliers2(Bool_t perdetector){
1630 // Remove the results too far from the mean and rms
1634 Int_t loop = (Int_t) fVectorFit2.GetEntriesFast();
1636 AliInfo("The Vector Fit is not complete!");
1639 Int_t detector = -1;
1641 Float_t value = 0.0;
1643 /////////////////////////////////
1644 // Calculate the mean values
1645 ////////////////////////////////
1647 ////////////////////////
1648 Double_t meanAll = 0.0;
1649 Double_t rmsAll = 0.0;
1654 for (Int_t k = 0; k < loop; k++) {
1655 detector = ((AliTRDFitInfo *) fVectorFit2.At(k))->GetDetector();
1656 sector = GetSector(detector);
1658 value = ((AliTRDFitInfo *) fVectorFit2.At(k))->GetCoef()[0];
1661 rmsAll += value*value;
1666 Int_t rowMax = fGeo->GetRowMax(GetLayer(detector),GetStack(detector),GetSector(detector));
1667 Int_t colMax = fGeo->GetColMax(GetLayer(detector));
1668 for (Int_t row = 0; row < rowMax; row++) {
1669 for (Int_t col = 0; col < colMax; col++) {
1670 value = ((AliTRDFitInfo *) fVectorFit2.At(k))->GetCoef()[(Int_t)(col*rowMax+row)];
1672 rmsAll += value*value;
1681 meanAll = meanAll/countAll;
1682 rmsAll = TMath::Sqrt(TMath::Abs(rmsAll/countAll - (meanAll*meanAll)));
1684 //printf("Remove outliers 2: meanAll %f, rmsAll %f\n",meanAll,rmsAll);
1685 /////////////////////////////////////////////////
1687 ////////////////////////////////////////////////
1688 for (Int_t k = 0; k < loop; k++) {
1689 detector = ((AliTRDFitInfo *) fVectorFit2.At(k))->GetDetector();
1690 sector = GetSector(detector);
1691 Int_t rowMax = fGeo->GetRowMax(GetLayer(detector),GetStack(detector),GetSector(detector));
1692 Int_t colMax = fGeo->GetColMax(GetLayer(detector));
1693 Float_t *coef = ((AliTRDFitInfo *) fVectorFit2.At(k))->GetCoef();
1695 // remove the results too far away
1696 for (Int_t row = 0; row < rowMax; row++) {
1697 for (Int_t col = 0; col < colMax; col++) {
1698 value = coef[(Int_t)(col*rowMax+row)];
1699 if((value < 70.0) && (rmsAll > 0.0) && (TMath::Abs(value-meanAll) > (2.5*rmsAll))) coef[(Int_t)(col*rowMax+row)] = 100.0;
1704 //______________________________________________________________________
1705 void AliTRDCalibraFit::PutMeanValueOtherVectorFit(Int_t ofwhat, Bool_t perdetector){
1707 // ofwhat is equaled to 0: mean value of all passing detectors
1708 // ofwhat is equaled to 1: mean value of the detector, otherwise supermodule, otherwise all
1711 Int_t loop = (Int_t) fVectorFit.GetEntriesFast();
1713 AliInfo("The Vector Fit is not complete!");
1716 Int_t detector = -1;
1718 Float_t value = 0.0;
1720 /////////////////////////////////
1721 // Calculate the mean values
1722 ////////////////////////////////
1724 ////////////////////////
1725 Double_t meanAll = 0.0;
1726 Double_t meanSupermodule[18];
1727 Double_t meanDetector[540];
1728 Double_t rmsAll = 0.0;
1729 Double_t rmsSupermodule[18];
1730 Double_t rmsDetector[540];
1732 Int_t countSupermodule[18];
1733 Int_t countDetector[540];
1734 for(Int_t sm = 0; sm < 18; sm++){
1735 rmsSupermodule[sm] = 0.0;
1736 meanSupermodule[sm] = 0.0;
1737 countSupermodule[sm] = 0;
1739 for(Int_t det = 0; det < 540; det++){
1740 rmsDetector[det] = 0.0;
1741 meanDetector[det] = 0.0;
1742 countDetector[det] = 0;
1747 for (Int_t k = 0; k < loop; k++) {
1748 detector = ((AliTRDFitInfo *) fVectorFit.At(k))->GetDetector();
1749 sector = GetSector(detector);
1751 value = ((AliTRDFitInfo *) fVectorFit.At(k))->GetCoef()[0];
1753 rmsDetector[detector] += value*value;
1754 meanDetector[detector] += value;
1755 countDetector[detector]++;
1756 rmsSupermodule[sector] += value*value;
1757 meanSupermodule[sector] += value;
1758 countSupermodule[sector]++;
1759 rmsAll += value*value;
1765 Int_t rowMax = fGeo->GetRowMax(GetLayer(detector),GetStack(detector),GetSector(detector));
1766 Int_t colMax = fGeo->GetColMax(GetLayer(detector));
1767 for (Int_t row = 0; row < rowMax; row++) {
1768 for (Int_t col = 0; col < colMax; col++) {
1769 value = ((AliTRDFitInfo *) fVectorFit.At(k))->GetCoef()[(Int_t)(col*rowMax+row)];
1771 rmsDetector[detector] += value*value;
1772 meanDetector[detector] += value;
1773 countDetector[detector]++;
1774 rmsSupermodule[sector] += value*value;
1775 meanSupermodule[sector] += value;
1776 countSupermodule[sector]++;
1777 rmsAll += value*value;
1787 meanAll = meanAll/countAll;
1788 rmsAll = TMath::Abs(rmsAll/countAll - (meanAll*meanAll));
1790 for(Int_t sm = 0; sm < 18; sm++){
1791 if(countSupermodule[sm] > 0) {
1792 meanSupermodule[sm] = meanSupermodule[sm]/countSupermodule[sm];
1793 rmsSupermodule[sm] = TMath::Abs(rmsSupermodule[sm]/countSupermodule[sm] - (meanSupermodule[sm]*meanSupermodule[sm]));
1796 for(Int_t det = 0; det < 540; det++){
1797 if(countDetector[det] > 0) {
1798 meanDetector[det] = meanDetector[det]/countDetector[det];
1799 rmsDetector[det] = TMath::Abs(rmsDetector[det]/countDetector[det] - (meanDetector[det]*meanDetector[det]));
1802 //printf("Put mean value, meanAll %f, rmsAll %f\n",meanAll,rmsAll);
1803 ///////////////////////////////////////////////
1804 // Put the mean value for the no-fitted
1805 /////////////////////////////////////////////
1806 for (Int_t k = 0; k < loop; k++) {
1807 detector = ((AliTRDFitInfo *) fVectorFit.At(k))->GetDetector();
1808 sector = GetSector(detector);
1809 Int_t rowMax = fGeo->GetRowMax(GetLayer(detector),GetStack(detector),GetSector(detector));
1810 Int_t colMax = fGeo->GetColMax(GetLayer(detector));
1811 Float_t *coef = ((AliTRDFitInfo *) fVectorFit.At(k))->GetCoef();
1813 for (Int_t row = 0; row < rowMax; row++) {
1814 for (Int_t col = 0; col < colMax; col++) {
1815 value = coef[(Int_t)(col*rowMax+row)];
1817 if((ofwhat == 0) && (meanAll > 0.0) && (countAll > 15)) coef[(Int_t)(col*rowMax+row)] = -TMath::Abs(meanAll);
1819 if((meanDetector[detector] > 0.0) && (countDetector[detector] > 20)) coef[(Int_t)(col*rowMax+row)] = -TMath::Abs(meanDetector[detector]);
1820 else if((meanSupermodule[sector] > 0.0) && (countSupermodule[sector] > 15)) coef[(Int_t)(col*rowMax+row)] = -TMath::Abs(meanSupermodule[sector]);
1821 else if((meanAll > 0.0) && (countAll > 15)) coef[(Int_t)(col*rowMax+row)] = -TMath::Abs(meanAll);
1825 if(fDebugLevel > 1){
1827 if ( !fDebugStreamer ) {
1829 TDirectory *backup = gDirectory;
1830 fDebugStreamer = new TTreeSRedirector("TRDDebugFit.root");
1831 if ( backup ) backup->cd(); //we don't want to be cd'd to the debug streamer
1834 Float_t coefnow = coef[(Int_t)(col*rowMax+row)];
1836 (* fDebugStreamer) << "PutMeanValueOtherVectorFit"<<
1837 "detector="<<detector<<
1849 //______________________________________________________________________
1850 void AliTRDCalibraFit::PutMeanValueOtherVectorFit2(Int_t ofwhat, Bool_t perdetector){
1852 // ofwhat is equaled to 0: mean value of all passing detectors
1853 // ofwhat is equaled to 1: mean value of the detector, otherwise supermodule, otherwise all
1856 Int_t loop = (Int_t) fVectorFit2.GetEntriesFast();
1858 AliInfo("The Vector Fit is not complete!");
1861 Int_t detector = -1;
1863 Float_t value = 0.0;
1865 /////////////////////////////////
1866 // Calculate the mean values
1867 ////////////////////////////////
1869 ////////////////////////
1870 Double_t meanAll = 0.0;
1871 Double_t rmsAll = 0.0;
1872 Double_t meanSupermodule[18];
1873 Double_t rmsSupermodule[18];
1874 Double_t meanDetector[540];
1875 Double_t rmsDetector[540];
1877 Int_t countSupermodule[18];
1878 Int_t countDetector[540];
1879 for(Int_t sm = 0; sm < 18; sm++){
1880 rmsSupermodule[sm] = 0.0;
1881 meanSupermodule[sm] = 0.0;
1882 countSupermodule[sm] = 0;
1884 for(Int_t det = 0; det < 540; det++){
1885 rmsDetector[det] = 0.0;
1886 meanDetector[det] = 0.0;
1887 countDetector[det] = 0;
1891 for (Int_t k = 0; k < loop; k++) {
1892 detector = ((AliTRDFitInfo *) fVectorFit2.At(k))->GetDetector();
1893 sector = GetSector(detector);
1895 value = ((AliTRDFitInfo *) fVectorFit2.At(k))->GetCoef()[0];
1897 rmsDetector[detector] += value*value;
1898 meanDetector[detector] += value;
1899 countDetector[detector]++;
1900 rmsSupermodule[sector] += value*value;
1901 meanSupermodule[sector] += value;
1902 countSupermodule[sector]++;
1904 rmsAll += value*value;
1909 Int_t rowMax = fGeo->GetRowMax(GetLayer(detector),GetStack(detector),GetSector(detector));
1910 Int_t colMax = fGeo->GetColMax(GetLayer(detector));
1911 for (Int_t row = 0; row < rowMax; row++) {
1912 for (Int_t col = 0; col < colMax; col++) {
1913 value = ((AliTRDFitInfo *) fVectorFit2.At(k))->GetCoef()[(Int_t)(col*rowMax+row)];
1915 rmsDetector[detector] += value*value;
1916 meanDetector[detector] += value;
1917 countDetector[detector]++;
1918 rmsSupermodule[sector] += value*value;
1919 meanSupermodule[sector] += value;
1920 countSupermodule[sector]++;
1921 rmsAll += value*value;
1931 meanAll = meanAll/countAll;
1932 rmsAll = TMath::Abs(rmsAll/countAll - (meanAll*meanAll));
1934 for(Int_t sm = 0; sm < 18; sm++){
1935 if(countSupermodule[sm] > 0) {
1936 meanSupermodule[sm] = meanSupermodule[sm]/countSupermodule[sm];
1937 rmsSupermodule[sm] = TMath::Abs(rmsSupermodule[sm]/countSupermodule[sm] - (meanSupermodule[sm]*meanSupermodule[sm]));
1940 for(Int_t det = 0; det < 540; det++){
1941 if(countDetector[det] > 0) {
1942 meanDetector[det] = meanDetector[det]/countDetector[det];
1943 rmsDetector[det] = TMath::Abs(rmsDetector[det]/countDetector[det] - (meanDetector[det]*meanDetector[det]));
1946 //printf("Put mean value 2: meanAll %f, rmsAll %f\n",meanAll,rmsAll);
1947 ////////////////////////////////////////////
1948 // Put the mean value for the no-fitted
1949 /////////////////////////////////////////////
1950 for (Int_t k = 0; k < loop; k++) {
1951 detector = ((AliTRDFitInfo *) fVectorFit2.At(k))->GetDetector();
1952 sector = GetSector(detector);
1953 Int_t rowMax = fGeo->GetRowMax(GetLayer(detector),GetStack(detector),GetSector(detector));
1954 Int_t colMax = fGeo->GetColMax(GetLayer(detector));
1955 Float_t *coef = ((AliTRDFitInfo *) fVectorFit2.At(k))->GetCoef();
1957 for (Int_t row = 0; row < rowMax; row++) {
1958 for (Int_t col = 0; col < colMax; col++) {
1959 value = coef[(Int_t)(col*rowMax+row)];
1961 if((ofwhat == 0) && (meanAll > -1.5) && (countAll > 15)) coef[(Int_t)(col*rowMax+row)] = meanAll+100.0;
1963 if((meanDetector[detector] > -1.5) && (countDetector[detector] > 20)) coef[(Int_t)(col*rowMax+row)] = meanDetector[detector]+100.0;
1964 else if((meanSupermodule[sector] > -1.5) && (countSupermodule[sector] > 15)) coef[(Int_t)(col*rowMax+row)] = meanSupermodule[sector]+100.0;
1965 else if((meanAll > -1.5) && (countAll > 15)) coef[(Int_t)(col*rowMax+row)] = meanAll+100.0;
1969 if(fDebugLevel > 1){
1971 if ( !fDebugStreamer ) {
1973 TDirectory *backup = gDirectory;
1974 fDebugStreamer = new TTreeSRedirector("TRDDebugFit.root");
1975 if ( backup ) backup->cd(); //we don't want to be cd'd to the debug streamer
1978 Float_t coefnow = coef[(Int_t)(col*rowMax+row)];
1980 (* fDebugStreamer) << "PutMeanValueOtherVectorFit2"<<
1981 "detector="<<detector<<
1994 //_____________________________________________________________________________
1995 AliTRDCalDet *AliTRDCalibraFit::CreateDetObjectVdrift(const TObjArray *vectorFit, Bool_t perdetector)
1998 // It creates the AliTRDCalDet object from the AliTRDFitInfo
1999 // It takes the mean value of the coefficients per detector
2000 // This object has to be written in the database
2003 // Create the DetObject
2004 AliTRDCalDet *object = new AliTRDCalDet("ChamberVdrift","TRD drift velocities (detector value)");
2006 Int_t loop = (Int_t) vectorFit->GetEntriesFast();
2007 if(loop != 540) AliInfo("The Vector Fit is not complete!");
2008 Int_t detector = -1;
2009 Float_t value = 0.0;
2012 for (Int_t k = 0; k < loop; k++) {
2013 detector = ((AliTRDFitInfo *) vectorFit->At(k))->GetDetector();
2016 mean = TMath::Abs(((AliTRDFitInfo *) vectorFit->At(k))->GetCoef()[0]);
2020 Int_t rowMax = fGeo->GetRowMax(GetLayer(detector),GetStack(detector),GetSector(detector));
2021 Int_t colMax = fGeo->GetColMax(GetLayer(detector));
2022 for (Int_t row = 0; row < rowMax; row++) {
2023 for (Int_t col = 0; col < colMax; col++) {
2024 value = ((AliTRDFitInfo *) vectorFit->At(k))->GetCoef()[(Int_t)(col*rowMax+row)];
2025 mean += TMath::Abs(value);
2029 if(count > 0) mean = mean/count;
2031 object->SetValue(detector,mean);
2036 //_____________________________________________________________________________
2037 AliTRDCalDet *AliTRDCalibraFit::CreateDetObjectGain(const TObjArray *vectorFit, Bool_t meanOtherBefore, Double_t scaleFitFactor, Bool_t perdetector)
2040 // It creates the AliTRDCalDet object from the AliTRDFitInfo
2041 // It takes the mean value of the coefficients per detector
2042 // This object has to be written in the database
2045 // Create the DetObject
2046 AliTRDCalDet *object = new AliTRDCalDet("ChamberGainFactor","GainFactor (detector value)");
2048 fScaleGain = scaleFitFactor;
2050 Int_t loop = (Int_t) vectorFit->GetEntriesFast();
2051 if(loop != 540) AliInfo("The Vector Fit is not complete!");
2052 Int_t detector = -1;
2053 Float_t value = 0.0;
2055 for (Int_t k = 0; k < loop; k++) {
2056 detector = ((AliTRDFitInfo *) vectorFit->At(k))->GetDetector();
2059 value = ((AliTRDFitInfo *) vectorFit->At(k))->GetCoef()[0];
2060 if(!meanOtherBefore){
2061 if(value > 0) value = value*scaleFitFactor;
2063 else value = value*scaleFitFactor;
2064 mean = TMath::Abs(value);
2068 Int_t rowMax = fGeo->GetRowMax(GetLayer(detector),GetStack(detector),GetSector(detector));
2069 Int_t colMax = fGeo->GetColMax(GetLayer(detector));
2070 for (Int_t row = 0; row < rowMax; row++) {
2071 for (Int_t col = 0; col < colMax; col++) {
2072 value = ((AliTRDFitInfo *) vectorFit->At(k))->GetCoef()[(Int_t)(col*rowMax+row)];
2073 if(!meanOtherBefore) {
2074 if(value > 0) value = value*scaleFitFactor;
2076 else value = value*scaleFitFactor;
2077 mean += TMath::Abs(value);
2081 if(count > 0) mean = mean/count;
2083 if(mean < 0.1) mean = 0.1;
2084 object->SetValue(detector,mean);
2089 //_____________________________________________________________________________
2090 AliTRDCalDet *AliTRDCalibraFit::CreateDetObjectT0(const TObjArray *vectorFit, Bool_t perdetector)
2093 // It creates the AliTRDCalDet object from the AliTRDFitInfo2
2094 // It takes the min value of the coefficients per detector
2095 // This object has to be written in the database
2098 // Create the DetObject
2099 AliTRDCalDet *object = new AliTRDCalDet("ChamberT0","T0 (detector value)");
2101 Int_t loop = (Int_t) vectorFit->GetEntriesFast();
2102 if(loop != 540) AliInfo("The Vector Fit is not complete!");
2103 Int_t detector = -1;
2104 Float_t value = 0.0;
2106 for (Int_t k = 0; k < loop; k++) {
2107 detector = ((AliTRDFitInfo *) vectorFit->At(k))->GetDetector();
2108 Float_t min = 100.0;
2110 value = ((AliTRDFitInfo *) vectorFit->At(k))->GetCoef()[0];
2112 if(value > 70.0) value = value-100.0;
2117 Int_t rowMax = fGeo->GetRowMax(GetLayer(detector),GetStack(detector),GetSector(detector));
2118 Int_t colMax = fGeo->GetColMax(GetLayer(detector));
2119 for (Int_t row = 0; row < rowMax; row++) {
2120 for (Int_t col = 0; col < colMax; col++) {
2121 value = ((AliTRDFitInfo *) vectorFit->At(k))->GetCoef()[(Int_t)(col*rowMax+row)];
2123 if(value > 70.0) value = value-100.0;
2125 if(min > value) min = value;
2129 object->SetValue(detector,min);
2135 //_____________________________________________________________________________
2136 AliTRDCalDet *AliTRDCalibraFit::CreateDetObjectLorentzAngle(const TObjArray *vectorFit)
2139 // It creates the AliTRDCalDet object from the AliTRDFitInfo2
2140 // It takes the min value of the coefficients per detector
2141 // This object has to be written in the database
2144 // Create the DetObject
2145 AliTRDCalDet *object = new AliTRDCalDet("tan(lorentzangle)","tan(lorentzangle) (detector value)");
2148 Int_t loop = (Int_t) vectorFit->GetEntriesFast();
2149 if(loop != 540) AliInfo("The Vector Fit is not complete!");
2150 Int_t detector = -1;
2151 Float_t value = 0.0;
2153 for (Int_t k = 0; k < loop; k++) {
2154 detector = ((AliTRDFitInfo *) vectorFit->At(k))->GetDetector();
2156 Int_t rowMax = fGeo->GetRowMax(GetLayer(detector),GetStack(detector),GetSector(detector));
2157 Int_t colMax = fGeo->GetColMax(GetLayer(detector));
2158 Float_t min = 100.0;
2159 for (Int_t row = 0; row < rowMax; row++) {
2160 for (Int_t col = 0; col < colMax; col++) {
2161 value = ((AliTRDFitInfo *) fVectorFit2.At(k))->GetCoef()[(Int_t)(col*rowMax+row)];
2162 mean += -TMath::Abs(value);
2166 if(count > 0) mean = mean/count;
2168 value = ((AliTRDFitInfo *) vectorFit->At(k))->GetCoef()[0];
2169 object->SetValue(detector,-TMath::Abs(value));
2175 //_____________________________________________________________________________
2176 TObject *AliTRDCalibraFit::CreatePadObjectGain(const TObjArray *vectorFit, Double_t scaleFitFactor, const AliTRDCalDet *detobject)
2179 // It Creates the AliTRDCalPad object from AliTRDFitInfo
2180 // You need first to create the object for the detectors,
2181 // where the mean value is put.
2182 // This object has to be written in the database
2185 // Create the DetObject
2186 AliTRDCalPad *object = new AliTRDCalPad("GainFactor","GainFactor (local variations)");
2189 for(Int_t k = 0; k < 540; k++){
2190 AliTRDCalROC *calROC = object->GetCalROC(k);
2191 Int_t nchannels = calROC->GetNchannels();
2192 for(Int_t ch = 0; ch < nchannels; ch++){
2193 calROC->SetValue(ch,1.0);
2199 Int_t loop = (Int_t) vectorFit->GetEntriesFast();
2200 if(loop != 540) AliInfo("The Vector Fit is not complete!");
2201 Int_t detector = -1;
2202 Float_t value = 0.0;
2204 for (Int_t k = 0; k < loop; k++) {
2205 detector = ((AliTRDFitInfo *) vectorFit->At(k))->GetDetector();
2206 AliTRDCalROC *calROC = object->GetCalROC(detector);
2207 Float_t mean = detobject->GetValue(detector);
2208 if(TMath::Abs(mean) <= 0.0000000001) continue;
2209 Int_t rowMax = calROC->GetNrows();
2210 Int_t colMax = calROC->GetNcols();
2211 for (Int_t row = 0; row < rowMax; row++) {
2212 for (Int_t col = 0; col < colMax; col++) {
2213 value = ((AliTRDFitInfo *) vectorFit->At(k))->GetCoef()[(Int_t)(col*rowMax+row)];
2214 if(value > 0) value = value*scaleFitFactor;
2215 calROC->SetValue(col,row,TMath::Abs(value)/mean);
2223 //_____________________________________________________________________________
2224 TObject *AliTRDCalibraFit::CreatePadObjectVdrift(const TObjArray *vectorFit, const AliTRDCalDet *detobject)
2227 // It Creates the AliTRDCalPad object from AliTRDFitInfo
2228 // You need first to create the object for the detectors,
2229 // where the mean value is put.
2230 // This object has to be written in the database
2233 // Create the DetObject
2234 AliTRDCalPad *object = new AliTRDCalPad("LocalVdrift","TRD drift velocities (local variations)");
2237 for(Int_t k = 0; k < 540; k++){
2238 AliTRDCalROC *calROC = object->GetCalROC(k);
2239 Int_t nchannels = calROC->GetNchannels();
2240 for(Int_t ch = 0; ch < nchannels; ch++){
2241 calROC->SetValue(ch,1.0);
2247 Int_t loop = (Int_t) vectorFit->GetEntriesFast();
2248 if(loop != 540) AliInfo("The Vector Fit is not complete!");
2249 Int_t detector = -1;
2250 Float_t value = 0.0;
2252 for (Int_t k = 0; k < loop; k++) {
2253 detector = ((AliTRDFitInfo *) vectorFit->At(k))->GetDetector();
2254 AliTRDCalROC *calROC = object->GetCalROC(detector);
2255 Float_t mean = detobject->GetValue(detector);
2256 if(mean == 0) continue;
2257 Int_t rowMax = calROC->GetNrows();
2258 Int_t colMax = calROC->GetNcols();
2259 for (Int_t row = 0; row < rowMax; row++) {
2260 for (Int_t col = 0; col < colMax; col++) {
2261 value = ((AliTRDFitInfo *) vectorFit->At(k))->GetCoef()[(Int_t)(col*rowMax+row)];
2262 calROC->SetValue(col,row,TMath::Abs(value)/mean);
2270 //_____________________________________________________________________________
2271 TObject *AliTRDCalibraFit::CreatePadObjectT0(const TObjArray *vectorFit, const AliTRDCalDet *detobject)
2274 // It Creates the AliTRDCalPad object from AliTRDFitInfo2
2275 // You need first to create the object for the detectors,
2276 // where the mean value is put.
2277 // This object has to be written in the database
2280 // Create the DetObject
2281 AliTRDCalPad *object = new AliTRDCalPad("LocalT0","T0 (local variations)");
2284 for(Int_t k = 0; k < 540; k++){
2285 AliTRDCalROC *calROC = object->GetCalROC(k);
2286 Int_t nchannels = calROC->GetNchannels();
2287 for(Int_t ch = 0; ch < nchannels; ch++){
2288 calROC->SetValue(ch,0.0);
2294 Int_t loop = (Int_t) vectorFit->GetEntriesFast();
2295 if(loop != 540) AliInfo("The Vector Fit is not complete!");
2296 Int_t detector = -1;
2297 Float_t value = 0.0;
2299 for (Int_t k = 0; k < loop; k++) {
2300 detector = ((AliTRDFitInfo *) vectorFit->At(k))->GetDetector();
2301 AliTRDCalROC *calROC = object->GetCalROC(detector);
2302 Float_t min = detobject->GetValue(detector);
2303 Int_t rowMax = calROC->GetNrows();
2304 Int_t colMax = calROC->GetNcols();
2305 for (Int_t row = 0; row < rowMax; row++) {
2306 for (Int_t col = 0; col < colMax; col++) {
2307 value = ((AliTRDFitInfo *) vectorFit->At(k))->GetCoef()[(Int_t)(col*rowMax+row)];
2309 if(value > 70.0) value = value - 100.0;
2311 calROC->SetValue(col,row,value-min);
2319 //_____________________________________________________________________________
2320 TObject *AliTRDCalibraFit::CreatePadObjectPRF(const TObjArray *vectorFit)
2323 // It Creates the AliTRDCalPad object from AliTRDFitInfo
2324 // This object has to be written in the database
2327 // Create the DetObject
2328 AliTRDCalPad *object = new AliTRDCalPad("PRFWidth","PRFWidth");
2330 Int_t loop = (Int_t) vectorFit->GetEntriesFast();
2331 if(loop != 540) AliInfo("The Vector Fit is not complete!");
2332 Int_t detector = -1;
2333 Float_t value = 0.0;
2335 for (Int_t k = 0; k < loop; k++) {
2336 detector = ((AliTRDFitInfo *) vectorFit->At(k))->GetDetector();
2337 AliTRDCalROC *calROC = object->GetCalROC(detector);
2338 Int_t rowMax = calROC->GetNrows();
2339 Int_t colMax = calROC->GetNcols();
2340 for (Int_t row = 0; row < rowMax; row++) {
2341 for (Int_t col = 0; col < colMax; col++) {
2342 value = ((AliTRDFitInfo *) vectorFit->At(k))->GetCoef()[(Int_t)(col*rowMax+row)];
2343 calROC->SetValue(col,row,TMath::Abs(value));
2351 //_____________________________________________________________________________
2352 AliTRDCalDet *AliTRDCalibraFit::MakeOutliersStatDet(const TObjArray *vectorFit, const char *name, Double_t &mean)
2355 // It Creates the AliTRDCalDet object from AliTRDFitInfo
2356 // 0 successful fit 1 not successful fit
2357 // mean is the mean value over the successful fit
2358 // do not use it for t0: no meaning
2361 // Create the CalObject
2362 AliTRDCalDet *object = new AliTRDCalDet(name,name);
2366 Int_t loop = (Int_t) vectorFit->GetEntriesFast();
2368 AliInfo("The Vector Fit is not complete! We initialise all outliers");
2369 for(Int_t k = 0; k < 540; k++){
2370 object->SetValue(k,1.0);
2373 Int_t detector = -1;
2374 Float_t value = 0.0;
2376 for (Int_t k = 0; k < loop; k++) {
2377 detector = ((AliTRDFitInfo *) vectorFit->At(k))->GetDetector();
2378 value = ((AliTRDFitInfo *) vectorFit->At(k))->GetCoef()[0];
2379 if(value <= 0) object->SetValue(detector,1.0);
2381 object->SetValue(detector,0.0);
2386 if(count > 0) mean /= count;
2389 //_____________________________________________________________________________
2390 TObject *AliTRDCalibraFit::MakeOutliersStatPad(const TObjArray *vectorFit, const char *name, Double_t &mean)
2393 // It Creates the AliTRDCalPad object from AliTRDFitInfo
2394 // 0 not successful fit 1 successful fit
2395 // mean mean value over the successful fit
2398 // Create the CalObject
2399 AliTRDCalPad *object = new AliTRDCalPad(name,name);
2403 Int_t loop = (Int_t) vectorFit->GetEntriesFast();
2405 AliInfo("The Vector Fit is not complete! We initialise all outliers");
2406 for(Int_t k = 0; k < 540; k++){
2407 AliTRDCalROC *calROC = object->GetCalROC(k);
2408 Int_t nchannels = calROC->GetNchannels();
2409 for(Int_t ch = 0; ch < nchannels; ch++){
2410 calROC->SetValue(ch,1.0);
2414 Int_t detector = -1;
2415 Float_t value = 0.0;
2417 for (Int_t k = 0; k < loop; k++) {
2418 detector = ((AliTRDFitInfo *) vectorFit->At(k))->GetDetector();
2419 AliTRDCalROC *calROC = object->GetCalROC(detector);
2420 Int_t nchannels = calROC->GetNchannels();
2421 for (Int_t ch = 0; ch < nchannels; ch++) {
2422 value = ((AliTRDFitInfo *) vectorFit->At(k))->GetCoef()[ch];
2423 if(value <= 0) calROC->SetValue(ch,1.0);
2425 calROC->SetValue(ch,0.0);
2431 if(count > 0) mean /= count;
2434 //_____________________________________________________________________________
2435 void AliTRDCalibraFit::SetPeriodeFitPH(Int_t periodeFitPH)
2438 // Set FitPH if 1 then each detector will be fitted
2441 if (periodeFitPH > 0) {
2442 fFitPHPeriode = periodeFitPH;
2445 AliInfo("periodeFitPH must be higher than 0!");
2449 //_____________________________________________________________________________
2450 void AliTRDCalibraFit::SetBeginFitCharge(Float_t beginFitCharge)
2453 // The fit of the deposited charge distribution begins at
2454 // histo->Mean()/beginFitCharge
2455 // You can here set beginFitCharge
2458 if (beginFitCharge > 0) {
2459 fBeginFitCharge = beginFitCharge;
2462 AliInfo("beginFitCharge must be strict positif!");
2467 //_____________________________________________________________________________
2468 void AliTRDCalibraFit::SetT0Shift0(Float_t t0Shift)
2471 // The t0 calculated with the maximum positif slope is shift from t0Shift0
2472 // You can here set t0Shift0
2476 fT0Shift0 = t0Shift;
2479 AliInfo("t0Shift0 must be strict positif!");
2484 //_____________________________________________________________________________
2485 void AliTRDCalibraFit::SetT0Shift1(Float_t t0Shift)
2488 // The t0 calculated with the maximum of the amplification region is shift from t0Shift1
2489 // You can here set t0Shift1
2493 fT0Shift1 = t0Shift;
2496 AliInfo("t0Shift must be strict positif!");
2501 //_____________________________________________________________________________
2502 void AliTRDCalibraFit::SetRangeFitPRF(Float_t rangeFitPRF)
2505 // The fit of the PRF is from -rangeFitPRF to rangeFitPRF
2506 // You can here set rangeFitPRF
2509 if ((rangeFitPRF > 0) &&
2510 (rangeFitPRF <= 1.5)) {
2511 fRangeFitPRF = rangeFitPRF;
2514 AliInfo("rangeFitPRF must be between 0 and 1.0");
2519 //_____________________________________________________________________________
2520 void AliTRDCalibraFit::SetMinEntries(Int_t minEntries)
2523 // Minimum entries for fitting
2526 if (minEntries > 0) {
2527 fMinEntries = minEntries;
2530 AliInfo("fMinEntries must be >= 0.");
2535 //_____________________________________________________________________________
2536 void AliTRDCalibraFit::SetRebin(Short_t rebin)
2539 // Rebin with rebin time less bins the Ch histo
2540 // You can set here rebin that should divide the number of bins of CH histo
2545 AliInfo("You have to be sure that fRebin divides fNumberBinCharge used!");
2548 AliInfo("You have to choose a positiv value!");
2552 //_____________________________________________________________________________
2553 Bool_t AliTRDCalibraFit::FillVectorFit()
2556 // For the Fit functions fill the vector Fit
2559 AliTRDFitInfo *fitInfo = new AliTRDFitInfo();
2562 if (GetStack(fCountDet) == 2) {
2569 //printf("For the detector %d , ntotal %d and fCoefCH[0] %f\n",countdet,ntotal,fCoefCH[0]);
2570 Float_t *coef = new Float_t[ntotal];
2571 for (Int_t i = 0; i < ntotal; i++) {
2572 coef[i] = fCurrentCoefDetector[i];
2575 Int_t detector = fCountDet;
2577 fitInfo->SetCoef(coef);
2578 fitInfo->SetDetector(detector);
2579 fVectorFit.Add((TObject *) fitInfo);
2584 //_____________________________________________________________________________
2585 Bool_t AliTRDCalibraFit::FillVectorFit2()
2588 // For the Fit functions fill the vector Fit
2591 AliTRDFitInfo *fitInfo = new AliTRDFitInfo();
2594 if (GetStack(fCountDet) == 2) {
2601 //printf("For the detector %d , ntotal %d and fCoefCH[0] %f\n",countdet,ntotal,fCoefCH[0]);
2602 Float_t *coef = new Float_t[ntotal];
2603 for (Int_t i = 0; i < ntotal; i++) {
2604 coef[i] = fCurrentCoefDetector2[i];
2607 Int_t detector = fCountDet;
2609 fitInfo->SetCoef(coef);
2610 fitInfo->SetDetector(detector);
2611 fVectorFit2.Add((TObject *) fitInfo);
2616 //____________Functions for initialising the AliTRDCalibraFit in the code_________
2617 Bool_t AliTRDCalibraFit::InitFit(Int_t nbins, Int_t i)
2620 // Init the number of expected bins and fDect1[i] fDect2[i]
2623 gStyle->SetPalette(1);
2624 gStyle->SetOptStat(1111);
2625 gStyle->SetPadBorderMode(0);
2626 gStyle->SetCanvasColor(10);
2627 gStyle->SetPadLeftMargin(0.13);
2628 gStyle->SetPadRightMargin(0.01);
2630 // Mode groups of pads: the total number of bins!
2631 CalculNumberOfBinsExpected(i);
2633 // Quick verification that we have the good pad calibration mode!
2634 if (fNumberOfBinsExpected != nbins) {
2635 AliInfo(Form("It doesn't correspond to the mode of pad group calibration: expected %d and seen %d!",fNumberOfBinsExpected,nbins));
2639 // Security for fDebug 3 and 4
2640 if ((fDebugLevel >= 3) &&
2644 AliInfo("This detector doesn't exit!");
2648 // Determine fDet1 and fDet2 and set the fNfragZ and fNfragRphi for debug 3 and 4
2649 CalculDect1Dect2(i);
2654 //____________Functions for initialising the AliTRDCalibraFit in the code_________
2655 Bool_t AliTRDCalibraFit::InitFitCH()
2658 // Init the fVectorFitCH for normalisation
2659 // Init the histo for debugging
2664 fScaleFitFactor = 0.0;
2665 fCurrentCoefDetector = new Float_t[2304];
2666 for (Int_t k = 0; k < 2304; k++) {
2667 fCurrentCoefDetector[k] = 0.0;
2669 fVectorFit.SetName("gainfactorscoefficients");
2671 // fDebug == 0 nothing
2672 // fDebug == 1 and fFitVoir no histo
2673 if (fDebugLevel == 1) {
2674 if(!CheckFitVoir()) return kFALSE;
2676 //Get the CalDet object
2678 AliTRDcalibDB *cal = AliTRDcalibDB::Instance();
2680 AliInfo("Could not get calibDB");
2683 if(fCalDet) delete fCalDet;
2684 fCalDet = new AliTRDCalDet(*(cal->GetGainFactorDet()));
2687 Float_t devalue = 1.0;
2688 if(fCalDet) delete fCalDet;
2689 fCalDet = new AliTRDCalDet("ChamberGainFactor","GainFactor (detector value)");
2690 for(Int_t k = 0; k < 540; k++){
2691 fCalDet->SetValue(k,devalue);
2697 //____________Functions for initialising the AliTRDCalibraFit in the code_________
2698 Bool_t AliTRDCalibraFit::InitFitPH()
2701 // Init the arrays of results
2702 // Init the histos for debugging
2706 fVectorFit.SetName("driftvelocitycoefficients");
2707 fVectorFit2.SetName("t0coefficients");
2709 fCurrentCoefDetector = new Float_t[2304];
2710 for (Int_t k = 0; k < 2304; k++) {
2711 fCurrentCoefDetector[k] = 0.0;
2714 fCurrentCoefDetector2 = new Float_t[2304];
2715 for (Int_t k = 0; k < 2304; k++) {
2716 fCurrentCoefDetector2[k] = 0.0;
2719 //fDebug == 0 nothing
2720 // fDebug == 1 and fFitVoir no histo
2721 if (fDebugLevel == 1) {
2722 if(!CheckFitVoir()) return kFALSE;
2724 //Get the CalDet object
2726 AliTRDcalibDB *cal = AliTRDcalibDB::Instance();
2728 AliInfo("Could not get calibDB");
2731 if(fCalDet) delete fCalDet;
2732 if(fCalDet2) delete fCalDet2;
2733 fCalDet = new AliTRDCalDet(*(cal->GetVdriftDet()));
2734 fCalDet2 = new AliTRDCalDet(*(cal->GetT0Det()));
2737 Float_t devalue = 1.5;
2738 Float_t devalue2 = 0.0;
2739 if(fCalDet) delete fCalDet;
2740 if(fCalDet2) delete fCalDet2;
2741 fCalDet = new AliTRDCalDet("ChamberVdrift","TRD drift velocities (detector value)");
2742 fCalDet2 = new AliTRDCalDet("ChamberT0","T0 (detector value)");
2743 for(Int_t k = 0; k < 540; k++){
2744 fCalDet->SetValue(k,devalue);
2745 fCalDet2->SetValue(k,devalue2);
2750 //____________Functions for initialising the AliTRDCalibraFit in the code_________
2751 Bool_t AliTRDCalibraFit::InitFitPRF()
2754 // Init the calibration mode (Nz, Nrphi), the histograms for
2755 // debugging the fit methods if fDebug > 0,
2759 fVectorFit.SetName("prfwidthcoefficients");
2761 fCurrentCoefDetector = new Float_t[2304];
2762 for (Int_t k = 0; k < 2304; k++) {
2763 fCurrentCoefDetector[k] = 0.0;
2766 // fDebug == 0 nothing
2767 // fDebug == 1 and fFitVoir no histo
2768 if (fDebugLevel == 1) {
2769 if(!CheckFitVoir()) return kFALSE;
2773 //____________Functions for initialising the AliTRDCalibraFit in the code_________
2774 Bool_t AliTRDCalibraFit::InitFitLinearFitter()
2777 // Init the fCalDet, fVectorFit fCurrentCoefDetector
2782 fCurrentCoefDetector = new Float_t[2304];
2783 fCurrentCoefDetector2 = new Float_t[2304];
2784 for (Int_t k = 0; k < 2304; k++) {
2785 fCurrentCoefDetector[k] = 0.0;
2786 fCurrentCoefDetector2[k] = 0.0;
2789 //printf("test0\n");
2791 AliTRDcalibDB *cal = AliTRDcalibDB::Instance();
2793 AliInfo("Could not get calibDB");
2797 //Get the CalDet object
2799 if(fCalDet) delete fCalDet;
2800 if(fCalDet2) delete fCalDet2;
2801 fCalDet = new AliTRDCalDet(*(cal->GetVdriftDet()));
2802 //printf("test1\n");
2803 fCalDet2 = new AliTRDCalDet("lorentz angle tan","lorentz angle tan (detector value)");
2804 //printf("test2\n");
2805 for(Int_t k = 0; k < 540; k++){
2806 fCalDet2->SetValue(k,AliTRDCommonParam::Instance()->GetOmegaTau(fCalDet->GetValue(k)));
2808 //printf("test3\n");
2811 Float_t devalue = 1.5;
2812 Float_t devalue2 = AliTRDCommonParam::Instance()->GetOmegaTau(1.5);
2813 if(fCalDet) delete fCalDet;
2814 if(fCalDet2) delete fCalDet2;
2815 //printf("test1\n");
2816 fCalDet = new AliTRDCalDet("ChamberVdrift","TRD drift velocities (detector value)");
2817 fCalDet2 = new AliTRDCalDet("lorentz angle tan","lorentz angle tan (detector value)");
2818 //printf("test2\n");
2819 for(Int_t k = 0; k < 540; k++){
2820 fCalDet->SetValue(k,devalue);
2821 fCalDet2->SetValue(k,devalue2);
2823 //printf("test3\n");
2828 //____________Functions for initialising the AliTRDCalibraFit in the code_________
2829 void AliTRDCalibraFit::InitfCountDetAndfCount(Int_t i)
2832 // Init the current detector where we are fCountDet and the
2833 // next fCount for the functions Fit...
2836 // Loop on the Xbins of ch!!
2837 fCountDet = -1; // Current detector
2838 fCount = 0; // To find the next detector
2841 if (fDebugLevel >= 3) {
2842 // Set countdet to the detector
2843 fCountDet = AliTRDgeometry::GetDetector(fDet[0],fDet[1],fDet[2]);
2844 // Set counter to write at the end of the detector
2846 // Get the right calib objects
2849 if(fDebugLevel == 1) {
2851 fCalibraMode->CalculXBins(fCountDet,i);
2852 if((fCalibraMode->GetNz(i)!=100) && (fCalibraMode->GetNrphi(i)!=100)){
2853 while(fCalibraMode->GetXbins(i) <=fFitVoir){
2855 fCalibraMode->CalculXBins(fCountDet,i);
2856 //printf("GetXBins %d\n",fCalibraMode->GetXbins(i));
2862 fCount = fCalibraMode->GetXbins(i);
2864 // Determination of fNnZ, fNnRphi, fNfragZ and fNfragRphi
2865 fCalibraMode->ModePadCalibration((Int_t) GetStack(fCountDet),i);
2866 fCalibraMode->ModePadFragmentation((Int_t) GetLayer(fCountDet)
2867 ,(Int_t) GetStack(fCountDet)
2868 ,(Int_t) GetSector(fCountDet),i);
2871 //_______________________________________________________________________________
2872 void AliTRDCalibraFit::CalculNumberOfBinsExpected(Int_t i)
2875 // Calculate the number of bins expected (calibration groups)
2878 fNumberOfBinsExpected = 0;
2880 if((fCalibraMode->GetNz(i) == 100) && (fCalibraMode->GetNrphi(i) == 100)){
2881 fNumberOfBinsExpected = 1;
2885 if((fCalibraMode->GetNz(i) == 10) && (fCalibraMode->GetNrphi(i) == 10)){
2886 fNumberOfBinsExpected = 18;
2890 fCalibraMode->ModePadCalibration(2,i);
2891 fCalibraMode->ModePadFragmentation(0,2,0,i);
2892 fCalibraMode->SetDetChamb2(i);
2893 if (fDebugLevel > 1) {
2894 AliInfo(Form("For the chamber 2: %d",fCalibraMode->GetDetChamb2(i)));
2896 fNumberOfBinsExpected += 6 * 18 * fCalibraMode->GetDetChamb2(i);
2897 fCalibraMode->ModePadCalibration(0,i);
2898 fCalibraMode->ModePadFragmentation(0,0,0,i);
2899 fCalibraMode->SetDetChamb0(i);
2900 if (fDebugLevel > 1) {
2901 AliInfo(Form("For the other chamber 0: %d",fCalibraMode->GetDetChamb0(i)));
2903 fNumberOfBinsExpected += 6 * 4 * 18 * fCalibraMode->GetDetChamb0(i);
2906 //_______________________________________________________________________________
2907 void AliTRDCalibraFit::CalculDect1Dect2(Int_t i)
2910 // Calculate the range of fits
2915 if (fDebugLevel == 1) {
2919 if ((fDebugLevel == 2) || (fDebugLevel == 0)) {
2921 fDect2 = fNumberOfBinsExpected;
2923 if (fDebugLevel >= 3) {
2924 fCountDet = AliTRDgeometry::GetDetector(fDet[0],fDet[1],fDet[2]);
2925 fCalibraMode->CalculXBins(fCountDet,i);
2926 fDect1 = fCalibraMode->GetXbins(i);
2927 // Determination of fNnZ, fNnRphi, fNfragZ and fNfragRphi
2928 fCalibraMode->ModePadCalibration((Int_t) GetStack(fCountDet),i);
2929 fCalibraMode->ModePadFragmentation((Int_t) GetLayer(fCountDet)
2930 ,(Int_t) GetStack(fCountDet)
2931 ,(Int_t) GetSector(fCountDet),i);
2932 // Set for the next detector
2933 fDect2 = fDect1 + fCalibraMode->GetNfragZ(i)*fCalibraMode->GetNfragRphi(i);
2936 //_______________________________________________________________________________
2937 Bool_t AliTRDCalibraFit::CheckFitVoir()
2940 // Check if fFitVoir is in the range
2943 if (fFitVoir < fNumberOfBinsExpected) {
2944 AliInfo(Form("We will see the fit of the object %d",fFitVoir));
2947 AliInfo("fFitVoir is out of range of the histo!");
2952 //____________Functions for initialising the AliTRDCalibraFit in the code_________
2953 void AliTRDCalibraFit::UpdatefCountDetAndfCount(Int_t idect, Int_t i)
2956 // See if we are in a new detector and update the
2957 // variables fNfragZ and fNfragRphi if yes
2958 // Will never happen for only one detector (3 and 4)
2959 // Doesn't matter for 2
2961 if (fCount == idect) {
2962 // On en est au detector (or first detector in the group)
2964 AliDebug(2,Form("We are at the detector %d\n",fCountDet));
2965 // Determination of fNnZ, fNnRphi, fNfragZ and fNfragRphi
2966 fCalibraMode->ModePadCalibration((Int_t) GetStack(fCountDet),i);
2967 fCalibraMode->ModePadFragmentation((Int_t) GetLayer(fCountDet)
2968 ,(Int_t) GetStack(fCountDet)
2969 ,(Int_t) GetSector(fCountDet),i);
2970 // Set for the next detector
2971 fCount += fCalibraMode->GetNfragZ(i)*fCalibraMode->GetNfragRphi(i);
2976 //____________Functions for initialising the AliTRDCalibraFit in the code_________
2977 void AliTRDCalibraFit::ReconstructFitRowMinRowMax(Int_t idect, Int_t i)
2980 // Reconstruct the min pad row, max pad row, min pad col and
2981 // max pad col of the calibration group for the Fit functions
2982 // idect is the calibration group inside the detector
2984 if (fDebugLevel != 1) {
2985 fCalibraMode->ReconstructionRowPadGroup((Int_t) (idect-(fCount-(fCalibraMode->GetNfragZ(i)*fCalibraMode->GetNfragRphi(i)))),i);
2987 AliDebug(2,Form("AliTRDCalibraFit::ReconstructFitRowMinRowMax: the local calibration group is %d",idect-(fCount-(fCalibraMode->GetNfragZ(i)*fCalibraMode->GetNfragRphi(i)))));
2988 AliDebug(2,Form("AliTRDCalibraFit::ReconstructFitRowMinRowMax: the number of group per detector is %d",fCalibraMode->GetNfragZ(i)*fCalibraMode->GetNfragRphi(i)));
2990 //____________Functions for initialising the AliTRDCalibraFit in the code_________
2991 Bool_t AliTRDCalibraFit::NotEnoughStatisticCH(Int_t idect)
2994 // For the case where there are not enough entries in the histograms
2995 // of the calibration group, the value present in the choosen database
2996 // will be put. A negativ sign enables to know that a fit was not possible.
2999 if (fDebugLevel == 1) {
3000 AliInfo("The element has not enough statistic to be fitted");
3002 else if (fNbDet > 0){
3003 Int_t firstdetector = fCountDet;
3004 Int_t lastdetector = fCountDet+fNbDet;
3005 AliInfo(Form("The element %d containing the detectors %d to %d has not enough statistic to be fitted"
3006 ,idect,firstdetector,lastdetector));
3007 // loop over detectors
3008 for(Int_t det = firstdetector; det < lastdetector; det++){
3010 //Set the calibration object again
3014 // Determination of fNnZ, fNnRphi, fNfragZ and fNfragRphi
3016 fCalibraMode->ModePadCalibration((Int_t) GetStack(fCountDet),0);
3017 fCalibraMode->ModePadFragmentation((Int_t) GetLayer(fCountDet)
3018 ,(Int_t) GetStack(fCountDet)
3019 ,(Int_t) GetSector(fCountDet),0);
3020 // Reconstruct row min row max
3021 ReconstructFitRowMinRowMax(idect,0);
3023 // Calcul the coef from the database choosen for the detector
3024 CalculChargeCoefMean(kFALSE);
3026 //stack 2, not stack 2
3028 if(GetStack(fCountDet) == 2) factor = 12;
3031 // Fill the fCurrentCoefDetector with negative value to say: not fitted
3032 for (Int_t k = fCalibraMode->GetRowMin(0); k < fCalibraMode->GetRowMax(0); k++) {
3033 for (Int_t j = fCalibraMode->GetColMin(0); j < fCalibraMode->GetColMax(0); j++) {
3034 fCurrentCoefDetector[(Int_t)(j*factor+k)] = -TMath::Abs(fCurrentCoef[1]);
3038 //Put default value negative
3039 fCurrentCoef[0] = -TMath::Abs(fCurrentCoef[1]);
3040 fCurrentCoefE = 0.0;
3045 if(fDebugLevel > 1){
3047 if ( !fDebugStreamer ) {
3049 TDirectory *backup = gDirectory;
3050 fDebugStreamer = new TTreeSRedirector("TRDDebugFitCH.root");
3051 if ( backup ) backup->cd(); //we don't want to be cd'd to the debug streamer
3054 Int_t detector = fCountDet;
3055 Int_t caligroup = idect;
3056 Short_t rowmin = fCalibraMode->GetRowMin(0);
3057 Short_t rowmax = fCalibraMode->GetRowMax(0);
3058 Short_t colmin = fCalibraMode->GetColMin(0);
3059 Short_t colmax = fCalibraMode->GetColMax(0);
3060 Float_t gf = fCurrentCoef[0];
3061 Float_t gfs = fCurrentCoef[1];
3062 Float_t gfE = fCurrentCoefE;
3064 (*fDebugStreamer) << "FillFillCH" <<
3065 "detector=" << detector <<
3066 "caligroup=" << caligroup <<
3067 "rowmin=" << rowmin <<
3068 "rowmax=" << rowmax <<
3069 "colmin=" << colmin <<
3070 "colmax=" << colmax <<
3078 for (Int_t k = 0; k < 2304; k++) {
3079 fCurrentCoefDetector[k] = 0.0;
3083 AliDebug(2,Form("Check the count now: fCountDet %d",fCountDet));
3087 AliInfo(Form("The element %d in this detector %d has not enough statistic to be fitted"
3088 ,idect-(fCount-(fCalibraMode->GetNfragZ(0)*fCalibraMode->GetNfragRphi(0))),fCountDet));
3090 // Calcul the coef from the database choosen
3091 CalculChargeCoefMean(kFALSE);
3093 //stack 2, not stack 2
3095 if(GetStack(fCountDet) == 2) factor = 12;
3098 // Fill the fCurrentCoefDetector with negative value to say: not fitted
3099 for (Int_t k = fCalibraMode->GetRowMin(0); k < fCalibraMode->GetRowMax(0); k++) {
3100 for (Int_t j = fCalibraMode->GetColMin(0); j < fCalibraMode->GetColMax(0); j++) {
3101 fCurrentCoefDetector[(Int_t)(j*factor+k)] = -TMath::Abs(fCurrentCoef[1]);
3105 //Put default value negative
3106 fCurrentCoef[0] = -TMath::Abs(fCurrentCoef[1]);
3107 fCurrentCoefE = 0.0;
3116 //____________Functions for initialising the AliTRDCalibraFit in the code_________
3117 Bool_t AliTRDCalibraFit::NotEnoughStatisticPH(Int_t idect,Double_t nentries)
3120 // For the case where there are not enough entries in the histograms
3121 // of the calibration group, the value present in the choosen database
3122 // will be put. A negativ sign enables to know that a fit was not possible.
3124 if (fDebugLevel == 1) {
3125 AliInfo("The element has not enough statistic to be fitted");
3127 else if (fNbDet > 0) {
3129 Int_t firstdetector = fCountDet;
3130 Int_t lastdetector = fCountDet+fNbDet;
3131 AliInfo(Form("The element %d containing the detectors %d to %d has not enough statistic to be fitted"
3132 ,idect,firstdetector,lastdetector));
3133 // loop over detectors
3134 for(Int_t det = firstdetector; det < lastdetector; det++){
3136 //Set the calibration object again
3140 // Determination of fNnZ, fNnRphi, fNfragZ and fNfragRphi
3142 fCalibraMode->ModePadCalibration((Int_t) GetStack(fCountDet),1);
3143 fCalibraMode->ModePadFragmentation((Int_t) GetLayer(fCountDet)
3144 ,(Int_t) GetStack(fCountDet)
3145 ,(Int_t) GetSector(fCountDet),1);
3146 // Reconstruct row min row max
3147 ReconstructFitRowMinRowMax(idect,1);
3149 // Calcul the coef from the database choosen for the detector
3150 CalculVdriftCoefMean();
3153 //stack 2, not stack 2
3155 if(GetStack(fCountDet) == 2) factor = 12;
3158 // Fill the fCurrentCoefDetector with negative value to say: not fitted
3159 for (Int_t k = fCalibraMode->GetRowMin(1); k < fCalibraMode->GetRowMax(1); k++) {
3160 for (Int_t j = fCalibraMode->GetColMin(1); j < fCalibraMode->GetColMax(1); j++) {
3161 fCurrentCoefDetector[(Int_t)(j*factor+k)] = -TMath::Abs(fCurrentCoef[1]);
3162 fCurrentCoefDetector2[(Int_t)(j*factor+k)] = fCurrentCoef2[1] + 100.0;
3166 //Put default value negative
3167 fCurrentCoef[0] = -TMath::Abs(fCurrentCoef[1]);
3168 fCurrentCoefE = 0.0;
3169 fCurrentCoef2[0] = fCurrentCoef2[1] + 100.0;
3170 fCurrentCoefE2 = 0.0;
3176 if(fDebugLevel > 1){
3178 if ( !fDebugStreamer ) {
3180 TDirectory *backup = gDirectory;
3181 fDebugStreamer = new TTreeSRedirector("TRDDebugFitPH.root");
3182 if ( backup ) backup->cd(); //we don't want to be cd'd to the debug streamer
3186 Int_t detector = fCountDet;
3187 Int_t caligroup = idect;
3188 Short_t rowmin = fCalibraMode->GetRowMin(1);
3189 Short_t rowmax = fCalibraMode->GetRowMax(1);
3190 Short_t colmin = fCalibraMode->GetColMin(1);
3191 Short_t colmax = fCalibraMode->GetColMax(1);
3192 Float_t vf = fCurrentCoef[0];
3193 Float_t vs = fCurrentCoef[1];
3194 Float_t vfE = fCurrentCoefE;
3195 Float_t t0f = fCurrentCoef2[0];
3196 Float_t t0s = fCurrentCoef2[1];
3197 Float_t t0E = fCurrentCoefE2;
3201 (* fDebugStreamer) << "FillFillPH"<<
3202 "detector="<<detector<<
3203 "nentries="<<nentries<<
3204 "caligroup="<<caligroup<<
3218 for (Int_t k = 0; k < 2304; k++) {
3219 fCurrentCoefDetector[k] = 0.0;
3220 fCurrentCoefDetector2[k] = 0.0;
3224 AliDebug(2,Form("Check the count now: fCountDet %d",fCountDet));
3228 AliInfo(Form("The element %d in this detector %d has not enough statistic to be fitted"
3229 ,idect-(fCount-(fCalibraMode->GetNfragZ(1)*fCalibraMode->GetNfragRphi(1))),fCountDet));
3231 CalculVdriftCoefMean();
3234 //stack 2 and not stack 2
3236 if(GetStack(fCountDet) == 2) factor = 12;
3240 // Fill the fCurrentCoefDetector 2
3241 for (Int_t k = fCalibraMode->GetRowMin(1); k < fCalibraMode->GetRowMax(1); k++) {
3242 for (Int_t j = fCalibraMode->GetColMin(1); j < fCalibraMode->GetColMax(1); j++) {
3243 fCurrentCoefDetector[(Int_t)(j*factor+k)] = -TMath::Abs(fCurrentCoef[1]);
3244 fCurrentCoefDetector2[(Int_t)(j*factor+k)] = fCurrentCoef2[1] + 100.0;
3248 // Put the default value
3249 fCurrentCoef[0] = -TMath::Abs(fCurrentCoef[1]);
3250 fCurrentCoefE = 0.0;
3251 fCurrentCoef2[0] = fCurrentCoef2[1] + 100.0;
3252 fCurrentCoefE2 = 0.0;
3254 FillFillPH(idect,nentries);
3263 //____________Functions for initialising the AliTRDCalibraFit in the code_________
3264 Bool_t AliTRDCalibraFit::NotEnoughStatisticPRF(Int_t idect)
3267 // For the case where there are not enough entries in the histograms
3268 // of the calibration group, the value present in the choosen database
3269 // will be put. A negativ sign enables to know that a fit was not possible.
3272 if (fDebugLevel == 1) {
3273 AliInfo("The element has not enough statistic to be fitted");
3275 else if (fNbDet > 0){
3277 Int_t firstdetector = fCountDet;
3278 Int_t lastdetector = fCountDet+fNbDet;
3279 AliInfo(Form("The element %d containing the detectors %d to %d has not enough statistic to be fitted"
3280 ,idect,firstdetector,lastdetector));
3282 // loop over detectors
3283 for(Int_t det = firstdetector; det < lastdetector; det++){
3285 //Set the calibration object again
3289 // Determination of fNnZ, fNnRphi, fNfragZ and fNfragRphi
3291 fCalibraMode->ModePadCalibration((Int_t) GetStack(fCountDet),2);
3292 fCalibraMode->ModePadFragmentation((Int_t) GetLayer(fCountDet)
3293 ,(Int_t) GetStack(fCountDet)
3294 ,(Int_t) GetSector(fCountDet),2);
3295 // Reconstruct row min row max
3296 ReconstructFitRowMinRowMax(idect,2);
3298 // Calcul the coef from the database choosen for the detector
3299 CalculPRFCoefMean();
3301 //stack 2, not stack 2
3303 if(GetStack(fCountDet) == 2) factor = 12;
3306 // Fill the fCurrentCoefDetector with negative value to say: not fitted
3307 for (Int_t k = fCalibraMode->GetRowMin(2); k < fCalibraMode->GetRowMax(2); k++) {
3308 for (Int_t j = fCalibraMode->GetColMin(2); j < fCalibraMode->GetColMax(2); j++) {
3309 fCurrentCoefDetector[(Int_t)(j*factor+k)] = -TMath::Abs(fCurrentCoef[1]);
3313 //Put default value negative
3314 fCurrentCoef[0] = -TMath::Abs(fCurrentCoef[1]);
3315 fCurrentCoefE = 0.0;
3320 if(fDebugLevel > 1){
3322 if ( !fDebugStreamer ) {
3324 TDirectory *backup = gDirectory;
3325 fDebugStreamer = new TTreeSRedirector("TRDDebugFitPRF.root");
3326 if ( backup ) backup->cd(); //we don't want to be cd'd to the debug streamer
3329 Int_t detector = fCountDet;
3330 Int_t layer = GetLayer(fCountDet);
3331 Int_t caligroup = idect;
3332 Short_t rowmin = fCalibraMode->GetRowMin(2);
3333 Short_t rowmax = fCalibraMode->GetRowMax(2);
3334 Short_t colmin = fCalibraMode->GetColMin(2);
3335 Short_t colmax = fCalibraMode->GetColMax(2);
3336 Float_t widf = fCurrentCoef[0];
3337 Float_t wids = fCurrentCoef[1];
3338 Float_t widfE = fCurrentCoefE;
3340 (* fDebugStreamer) << "FillFillPRF"<<
3341 "detector="<<detector<<
3343 "caligroup="<<caligroup<<
3354 for (Int_t k = 0; k < 2304; k++) {
3355 fCurrentCoefDetector[k] = 0.0;
3359 AliDebug(2,Form("Check the count now: fCountDet %d",fCountDet));
3363 AliInfo(Form("The element %d in this detector %d has not enough statistic to be fitted"
3364 ,idect-(fCount-(fCalibraMode->GetNfragZ(2)*fCalibraMode->GetNfragRphi(2))),fCountDet));
3366 CalculPRFCoefMean();
3368 // stack 2 and not stack 2
3370 if(GetStack(fCountDet) == 2) factor = 12;
3374 // Fill the fCurrentCoefDetector
3375 for (Int_t k = fCalibraMode->GetRowMin(2); k < fCalibraMode->GetRowMax(2); k++) {
3376 for (Int_t j = fCalibraMode->GetColMin(2); j < fCalibraMode->GetColMax(2); j++) {
3377 fCurrentCoefDetector[(Int_t)(j*factor+k)] = -TMath::Abs(fCurrentCoef[1]);
3381 // Put the default value
3382 fCurrentCoef[0] = -TMath::Abs(fCurrentCoef[1]);
3383 fCurrentCoefE = 0.0;
3391 //____________Functions for initialising the AliTRDCalibraFit in the code_________
3392 Bool_t AliTRDCalibraFit::NotEnoughStatisticLinearFitter()
3395 // For the case where there are not enough entries in the histograms
3396 // of the calibration group, the value present in the choosen database
3397 // will be put. A negativ sign enables to know that a fit was not possible.
3400 // Calcul the coef from the database choosen
3401 CalculVdriftLorentzCoef();
3404 if(GetStack(fCountDet) == 2) factor = 1728;
3408 // Fill the fCurrentCoefDetector
3409 for (Int_t k = 0; k < factor; k++) {
3410 fCurrentCoefDetector[k] = -TMath::Abs(fCurrentCoef[1]);
3411 // should be negative
3412 fCurrentCoefDetector2[k] = +TMath::Abs(fCurrentCoef2[1]);
3416 //Put default opposite sign
3417 fCurrentCoef[0] = -TMath::Abs(fCurrentCoef[1]);
3418 fCurrentCoefE = 0.0;
3419 fCurrentCoef2[0] = +TMath::Abs(fCurrentCoef2[1]);
3420 fCurrentCoefE2 = 0.0;
3422 FillFillLinearFitter();
3427 //____________Functions for initialising the AliTRDCalibraFit in the code_________
3428 Bool_t AliTRDCalibraFit::FillInfosFitCH(Int_t idect)
3431 // Fill the coefficients found with the fits or other
3432 // methods from the Fit functions
3435 if (fDebugLevel != 1) {
3437 Int_t firstdetector = fCountDet;
3438 Int_t lastdetector = fCountDet+fNbDet;
3439 AliInfo(Form("The element %d containing the detectors %d to %d has been fitted"
3440 ,idect,firstdetector,lastdetector));
3441 // loop over detectors
3442 for(Int_t det = firstdetector; det < lastdetector; det++){
3444 //Set the calibration object again
3448 // Determination of fNnZ, fNnRphi, fNfragZ and fNfragRphi
3450 fCalibraMode->ModePadCalibration((Int_t) GetStack(fCountDet),0);
3451 fCalibraMode->ModePadFragmentation((Int_t) GetLayer(fCountDet)
3452 ,(Int_t) GetStack(fCountDet)
3453 ,(Int_t) GetSector(fCountDet),0);
3454 // Reconstruct row min row max
3455 ReconstructFitRowMinRowMax(idect,0);
3457 // Calcul the coef from the database choosen for the detector
3458 if(fCurrentCoef[0] < 0.0) CalculChargeCoefMean(kFALSE);
3459 else CalculChargeCoefMean(kTRUE);
3461 //stack 2, not stack 2
3463 if(GetStack(fCountDet) == 2) factor = 12;
3466 // Fill the fCurrentCoefDetector with negative value to say: not fitted
3467 Double_t coeftoput = 1.0;
3468 if(fCurrentCoef[0] < 0.0) coeftoput = - TMath::Abs(fCurrentCoef[1]);
3469 else coeftoput = fCurrentCoef[0];
3470 for (Int_t k = fCalibraMode->GetRowMin(0); k < fCalibraMode->GetRowMax(0); k++) {
3471 for (Int_t j = fCalibraMode->GetColMin(0); j < fCalibraMode->GetColMax(0); j++) {
3472 fCurrentCoefDetector[(Int_t)(j*factor+k)] = coeftoput;
3479 if(fDebugLevel > 1){
3481 if ( !fDebugStreamer ) {
3483 TDirectory *backup = gDirectory;
3484 fDebugStreamer = new TTreeSRedirector("TRDDebugFitCH.root");
3485 if ( backup ) backup->cd(); //we don't want to be cd'd to the debug streamer
3488 Int_t detector = fCountDet;
3489 Int_t caligroup = idect;
3490 Short_t rowmin = fCalibraMode->GetRowMin(0);
3491 Short_t rowmax = fCalibraMode->GetRowMax(0);
3492 Short_t colmin = fCalibraMode->GetColMin(0);
3493 Short_t colmax = fCalibraMode->GetColMax(0);
3494 Float_t gf = fCurrentCoef[0];
3495 Float_t gfs = fCurrentCoef[1];
3496 Float_t gfE = fCurrentCoefE;
3498 (*fDebugStreamer) << "FillFillCH" <<
3499 "detector=" << detector <<
3500 "caligroup=" << caligroup <<
3501 "rowmin=" << rowmin <<
3502 "rowmax=" << rowmax <<
3503 "colmin=" << colmin <<
3504 "colmax=" << colmax <<
3512 for (Int_t k = 0; k < 2304; k++) {
3513 fCurrentCoefDetector[k] = 0.0;
3517 //printf("Check the count now: fCountDet %d\n",fCountDet);
3522 if(GetStack(fCountDet) == 2) factor = 12;
3525 for (Int_t k = fCalibraMode->GetRowMin(0); k < fCalibraMode->GetRowMax(0); k++) {
3526 for (Int_t j = fCalibraMode->GetColMin(0); j < fCalibraMode->GetColMax(0); j++) {
3527 fCurrentCoefDetector[(Int_t)(j*factor+k)] = fCurrentCoef[0];
3538 //____________Functions for initialising the AliTRDCalibraFit in the code_________
3539 Bool_t AliTRDCalibraFit::FillInfosFitPH(Int_t idect,Double_t nentries)
3542 // Fill the coefficients found with the fits or other
3543 // methods from the Fit functions
3546 if (fDebugLevel != 1) {
3549 Int_t firstdetector = fCountDet;
3550 Int_t lastdetector = fCountDet+fNbDet;
3551 AliInfo(Form("The element %d containing the detectors %d to %d has been fitted"
3552 ,idect,firstdetector,lastdetector));
3554 // loop over detectors
3555 for(Int_t det = firstdetector; det < lastdetector; det++){
3557 //Set the calibration object again
3561 // Determination of fNnZ, fNnRphi, fNfragZ and fNfragRphi
3563 fCalibraMode->ModePadCalibration((Int_t) GetStack(fCountDet),1);
3564 fCalibraMode->ModePadFragmentation((Int_t) GetLayer(fCountDet)
3565 ,(Int_t) GetStack(fCountDet)
3566 ,(Int_t) GetSector(fCountDet),1);
3567 // Reconstruct row min row max
3568 ReconstructFitRowMinRowMax(idect,1);
3570 // Calcul the coef from the database choosen for the detector
3571 CalculVdriftCoefMean();
3574 //stack 2, not stack 2
3576 if(GetStack(fCountDet) == 2) factor = 12;
3579 // Fill the fCurrentCoefDetector with negative value to say: not fitted
3580 Double_t coeftoput = 1.5;
3581 Double_t coeftoput2 = 0.0;
3583 if(fCurrentCoef[0] < 0.0) coeftoput = - TMath::Abs(fCurrentCoef[1]);
3584 else coeftoput = fCurrentCoef[0];
3586 if(fCurrentCoef2[0] > 70.0) coeftoput2 = fCurrentCoef2[1] + 100.0;
3587 else coeftoput2 = fCurrentCoef2[0];
3589 for (Int_t k = fCalibraMode->GetRowMin(1); k < fCalibraMode->GetRowMax(1); k++) {
3590 for (Int_t j = fCalibraMode->GetColMin(1); j < fCalibraMode->GetColMax(1); j++) {
3591 fCurrentCoefDetector[(Int_t)(j*factor+k)] = coeftoput;
3592 fCurrentCoefDetector2[(Int_t)(j*factor+k)] = coeftoput2;
3600 if(fDebugLevel > 1){
3602 if ( !fDebugStreamer ) {
3604 TDirectory *backup = gDirectory;
3605 fDebugStreamer = new TTreeSRedirector("TRDDebugFitPH.root");
3606 if ( backup ) backup->cd(); //we don't want to be cd'd to the debug streamer
3610 Int_t detector = fCountDet;
3611 Int_t caligroup = idect;
3612 Short_t rowmin = fCalibraMode->GetRowMin(1);
3613 Short_t rowmax = fCalibraMode->GetRowMax(1);
3614 Short_t colmin = fCalibraMode->GetColMin(1);
3615 Short_t colmax = fCalibraMode->GetColMax(1);
3616 Float_t vf = fCurrentCoef[0];
3617 Float_t vs = fCurrentCoef[1];
3618 Float_t vfE = fCurrentCoefE;
3619 Float_t t0f = fCurrentCoef2[0];
3620 Float_t t0s = fCurrentCoef2[1];
3621 Float_t t0E = fCurrentCoefE2;
3625 (* fDebugStreamer) << "FillFillPH"<<
3626 "detector="<<detector<<
3627 "nentries="<<nentries<<
3628 "caligroup="<<caligroup<<
3642 for (Int_t k = 0; k < 2304; k++) {
3643 fCurrentCoefDetector[k] = 0.0;
3644 fCurrentCoefDetector2[k] = 0.0;
3648 //printf("Check the count now: fCountDet %d\n",fCountDet);
3653 if(GetStack(fCountDet) == 2) factor = 12;
3656 for (Int_t k = fCalibraMode->GetRowMin(1); k < fCalibraMode->GetRowMax(1); k++) {
3657 for (Int_t j = fCalibraMode->GetColMin(1); j < fCalibraMode->GetColMax(1); j++) {
3658 fCurrentCoefDetector[(Int_t)(j*factor+k)] = fCurrentCoef[0];
3659 fCurrentCoefDetector2[(Int_t)(j*factor+k)] = fCurrentCoef2[0];
3663 FillFillPH(idect,nentries);
3668 //____________Functions for initialising the AliTRDCalibraFit in the code_________
3669 Bool_t AliTRDCalibraFit::FillInfosFitPRF(Int_t idect)
3672 // Fill the coefficients found with the fits or other
3673 // methods from the Fit functions
3676 if (fDebugLevel != 1) {
3679 Int_t firstdetector = fCountDet;
3680 Int_t lastdetector = fCountDet+fNbDet;
3681 AliInfo(Form("The element %d containing the detectors %d to %d has been fitted"
3682 ,idect,firstdetector,lastdetector));
3684 // loop over detectors
3685 for(Int_t det = firstdetector; det < lastdetector; det++){
3687 //Set the calibration object again
3691 // Determination of fNnZ, fNnRphi, fNfragZ and fNfragRphi
3693 fCalibraMode->ModePadCalibration((Int_t) GetStack(fCountDet),2);
3694 fCalibraMode->ModePadFragmentation((Int_t) GetLayer(fCountDet)
3695 ,(Int_t) GetStack(fCountDet)
3696 ,(Int_t) GetSector(fCountDet),2);
3697 // Reconstruct row min row max
3698 ReconstructFitRowMinRowMax(idect,2);
3700 // Calcul the coef from the database choosen for the detector
3701 CalculPRFCoefMean();
3703 //stack 2, not stack 2
3705 if(GetStack(fCountDet) == 2) factor = 12;
3708 // Fill the fCurrentCoefDetector with negative value to say: not fitted
3709 Double_t coeftoput = 1.0;
3710 if(fCurrentCoef[0] < 0.0) coeftoput = - TMath::Abs(fCurrentCoef[1]);
3711 else coeftoput = fCurrentCoef[0];
3712 for (Int_t k = fCalibraMode->GetRowMin(2); k < fCalibraMode->GetRowMax(2); k++) {
3713 for (Int_t j = fCalibraMode->GetColMin(2); j < fCalibraMode->GetColMax(2); j++) {
3714 fCurrentCoefDetector[(Int_t)(j*factor+k)] = coeftoput;
3721 if(fDebugLevel > 1){
3723 if ( !fDebugStreamer ) {
3725 TDirectory *backup = gDirectory;
3726 fDebugStreamer = new TTreeSRedirector("TRDDebugFitPRF.root");
3727 if ( backup ) backup->cd(); //we don't want to be cd'd to the debug streamer
3730 Int_t detector = fCountDet;
3731 Int_t layer = GetLayer(fCountDet);
3732 Int_t caligroup = idect;
3733 Short_t rowmin = fCalibraMode->GetRowMin(2);
3734 Short_t rowmax = fCalibraMode->GetRowMax(2);
3735 Short_t colmin = fCalibraMode->GetColMin(2);
3736 Short_t colmax = fCalibraMode->GetColMax(2);
3737 Float_t widf = fCurrentCoef[0];
3738 Float_t wids = fCurrentCoef[1];
3739 Float_t widfE = fCurrentCoefE;
3741 (* fDebugStreamer) << "FillFillPRF"<<
3742 "detector="<<detector<<
3744 "caligroup="<<caligroup<<
3755 for (Int_t k = 0; k < 2304; k++) {
3756 fCurrentCoefDetector[k] = 0.0;
3760 //printf("Check the count now: fCountDet %d\n",fCountDet);
3765 if(GetStack(fCountDet) == 2) factor = 12;
3768 // Pointer to the branch
3769 for (Int_t k = fCalibraMode->GetRowMin(2); k < fCalibraMode->GetRowMax(2); k++) {
3770 for (Int_t j = fCalibraMode->GetColMin(2); j < fCalibraMode->GetColMax(2); j++) {
3771 fCurrentCoefDetector[(Int_t)(j*factor+k)] = fCurrentCoef[0];
3781 //____________Functions for initialising the AliTRDCalibraFit in the code_________
3782 Bool_t AliTRDCalibraFit::FillInfosFitLinearFitter()
3785 // Fill the coefficients found with the fits or other
3786 // methods from the Fit functions
3790 if(GetStack(fCountDet) == 2) factor = 1728;
3793 // Pointer to the branch
3794 for (Int_t k = 0; k < factor; k++) {
3795 fCurrentCoefDetector[k] = fCurrentCoef[0];
3796 fCurrentCoefDetector2[k] = fCurrentCoef2[0];
3799 FillFillLinearFitter();
3804 //________________________________________________________________________________
3805 void AliTRDCalibraFit::FillFillCH(Int_t idect)
3808 // DebugStream and fVectorFit
3811 // End of one detector
3812 if ((idect == (fCount-1))) {
3815 for (Int_t k = 0; k < 2304; k++) {
3816 fCurrentCoefDetector[k] = 0.0;
3820 if(fDebugLevel > 1){
3822 if ( !fDebugStreamer ) {
3824 TDirectory *backup = gDirectory;
3825 fDebugStreamer = new TTreeSRedirector("TRDDebugFitCH.root");
3826 if ( backup ) backup->cd(); //we don't want to be cd'd to the debug streamer
3829 Int_t detector = fCountDet;
3830 Int_t caligroup = idect;
3831 Short_t rowmin = fCalibraMode->GetRowMin(0);
3832 Short_t rowmax = fCalibraMode->GetRowMax(0);
3833 Short_t colmin = fCalibraMode->GetColMin(0);
3834 Short_t colmax = fCalibraMode->GetColMax(0);
3835 Float_t gf = fCurrentCoef[0];
3836 Float_t gfs = fCurrentCoef[1];
3837 Float_t gfE = fCurrentCoefE;
3839 (*fDebugStreamer) << "FillFillCH" <<
3840 "detector=" << detector <<
3841 "caligroup=" << caligroup <<
3842 "rowmin=" << rowmin <<
3843 "rowmax=" << rowmax <<
3844 "colmin=" << colmin <<
3845 "colmax=" << colmax <<
3853 //________________________________________________________________________________
3854 void AliTRDCalibraFit::FillFillPH(Int_t idect,Double_t nentries)
3857 // DebugStream and fVectorFit and fVectorFit2
3860 // End of one detector
3861 if ((idect == (fCount-1))) {
3865 for (Int_t k = 0; k < 2304; k++) {
3866 fCurrentCoefDetector[k] = 0.0;
3867 fCurrentCoefDetector2[k] = 0.0;
3871 if(fDebugLevel > 1){
3873 if ( !fDebugStreamer ) {
3875 TDirectory *backup = gDirectory;
3876 fDebugStreamer = new TTreeSRedirector("TRDDebugFitPH.root");
3877 if ( backup ) backup->cd(); //we don't want to be cd'd to the debug streamer
3881 Int_t detector = fCountDet;
3882 Int_t caligroup = idect;
3883 Short_t rowmin = fCalibraMode->GetRowMin(1);
3884 Short_t rowmax = fCalibraMode->GetRowMax(1);
3885 Short_t colmin = fCalibraMode->GetColMin(1);
3886 Short_t colmax = fCalibraMode->GetColMax(1);
3887 Float_t vf = fCurrentCoef[0];
3888 Float_t vs = fCurrentCoef[1];
3889 Float_t vfE = fCurrentCoefE;
3890 Float_t t0f = fCurrentCoef2[0];
3891 Float_t t0s = fCurrentCoef2[1];
3892 Float_t t0E = fCurrentCoefE2;
3896 (* fDebugStreamer) << "FillFillPH"<<
3897 "detector="<<detector<<
3898 "nentries="<<nentries<<
3899 "caligroup="<<caligroup<<
3914 //________________________________________________________________________________
3915 void AliTRDCalibraFit::FillFillPRF(Int_t idect)
3918 // DebugStream and fVectorFit
3921 // End of one detector
3922 if ((idect == (fCount-1))) {
3925 for (Int_t k = 0; k < 2304; k++) {
3926 fCurrentCoefDetector[k] = 0.0;
3931 if(fDebugLevel > 1){
3933 if ( !fDebugStreamer ) {
3935 TDirectory *backup = gDirectory;
3936 fDebugStreamer = new TTreeSRedirector("TRDDebugFitPRF.root");
3937 if ( backup ) backup->cd(); //we don't want to be cd'd to the debug streamer
3940 Int_t detector = fCountDet;
3941 Int_t layer = GetLayer(fCountDet);
3942 Int_t caligroup = idect;
3943 Short_t rowmin = fCalibraMode->GetRowMin(2);
3944 Short_t rowmax = fCalibraMode->GetRowMax(2);
3945 Short_t colmin = fCalibraMode->GetColMin(2);
3946 Short_t colmax = fCalibraMode->GetColMax(2);
3947 Float_t widf = fCurrentCoef[0];
3948 Float_t wids = fCurrentCoef[1];
3949 Float_t widfE = fCurrentCoefE;
3951 (* fDebugStreamer) << "FillFillPRF"<<
3952 "detector="<<detector<<
3954 "caligroup="<<caligroup<<
3966 //________________________________________________________________________________
3967 void AliTRDCalibraFit::FillFillLinearFitter()
3970 // DebugStream and fVectorFit
3973 // End of one detector
3979 for (Int_t k = 0; k < 2304; k++) {
3980 fCurrentCoefDetector[k] = 0.0;
3981 fCurrentCoefDetector2[k] = 0.0;
3985 if(fDebugLevel > 1){
3987 if ( !fDebugStreamer ) {
3989 TDirectory *backup = gDirectory;
3990 fDebugStreamer = new TTreeSRedirector("TRDDebugFitLinearFitter.root");
3991 if ( backup ) backup->cd(); //we don't want to be cd'd to the debug streamer
3994 //Debug: comparaison of the different methods (okey for first time but not for iterative procedure)
3995 AliTRDpadPlane *padplane = fGeo->GetPadPlane(GetLayer(fCountDet),GetStack(fCountDet));
3996 Float_t rowmd = (padplane->GetRow0()+padplane->GetRowEnd())/2.;
3997 Float_t r = AliTRDgeometry::GetTime0(GetLayer(fCountDet));
3998 Float_t tiltangle = padplane->GetTiltingAngle();
3999 Int_t detector = fCountDet;
4000 Int_t stack = GetStack(fCountDet);
4001 Int_t layer = GetLayer(fCountDet);
4002 Float_t vf = fCurrentCoef[0];
4003 Float_t vs = fCurrentCoef[1];
4004 Float_t vfE = fCurrentCoefE;
4005 Float_t lorentzangler = fCurrentCoef2[0];
4006 Float_t elorentzangler = fCurrentCoefE2;
4007 Float_t lorentzangles = fCurrentCoef2[1];
4009 (* fDebugStreamer) << "FillFillLinearFitter"<<
4010 "detector="<<detector<<
4015 "tiltangle="<<tiltangle<<
4019 "lorentzangler="<<lorentzangler<<
4020 "Elorentzangler="<<elorentzangler<<
4021 "lorentzangles="<<lorentzangles<<
4027 //____________Calcul Coef Mean_________________________________________________
4029 //_____________________________________________________________________________
4030 Bool_t AliTRDCalibraFit::CalculT0CoefMean()
4033 // For the detector Dect calcul the mean time 0
4034 // for the calibration group idect from the choosen database
4037 fCurrentCoef2[1] = 0.0;
4038 if(fDebugLevel != 1){
4039 if(((fCalibraMode->GetNz(1) > 0) ||
4040 (fCalibraMode->GetNrphi(1) > 0)) && ((fCalibraMode->GetNz(1) != 10) && (fCalibraMode->GetNz(1) != 100))) {
4042 for (Int_t row = fCalibraMode->GetRowMin(1); row < fCalibraMode->GetRowMax(1); row++) {
4043 for (Int_t col = fCalibraMode->GetColMin(1); col < fCalibraMode->GetColMax(1); col++) {
4044 fCurrentCoef2[1] += (Float_t) (fCalROC2->GetValue(col,row)+fCalDet2->GetValue(fCountDet));
4048 fCurrentCoef2[1] = fCurrentCoef2[1] / ((fCalibraMode->GetColMax(1)-fCalibraMode->GetColMin(1))*(fCalibraMode->GetRowMax(1)-fCalibraMode->GetRowMin(1)));
4054 fCurrentCoef2[1] = fCalDet2->GetValue(fCountDet);
4058 for(Int_t row = 0; row < fGeo->GetRowMax(GetLayer(fCountDet),GetStack(fCountDet),GetSector(fCountDet)); row++){
4059 for(Int_t col = 0; col < fGeo->GetColMax(GetLayer(fCountDet)); col++){
4060 fCurrentCoef2[1] += (Float_t) (fCalROC2->GetValue(col,row)+fCalDet2->GetValue(fCountDet));
4063 fCurrentCoef2[1] = fCurrentCoef2[1] / ((fGeo->GetRowMax(GetLayer(fCountDet),GetStack(fCountDet),GetSector(fCountDet)))*(fGeo->GetColMax(GetLayer(fCountDet))));
4071 //_____________________________________________________________________________
4072 Bool_t AliTRDCalibraFit::CalculChargeCoefMean(Bool_t vrai)
4075 // For the detector Dect calcul the mean gain factor
4076 // for the calibration group idect from the choosen database
4079 fCurrentCoef[1] = 0.0;
4080 if(fDebugLevel != 1){
4081 if (((fCalibraMode->GetNz(0) > 0) ||
4082 (fCalibraMode->GetNrphi(0) > 0)) && ((fCalibraMode->GetNz(0) != 10) && (fCalibraMode->GetNz(0) != 100))) {
4083 for (Int_t row = fCalibraMode->GetRowMin(0); row < fCalibraMode->GetRowMax(0); row++) {
4084 for (Int_t col = fCalibraMode->GetColMin(0); col < fCalibraMode->GetColMax(0); col++) {
4085 fCurrentCoef[1] += (Float_t) (fCalROC->GetValue(col,row)*fCalDet->GetValue(fCountDet));
4086 if (vrai) fScaleFitFactor += (Float_t) (fCalROC->GetValue(col,row)*fCalDet->GetValue(fCountDet));
4089 fCurrentCoef[1] = fCurrentCoef[1] / ((fCalibraMode->GetColMax(0)-fCalibraMode->GetColMin(0))*(fCalibraMode->GetRowMax(0)-fCalibraMode->GetRowMin(0)));
4093 fCurrentCoef[1] = (Float_t) fCalDet->GetValue(fCountDet);
4094 if (vrai) fScaleFitFactor += ((Float_t) fCalDet->GetValue(fCountDet))*(fCalibraMode->GetColMax(0)-fCalibraMode->GetColMin(0))*(fCalibraMode->GetRowMax(0)-fCalibraMode->GetRowMin(0));
4099 //_____________________________________________________________________________
4100 Bool_t AliTRDCalibraFit::CalculPRFCoefMean()
4103 // For the detector Dect calcul the mean sigma of pad response
4104 // function for the calibration group idect from the choosen database
4107 fCurrentCoef[1] = 0.0;
4108 if(fDebugLevel != 1){
4109 for (Int_t row = fCalibraMode->GetRowMin(2); row < fCalibraMode->GetRowMax(2); row++) {
4110 for (Int_t col = fCalibraMode->GetColMin(2); col < fCalibraMode->GetColMax(2); col++) {
4111 fCurrentCoef[1] += (Float_t) fCalROC->GetValue(col,row);
4114 fCurrentCoef[1] = fCurrentCoef[1] / ((fCalibraMode->GetColMax(2)-fCalibraMode->GetColMin(2))*(fCalibraMode->GetRowMax(2)-fCalibraMode->GetRowMin(2)));
4118 //_____________________________________________________________________________
4119 Bool_t AliTRDCalibraFit::CalculVdriftCoefMean()
4122 // For the detector dect calcul the mean drift velocity for the
4123 // calibration group idect from the choosen database
4126 fCurrentCoef[1] = 0.0;
4127 if(fDebugLevel != 1){
4128 if (((fCalibraMode->GetNz(1) > 0) ||
4129 (fCalibraMode->GetNrphi(1) > 0)) && ((fCalibraMode->GetNz(1) != 10) && (fCalibraMode->GetNz(1) != 100))) {
4131 for (Int_t row = fCalibraMode->GetRowMin(1); row < fCalibraMode->GetRowMax(1); row++) {
4132 for (Int_t col = fCalibraMode->GetColMin(1); col < fCalibraMode->GetColMax(1); col++) {
4133 fCurrentCoef[1] += (Float_t) (fCalROC->GetValue(col,row)*fCalDet->GetValue(fCountDet));
4137 fCurrentCoef[1] = fCurrentCoef[1] / ((fCalibraMode->GetColMax(1)-fCalibraMode->GetColMin(1))*(fCalibraMode->GetRowMax(1)-fCalibraMode->GetRowMin(1)));
4142 fCurrentCoef[1] = (Float_t) fCalDet->GetValue(fCountDet);
4147 //_____________________________________________________________________________
4148 Bool_t AliTRDCalibraFit::CalculVdriftLorentzCoef()
4151 // For the detector fCountDet, mean drift velocity and tan lorentzangle
4154 fCurrentCoef[1] = fCalDet->GetValue(fCountDet);
4155 fCurrentCoef2[1] = fCalDet2->GetValue(fCountDet);
4159 //_____________________________________________________________________________
4160 Float_t AliTRDCalibraFit::GetPRFDefault(Int_t layer) const
4163 // Default width of the PRF if there is no database as reference
4168 //case 0: return 0.515;
4169 //case 1: return 0.502;
4170 //case 2: return 0.491;
4171 //case 3: return 0.481;
4172 //case 4: return 0.471;
4173 //case 5: return 0.463;
4174 //default: return 0.0;
4177 case 0: return 0.538429;
4178 case 1: return 0.524302;
4179 case 2: return 0.511591;
4180 case 3: return 0.500140;
4181 case 4: return 0.489821;
4182 case 5: return 0.480524;
4183 default: return 0.0;
4186 //________________________________________________________________________________
4187 void AliTRDCalibraFit::SetCalROC(Int_t i)
4190 // Set the calib object for fCountDet
4193 Float_t value = 0.0;
4195 //Get the CalDet object
4197 AliTRDcalibDB *cal = AliTRDcalibDB::Instance();
4199 AliInfo("Could not get calibDB");
4206 fCalROC->~AliTRDCalROC();
4207 new(fCalROC) AliTRDCalROC(*(cal->GetGainFactorROC(fCountDet)));
4208 }else fCalROC = new AliTRDCalROC(*(cal->GetGainFactorROC(fCountDet)));
4212 fCalROC->~AliTRDCalROC();
4213 new(fCalROC) AliTRDCalROC(*(cal->GetVdriftROC(fCountDet)));
4214 }else fCalROC = new AliTRDCalROC(*(cal->GetVdriftROC(fCountDet)));
4216 fCalROC2->~AliTRDCalROC();
4217 new(fCalROC2) AliTRDCalROC(*(cal->GetT0ROC(fCountDet)));
4218 }else fCalROC2 = new AliTRDCalROC(*(cal->GetT0ROC(fCountDet)));
4222 fCalROC->~AliTRDCalROC();
4223 new(fCalROC) AliTRDCalROC(*(cal->GetPRFROC(fCountDet)));
4224 }else fCalROC = new AliTRDCalROC(*(cal->GetPRFROC(fCountDet)));
4233 if(fCalROC) delete fCalROC;
4234 fCalROC = new AliTRDCalROC(GetLayer(fCountDet),GetStack(fCountDet));
4235 for(Int_t k = 0; k < fCalROC->GetNchannels(); k++){
4236 fCalROC->SetValue(k,1.0);
4240 if(fCalROC) delete fCalROC;
4241 if(fCalROC2) delete fCalROC2;
4242 fCalROC = new AliTRDCalROC(GetLayer(fCountDet),GetStack(fCountDet));
4243 fCalROC2 = new AliTRDCalROC(GetLayer(fCountDet),GetStack(fCountDet));
4244 for(Int_t k = 0; k < fCalROC->GetNchannels(); k++){
4245 fCalROC->SetValue(k,1.0);
4246 fCalROC2->SetValue(k,0.0);
4250 if(fCalROC) delete fCalROC;
4251 value = GetPRFDefault(GetLayer(fCountDet));
4252 fCalROC = new AliTRDCalROC(GetLayer(fCountDet),GetStack(fCountDet));
4253 for(Int_t k = 0; k < fCalROC->GetNchannels(); k++){
4254 fCalROC->SetValue(k,value);
4262 //____________Fit Methods______________________________________________________
4264 //_____________________________________________________________________________
4265 void AliTRDCalibraFit::FitPente(TH1* projPH)
4268 // Slope methode for the drift velocity
4272 const Float_t kDrWidth = AliTRDgeometry::DrThick();
4279 fCurrentCoefE = 0.0;
4280 fCurrentCoefE2 = 0.0;
4281 fCurrentCoef[0] = 0.0;
4282 fCurrentCoef2[0] = 0.0;
4283 TLine *line = new TLine();
4286 TAxis *xpph = projPH->GetXaxis();
4287 Int_t nbins = xpph->GetNbins();
4288 Double_t lowedge = xpph->GetBinLowEdge(1);
4289 Double_t upedge = xpph->GetBinUpEdge(xpph->GetNbins());
4290 Double_t widbins = (upedge - lowedge) / nbins;
4291 Double_t limit = upedge + 0.5 * widbins;
4294 // Beginning of the signal
4295 TH1D *pentea = new TH1D("pentea","pentea",projPH->GetNbinsX(),0,(Float_t) limit);
4296 for (Int_t k = 1; k < projPH->GetNbinsX(); k++) {
4297 pentea->SetBinContent(k,(Double_t) (projPH->GetBinContent(k+1) - projPH->GetBinContent(k)));
4299 binmax = (Int_t) pentea->GetMaximumBin();
4302 AliInfo("Put the binmax from 1 to 2 to enable the fit");
4304 if (binmax >= nbins) {
4307 AliInfo("Put the binmax from nbins-1 to nbins-2 to enable the fit");
4309 pentea->Fit("pol2","0MR","",TMath::Max(pentea->GetBinCenter(binmax-1),0.0),pentea->GetBinCenter(binmax+1));
4310 Float_t l3P1am = pentea->GetFunction("pol2")->GetParameter(1);
4311 Float_t l3P2am = pentea->GetFunction("pol2")->GetParameter(2);
4312 Float_t l3P1amE = pentea->GetFunction("pol2")->GetParError(1);
4313 Float_t l3P2amE = pentea->GetFunction("pol2")->GetParError(2);
4314 if (TMath::Abs(l3P2am) > 0.00000001) {
4315 fPhd[0] = -(l3P1am / (2 * l3P2am));
4318 if((TMath::Abs(l3P1am) > 0.0000000001) && (TMath::Abs(l3P2am) > 0.00000000001)){
4319 fCurrentCoefE2 = (l3P1amE/l3P1am + l3P2amE/l3P2am)*fPhd[0];
4322 // Amplification region
4325 for (Int_t kbin = 1; kbin < projPH->GetNbinsX(); kbin ++) {
4326 if (((projPH->GetBinContent(kbin+1) - projPH->GetBinContent(kbin)) <= 0.0) && (ju == 0) && (kbin > (fPhd[0]/widbins))) {
4333 AliInfo("Put the binmax from 1 to 2 to enable the fit");
4335 if (binmax >= nbins) {
4338 AliInfo("Put the binmax from nbins-1 to nbins-2 to enable the fit");
4340 projPH->Fit("pol2","0MR","",TMath::Max(projPH->GetBinCenter(binmax-1),0.0),projPH->GetBinCenter(binmax+1));
4341 Float_t l3P1amf = projPH->GetFunction("pol2")->GetParameter(1);
4342 Float_t l3P2amf = projPH->GetFunction("pol2")->GetParameter(2);
4343 Float_t l3P1amfE = projPH->GetFunction("pol2")->GetParError(1);
4344 Float_t l3P2amfE = projPH->GetFunction("pol2")->GetParError(2);
4345 if (TMath::Abs(l3P2amf) > 0.00000000001) {
4346 fPhd[1] = -(l3P1amf / (2 * l3P2amf));
4348 if((TMath::Abs(l3P1amf) > 0.0000000001) && (TMath::Abs(l3P2amf) > 0.000000001)){
4349 fCurrentCoefE = (l3P1amfE/l3P1amf + l3P2amfE/l3P2amf)*fPhd[1];
4352 fCurrentCoefE2 = fCurrentCoefE;
4355 TH1D *pente = new TH1D("pente","pente",projPH->GetNbinsX(),0,(Float_t) limit);
4356 for (Int_t k = TMath::Min(binmax+4,projPH->GetNbinsX()); k < projPH->GetNbinsX(); k++) {
4357 pente->SetBinContent(k,(Double_t) (projPH->GetBinContent(k+1) - projPH->GetBinContent(k)));
4360 if(pente->GetEntries() > 0) binmin = (Int_t) pente->GetMinimumBin();
4363 AliInfo("Put the binmax from 1 to 2 to enable the fit");
4365 if (binmin >= nbins) {
4368 AliInfo("Put the binmax from nbins-1 to nbins-2 to enable the fit");
4373 ,TMath::Max(pente->GetBinCenter(binmin-1), 0.0)
4374 ,TMath::Min(pente->GetBinCenter(binmin+1),(Double_t) limit));
4375 Float_t l3P1dr = pente->GetFunction("pol2")->GetParameter(1);
4376 Float_t l3P2dr = pente->GetFunction("pol2")->GetParameter(2);
4377 Float_t l3P1drE = pente->GetFunction("pol2")->GetParError(1);
4378 Float_t l3P2drE = pente->GetFunction("pol2")->GetParError(2);
4379 if (TMath::Abs(l3P2dr) > 0.00000001) {
4380 fPhd[2] = -(l3P1dr / (2 * l3P2dr));
4382 if((TMath::Abs(l3P1dr) > 0.0000000001) && (TMath::Abs(l3P2dr) > 0.00000000001)){
4383 fCurrentCoefE += (l3P1drE/l3P1dr + l3P2drE/l3P2dr)*fPhd[2];
4385 Float_t fPhdt0 = 0.0;
4386 Float_t t0Shift = 0.0;
4389 t0Shift = fT0Shift1;
4393 t0Shift = fT0Shift0;
4396 if ((fPhd[2] > fPhd[0]) &&
4397 (fPhd[2] > fPhd[1]) &&
4398 (fPhd[1] > fPhd[0]) &&
4400 fCurrentCoef[0] = (kDrWidth) / (fPhd[2]-fPhd[1]);
4401 fNumberFitSuccess++;
4403 if (fPhdt0 >= 0.0) {
4404 fCurrentCoef2[0] = (fPhdt0 - t0Shift) / widbins;
4405 if (fCurrentCoef2[0] < -1.0) {
4406 fCurrentCoef2[0] = fCurrentCoef2[1] + 100.0;
4410 fCurrentCoef2[0] = fCurrentCoef2[1] + 100.0;
4415 fCurrentCoef[0] = -TMath::Abs(fCurrentCoef[1]);
4416 fCurrentCoef2[0] = fCurrentCoef2[1] + 100.0;
4419 if (fDebugLevel == 1) {
4420 TCanvas *cpentei = new TCanvas("cpentei","cpentei",50,50,600,800);
4423 line->SetLineColor(2);
4424 line->DrawLine(fPhd[0],0,fPhd[0],projPH->GetMaximum());
4425 line->DrawLine(fPhd[1],0,fPhd[1],projPH->GetMaximum());
4426 line->DrawLine(fPhd[2],0,fPhd[2],projPH->GetMaximum());
4427 AliInfo(Form("fPhd[0] (beginning of the signal): %f" ,(Float_t) fPhd[0]));
4428 AliInfo(Form("fPhd[1] (end of the amplification region): %f" ,(Float_t) fPhd[1]));
4429 AliInfo(Form("fPhd[2] (end of the drift region): %f" ,(Float_t) fPhd[2]));
4430 AliInfo(Form("fVriftCoef[1] (with only the drift region(default)): %f",(Float_t) fCurrentCoef[0]));
4431 TCanvas *cpentei2 = new TCanvas("cpentei2","cpentei2",50,50,600,800);
4434 TCanvas *cpentei3 = new TCanvas("cpentei3","cpentei3",50,50,600,800);
4443 //_____________________________________________________________________________
4444 void AliTRDCalibraFit::FitLagrangePoly(TH1* projPH)
4447 // Slope methode but with polynomes de Lagrange
4451 const Float_t kDrWidth = AliTRDgeometry::DrThick();
4454 //Double_t *x = new Double_t[5];
4455 //Double_t *y = new Double_t[5];
4472 fCurrentCoefE = 0.0;
4473 fCurrentCoefE2 = 1.0;
4474 fCurrentCoef[0] = 0.0;
4475 fCurrentCoef2[0] = 0.0;
4476 TLine *line = new TLine();
4477 TF1 * polynome = 0x0;
4478 TF1 * polynomea = 0x0;
4479 TF1 * polynomeb = 0x0;
4487 TAxis *xpph = projPH->GetXaxis();
4488 Int_t nbins = xpph->GetNbins();
4489 Double_t lowedge = xpph->GetBinLowEdge(1);
4490 Double_t upedge = xpph->GetBinUpEdge(xpph->GetNbins());
4491 Double_t widbins = (upedge - lowedge) / nbins;
4492 Double_t limit = upedge + 0.5 * widbins;
4497 // Beginning of the signal
4498 TH1D *pentea = new TH1D("pentea","pentea",projPH->GetNbinsX(),0,(Float_t) limit);
4499 for (Int_t k = 1; k < projPH->GetNbinsX(); k++) {
4500 pentea->SetBinContent(k,(Double_t) (projPH->GetBinContent(k+1) - projPH->GetBinContent(k)));
4503 binmax = (Int_t) pentea->GetMaximumBin();
4505 Double_t minnn = 0.0;
4506 Double_t maxxx = 0.0;
4508 Int_t kase = nbins-binmax;
4516 minnn = pentea->GetBinCenter(binmax-2);
4517 maxxx = pentea->GetBinCenter(binmax);
4518 x[0] = pentea->GetBinCenter(binmax-2);
4519 x[1] = pentea->GetBinCenter(binmax-1);
4520 x[2] = pentea->GetBinCenter(binmax);
4521 y[0] = pentea->GetBinContent(binmax-2);
4522 y[1] = pentea->GetBinContent(binmax-1);
4523 y[2] = pentea->GetBinContent(binmax);
4524 CalculPolynomeLagrange2(x,y,c0,c1,c2,c3,c4);
4525 AliInfo("At the limit for beginning!");
4528 minnn = pentea->GetBinCenter(binmax-2);
4529 maxxx = pentea->GetBinCenter(binmax+1);
4530 x[0] = pentea->GetBinCenter(binmax-2);
4531 x[1] = pentea->GetBinCenter(binmax-1);
4532 x[2] = pentea->GetBinCenter(binmax);
4533 x[3] = pentea->GetBinCenter(binmax+1);
4534 y[0] = pentea->GetBinContent(binmax-2);
4535 y[1] = pentea->GetBinContent(binmax-1);
4536 y[2] = pentea->GetBinContent(binmax);
4537 y[3] = pentea->GetBinContent(binmax+1);
4538 CalculPolynomeLagrange3(x,y,c0,c1,c2,c3,c4);
4546 minnn = pentea->GetBinCenter(binmax);
4547 maxxx = pentea->GetBinCenter(binmax+2);
4548 x[0] = pentea->GetBinCenter(binmax);
4549 x[1] = pentea->GetBinCenter(binmax+1);
4550 x[2] = pentea->GetBinCenter(binmax+2);
4551 y[0] = pentea->GetBinContent(binmax);
4552 y[1] = pentea->GetBinContent(binmax+1);
4553 y[2] = pentea->GetBinContent(binmax+2);
4554 CalculPolynomeLagrange2(x,y,c0,c1,c2,c3,c4);
4557 minnn = pentea->GetBinCenter(binmax-1);
4558 maxxx = pentea->GetBinCenter(binmax+2);
4559 x[0] = pentea->GetBinCenter(binmax-1);
4560 x[1] = pentea->GetBinCenter(binmax);
4561 x[2] = pentea->GetBinCenter(binmax+1);
4562 x[3] = pentea->GetBinCenter(binmax+2);
4563 y[0] = pentea->GetBinContent(binmax-1);
4564 y[1] = pentea->GetBinContent(binmax);
4565 y[2] = pentea->GetBinContent(binmax+1);
4566 y[3] = pentea->GetBinContent(binmax+2);
4567 CalculPolynomeLagrange3(x,y,c0,c1,c2,c3,c4);
4570 minnn = pentea->GetBinCenter(binmax-2);
4571 maxxx = pentea->GetBinCenter(binmax+2);
4572 x[0] = pentea->GetBinCenter(binmax-2);
4573 x[1] = pentea->GetBinCenter(binmax-1);
4574 x[2] = pentea->GetBinCenter(binmax);
4575 x[3] = pentea->GetBinCenter(binmax+1);
4576 x[4] = pentea->GetBinCenter(binmax+2);
4577 y[0] = pentea->GetBinContent(binmax-2);
4578 y[1] = pentea->GetBinContent(binmax-1);
4579 y[2] = pentea->GetBinContent(binmax);
4580 y[3] = pentea->GetBinContent(binmax+1);
4581 y[4] = pentea->GetBinContent(binmax+2);
4582 CalculPolynomeLagrange4(x,y,c0,c1,c2,c3,c4);
4590 polynomeb = new TF1("polb","[0]+[1]*x+[2]*x*x+[3]*x*x*x+[4]*x*x*x*x",minnn,maxxx);
4591 polynomeb->SetParameters(c0,c1,c2,c3,c4);
4593 Double_t step = (maxxx-minnn)/10000;
4595 Double_t maxvalue = 0.0;
4596 Double_t placemaximum = minnn;
4597 for(Int_t o = 0; o < 10000; o++){
4598 if(o == 0) maxvalue = polynomeb->Eval(l);
4599 if(maxvalue < (polynomeb->Eval(l))){
4600 maxvalue = polynomeb->Eval(l);
4605 fPhd[0] = placemaximum;
4608 // Amplification region
4611 for (Int_t kbin = 1; kbin < projPH->GetNbinsX(); kbin ++) {
4612 if (((projPH->GetBinContent(kbin+1) - projPH->GetBinContent(kbin)) <= 0.0) && (ju == 0) && (kbin > (fPhd[0]/widbins))) {
4618 Double_t minn = 0.0;
4619 Double_t maxx = 0.0;
4631 Int_t kase1 = nbins - binmax;
4633 //Determination of minn and maxx
4634 //case binmax = nbins
4639 minn = projPH->GetBinCenter(binmax-2);
4640 maxx = projPH->GetBinCenter(binmax);
4641 x[0] = projPH->GetBinCenter(binmax-2);
4642 x[1] = projPH->GetBinCenter(binmax-1);
4643 x[2] = projPH->GetBinCenter(binmax);
4644 y[0] = projPH->GetBinContent(binmax-2);
4645 y[1] = projPH->GetBinContent(binmax-1);
4646 y[2] = projPH->GetBinContent(binmax);
4647 CalculPolynomeLagrange2(x,y,c0,c1,c2,c3,c4);
4648 //AliInfo("At the limit for the drift!");
4651 minn = projPH->GetBinCenter(binmax-2);
4652 maxx = projPH->GetBinCenter(binmax+1);
4653 x[0] = projPH->GetBinCenter(binmax-2);
4654 x[1] = projPH->GetBinCenter(binmax-1);
4655 x[2] = projPH->GetBinCenter(binmax);
4656 x[3] = projPH->GetBinCenter(binmax+1);
4657 y[0] = projPH->GetBinContent(binmax-2);
4658 y[1] = projPH->GetBinContent(binmax-1);
4659 y[2] = projPH->GetBinContent(binmax);
4660 y[3] = projPH->GetBinContent(binmax+1);
4661 CalculPolynomeLagrange3(x,y,c0,c1,c2,c3,c4);
4670 minn = projPH->GetBinCenter(binmax);
4671 maxx = projPH->GetBinCenter(binmax+2);
4672 x[0] = projPH->GetBinCenter(binmax);
4673 x[1] = projPH->GetBinCenter(binmax+1);
4674 x[2] = projPH->GetBinCenter(binmax+2);
4675 y[0] = projPH->GetBinContent(binmax);
4676 y[1] = projPH->GetBinContent(binmax+1);
4677 y[2] = projPH->GetBinContent(binmax+2);
4678 CalculPolynomeLagrange2(x,y,c0,c1,c2,c3,c4);
4681 minn = projPH->GetBinCenter(binmax-1);
4682 maxx = projPH->GetBinCenter(binmax+2);
4683 x[0] = projPH->GetBinCenter(binmax-1);
4684 x[1] = projPH->GetBinCenter(binmax);
4685 x[2] = projPH->GetBinCenter(binmax+1);
4686 x[3] = projPH->GetBinCenter(binmax+2);
4687 y[0] = projPH->GetBinContent(binmax-1);
4688 y[1] = projPH->GetBinContent(binmax);
4689 y[2] = projPH->GetBinContent(binmax+1);
4690 y[3] = projPH->GetBinContent(binmax+2);
4691 CalculPolynomeLagrange3(x,y,c0,c1,c2,c3,c4);
4694 minn = projPH->GetBinCenter(binmax-2);
4695 maxx = projPH->GetBinCenter(binmax+2);
4696 x[0] = projPH->GetBinCenter(binmax-2);
4697 x[1] = projPH->GetBinCenter(binmax-1);
4698 x[2] = projPH->GetBinCenter(binmax);
4699 x[3] = projPH->GetBinCenter(binmax+1);
4700 x[4] = projPH->GetBinCenter(binmax+2);
4701 y[0] = projPH->GetBinContent(binmax-2);
4702 y[1] = projPH->GetBinContent(binmax-1);
4703 y[2] = projPH->GetBinContent(binmax);
4704 y[3] = projPH->GetBinContent(binmax+1);
4705 y[4] = projPH->GetBinContent(binmax+2);
4706 CalculPolynomeLagrange4(x,y,c0,c1,c2,c3,c4);
4713 polynomea = new TF1("pola","[0]+[1]*x+[2]*x*x+[3]*x*x*x+[4]*x*x*x*x",minn,maxx);
4714 polynomea->SetParameters(c0,c1,c2,c3,c4);
4716 Double_t step = (maxx-minn)/1000;
4718 Double_t maxvalue = 0.0;
4719 Double_t placemaximum = minn;
4720 for(Int_t o = 0; o < 1000; o++){
4721 if(o == 0) maxvalue = polynomea->Eval(l);
4722 if(maxvalue < (polynomea->Eval(l))){
4723 maxvalue = polynomea->Eval(l);
4728 fPhd[1] = placemaximum;
4732 TH1D *pente = new TH1D("pente","pente", projPH->GetNbinsX(),0,(Float_t) limit);
4733 for (Int_t k = TMath::Min(binmax+4, projPH->GetNbinsX()); k < projPH->GetNbinsX(); k++) {
4734 pente->SetBinContent(k,(Double_t) (projPH->GetBinContent(k+1) - projPH->GetBinContent(k)));
4737 if(pente->GetEntries() > 0) binmin = (Int_t) pente->GetMinimumBin();
4743 AliInfo("Put the binmax from 1 to 2 to enable the fit");
4747 if((projPH->GetBinContent(binmin)-projPH->GetBinError(binmin)) < (projPH->GetBinContent(binmin+1))) {
4748 AliInfo("Too many fluctuations at the end!");
4751 if((projPH->GetBinContent(binmin)+projPH->GetBinError(binmin)) > (projPH->GetBinContent(binmin-1))) {
4752 AliInfo("Too many fluctuations at the end!");
4755 if(TMath::Abs(pente->GetBinContent(binmin+1)) <= 0.0000000000001){
4756 AliInfo("No entries for the next bin!");
4757 pente->SetBinContent(binmin,0);
4758 if(pente->GetEntries() > 0) binmin = (Int_t) pente->GetMinimumBin();
4774 Bool_t case1 = kFALSE;
4775 Bool_t case2 = kFALSE;
4776 Bool_t case4 = kFALSE;
4778 //Determination of min and max
4779 //case binmin <= nbins-3
4781 if((binmin <= (nbins-3)) && ((binmin-2) >= TMath::Min(binmax+4, projPH->GetNbinsX()))){
4782 min = pente->GetBinCenter(binmin-2);
4783 max = pente->GetBinCenter(binmin+2);
4784 x[0] = pente->GetBinCenter(binmin-2);
4785 x[1] = pente->GetBinCenter(binmin-1);
4786 x[2] = pente->GetBinCenter(binmin);
4787 x[3] = pente->GetBinCenter(binmin+1);
4788 x[4] = pente->GetBinCenter(binmin+2);
4789 y[0] = pente->GetBinContent(binmin-2);
4790 y[1] = pente->GetBinContent(binmin-1);
4791 y[2] = pente->GetBinContent(binmin);
4792 y[3] = pente->GetBinContent(binmin+1);
4793 y[4] = pente->GetBinContent(binmin+2);
4794 //Calcul the polynome de Lagrange
4795 CalculPolynomeLagrange4(x,y,c0,c1,c2,c3,c4);
4797 if((pente->GetBinContent(binmin+2) <= pente->GetBinContent(binmin+1)) &&
4798 (pente->GetBinContent(binmin-2) <= pente->GetBinContent(binmin-1))) {
4799 //AliInfo("polynome 4 false 1");
4802 if(((binmin+3) <= (nbins-1)) &&
4803 (pente->GetBinContent(binmin+3) <= pente->GetBinContent(binmin+2)) &&
4804 ((binmin-3) >= TMath::Min(binmax+4, projPH->GetNbinsX())) &&
4805 (pente->GetBinContent(binmin-3) <= pente->GetBinContent(binmin-2))) {
4806 AliInfo("polynome 4 false 2");
4810 if((pente->GetBinContent(binmin+2) <= pente->GetBinContent(binmin+1)) &&
4811 (pente->GetBinContent(binmin-2) > pente->GetBinContent(binmin-1))) {
4812 //AliInfo("polynome 4 case 1");
4815 if((pente->GetBinContent(binmin+2) > pente->GetBinContent(binmin+1)) &&
4816 (pente->GetBinContent(binmin-2) <= pente->GetBinContent(binmin-1))) {
4817 //AliInfo("polynome 4 case 4");
4822 //case binmin = nbins-2
4824 if(((binmin == (nbins-2)) && ((binmin-2) >= TMath::Min(binmax+4, projPH->GetNbinsX()))) ||
4826 min = pente->GetBinCenter(binmin-2);
4827 max = pente->GetBinCenter(binmin+1);
4828 x[0] = pente->GetBinCenter(binmin-2);
4829 x[1] = pente->GetBinCenter(binmin-1);
4830 x[2] = pente->GetBinCenter(binmin);
4831 x[3] = pente->GetBinCenter(binmin+1);
4832 y[0] = pente->GetBinContent(binmin-2);
4833 y[1] = pente->GetBinContent(binmin-1);
4834 y[2] = pente->GetBinContent(binmin);
4835 y[3] = pente->GetBinContent(binmin+1);
4836 //Calcul the polynome de Lagrange
4837 CalculPolynomeLagrange3(x,y,c0,c1,c2,c3,c4);
4838 //richtung +: nothing
4840 if((pente->GetBinContent(binmin-2) <= pente->GetBinContent(binmin-1))) {
4841 //AliInfo("polynome 3- case 2");
4846 if(((binmin <= (nbins-3)) && ((binmin-1) == TMath::Min(binmax+4, projPH->GetNbinsX()))) ||
4848 min = pente->GetBinCenter(binmin-1);
4849 max = pente->GetBinCenter(binmin+2);
4850 x[0] = pente->GetBinCenter(binmin-1);
4851 x[1] = pente->GetBinCenter(binmin);
4852 x[2] = pente->GetBinCenter(binmin+1);
4853 x[3] = pente->GetBinCenter(binmin+2);
4854 y[0] = pente->GetBinContent(binmin-1);
4855 y[1] = pente->GetBinContent(binmin);
4856 y[2] = pente->GetBinContent(binmin+1);
4857 y[3] = pente->GetBinContent(binmin+2);
4858 //Calcul the polynome de Lagrange
4859 CalculPolynomeLagrange3(x,y,c0,c1,c2,c3,c4);
4861 if((pente->GetBinContent(binmin+2) <= pente->GetBinContent(binmin+1))) {
4862 //AliInfo("polynome 3+ case 2");
4867 if((binmin <= (nbins-3)) && (binmin == TMath::Min(binmax+4, projPH->GetNbinsX()))){
4868 min = pente->GetBinCenter(binmin);
4869 max = pente->GetBinCenter(binmin+2);
4870 x[0] = pente->GetBinCenter(binmin);
4871 x[1] = pente->GetBinCenter(binmin+1);
4872 x[2] = pente->GetBinCenter(binmin+2);
4873 y[0] = pente->GetBinContent(binmin);
4874 y[1] = pente->GetBinContent(binmin+1);
4875 y[2] = pente->GetBinContent(binmin+2);
4876 //Calcul the polynome de Lagrange
4877 CalculPolynomeLagrange2(x,y,c0,c1,c2,c3,c4);
4879 if((pente->GetBinContent(binmin+2) <= pente->GetBinContent(binmin+1))) {
4880 //AliInfo("polynome 2+ false");
4885 if(((binmin == (nbins-2)) && ((binmin-1) == TMath::Min(binmax+4, projPH->GetNbinsX()))) ||
4887 min = pente->GetBinCenter(binmin-1);
4888 max = pente->GetBinCenter(binmin+1);
4889 x[0] = pente->GetBinCenter(binmin-1);
4890 x[1] = pente->GetBinCenter(binmin);
4891 x[2] = pente->GetBinCenter(binmin+1);
4892 y[0] = pente->GetBinContent(binmin-1);
4893 y[1] = pente->GetBinContent(binmin);
4894 y[2] = pente->GetBinContent(binmin+1);
4895 //Calcul the polynome de Lagrange
4896 CalculPolynomeLagrange2(x,y,c0,c1,c2,c3,c4);
4897 //richtung +: nothing
4898 //richtung -: nothing
4900 //case binmin = nbins-1
4902 if((binmin == (nbins-1)) && ((binmin-2) >= TMath::Min(binmax+4, projPH->GetNbinsX()))){
4903 min = pente->GetBinCenter(binmin-2);
4904 max = pente->GetBinCenter(binmin);
4905 x[0] = pente->GetBinCenter(binmin-2);
4906 x[1] = pente->GetBinCenter(binmin-1);
4907 x[2] = pente->GetBinCenter(binmin);
4908 y[0] = pente->GetBinContent(binmin-2);
4909 y[1] = pente->GetBinContent(binmin-1);
4910 y[2] = pente->GetBinContent(binmin);
4911 //Calcul the polynome de Lagrange
4912 CalculPolynomeLagrange2(x,y,c0,c1,c2,c3,c4);
4913 //AliInfo("At the limit for the drift!");
4914 //fluctuation too big!
4915 //richtung +: nothing
4917 if((pente->GetBinContent(binmin-2) <= pente->GetBinContent(binmin-1))) {
4918 //AliInfo("polynome 2- false ");
4922 if((binmin == (nbins-1)) && ((binmin-2) < TMath::Min(binmax+4, projPH->GetNbinsX()))) {
4924 AliInfo("At the limit for the drift and not usable!");
4928 if((binmin == (nbins-2)) && ((binmin-1) < TMath::Min(binmax+4, projPH->GetNbinsX()))){
4930 AliInfo("For the drift...problem!");
4932 //pass but should not happen
4933 if((binmin <= (nbins-3)) && (binmin < TMath::Min(binmax+6, projPH->GetNbinsX()))){
4935 AliInfo("For the drift...problem!");
4939 polynome = new TF1("pol","[0]+[1]*x+[2]*x*x+[3]*x*x*x+[4]*x*x*x*x",min,max);
4940 polynome->SetParameters(c0,c1,c2,c3,c4);
4941 //AliInfo(Form("GetMinimum of the function %f",polynome->GetMinimumX()));
4942 Double_t step = (max-min)/1000;
4944 Double_t minvalue = 0.0;
4945 Double_t placeminimum = min;
4946 for(Int_t o = 0; o < 1000; o++){
4947 if(o == 0) minvalue = polynome->Eval(l);
4948 if(minvalue > (polynome->Eval(l))){
4949 minvalue = polynome->Eval(l);
4954 fPhd[2] = placeminimum;
4956 //printf("La fin %d\n",((Int_t)(fPhd[2]*10.0))+2);
4957 if((((Int_t)(fPhd[2]*10.0))+2) >= projPH->GetNbinsX()) fPhd[2] = 0.0;
4958 if(((((Int_t)(fPhd[2]*10.0))+2) < projPH->GetNbinsX()) && (projPH->GetBinContent(((Int_t)(fPhd[2]*10.0))+2)==0)) fPhd[2] = 0.0;
4960 Float_t fPhdt0 = 0.0;
4961 Float_t t0Shift = 0.0;
4964 t0Shift = fT0Shift1;
4968 t0Shift = fT0Shift0;
4971 if ((fPhd[2] > fPhd[0]) &&
4972 (fPhd[2] > fPhd[1]) &&
4973 (fPhd[1] > fPhd[0]) &&
4975 fCurrentCoef[0] = (kDrWidth) / (fPhd[2]-fPhd[1]);
4976 if(fCurrentCoef[0] > 2.5) fCurrentCoef[0] = -TMath::Abs(fCurrentCoef[1]);
4977 else fNumberFitSuccess++;
4978 if (fPhdt0 >= 0.0) {
4979 fCurrentCoef2[0] = (fPhdt0 - t0Shift) / widbins;
4980 if (fCurrentCoef2[0] < -1.0) {
4981 fCurrentCoef2[0] = fCurrentCoef2[1] + 100.0;
4985 fCurrentCoef2[0] = fCurrentCoef2[1] + 100.0;
4989 //printf("Put default %f\n",-TMath::Abs(fCurrentCoef[1]));
4990 fCurrentCoef[0] = -TMath::Abs(fCurrentCoef[1]);
4992 if((fPhd[1] > fPhd[0]) &&
4994 if (fPhdt0 >= 0.0) {
4995 fCurrentCoef2[0] = (fPhdt0 - t0Shift) / widbins;
4996 if (fCurrentCoef2[0] < -1.0) {
4997 fCurrentCoef2[0] = fCurrentCoef2[1] + 100.0;
5001 fCurrentCoef2[0] = fCurrentCoef2[1] + 100.0;
5005 fCurrentCoef2[0] = fCurrentCoef2[1] + 100.0;
5006 //printf("Fit failed!\n");
5010 if (fDebugLevel == 1) {
5011 TCanvas *cpentei = new TCanvas("cpentei","cpentei",50,50,600,800);
5014 line->SetLineColor(2);
5015 line->DrawLine(fPhd[0],0,fPhd[0],projPH->GetMaximum());
5016 line->DrawLine(fPhd[1],0,fPhd[1],projPH->GetMaximum());
5017 line->DrawLine(fPhd[2],0,fPhd[2],projPH->GetMaximum());
5018 AliInfo(Form("fPhd[0] (beginning of the signal): %f" ,(Float_t) fPhd[0]));
5019 AliInfo(Form("fPhd[1] (end of the amplification region): %f" ,(Float_t) fPhd[1]));
5020 AliInfo(Form("fPhd[2] (end of the drift region): %f" ,(Float_t) fPhd[2]));
5021 AliInfo(Form("fVriftCoef[3] (with only the drift region(default)): %f",(Float_t) fCurrentCoef[0]));
5022 TCanvas *cpentei2 = new TCanvas("cpentei2","cpentei2",50,50,600,800);
5025 TCanvas *cpentei3 = new TCanvas("cpentei3","cpentei3",50,50,600,800);
5032 if(polynome) delete polynome;
5033 if(polynomea) delete polynomea;
5034 if(polynomeb) delete polynomeb;
5035 //if(x) delete [] x;
5036 //if(y) delete [] y;
5037 if(line) delete line;
5042 //if(fCurrentCoef[0] > 1.7) fCurrentCoef[0] = -TMath::Abs(fCurrentCoef[1]);
5043 //if((fCurrentCoef2[0] > 2.6) || (fCurrentCoef2[0] < 2.1)) fCurrentCoef2[0] = fCurrentCoef2[1] + 100.0;
5045 projPH->SetDirectory(0);
5049 //_____________________________________________________________________________
5050 void AliTRDCalibraFit::FitPH(TH1* projPH, Int_t idect)
5053 // Fit methode for the drift velocity
5057 const Float_t kDrWidth = AliTRDgeometry::DrThick();
5060 TAxis *xpph = projPH->GetXaxis();
5061 Double_t upedge = xpph->GetBinUpEdge(xpph->GetNbins());
5063 TF1 *fPH = new TF1("fPH",AliTRDCalibraFit::PH,-0.05,3.2,6);
5064 fPH->SetParameter(0,0.469); // Scaling
5065 fPH->SetParameter(1,0.18); // Start
5066 fPH->SetParameter(2,0.0857325); // AR
5067 fPH->SetParameter(3,1.89); // DR
5068 fPH->SetParameter(4,0.08); // QA/QD
5069 fPH->SetParameter(5,0.0); // Baseline
5071 TLine *line = new TLine();
5073 fCurrentCoef[0] = 0.0;
5074 fCurrentCoef2[0] = 0.0;
5075 fCurrentCoefE = 0.0;
5076 fCurrentCoefE2 = 0.0;
5078 if (idect%fFitPHPeriode == 0) {
5080 AliInfo(Form("The detector %d will be fitted",idect));
5081 fPH->SetParameter(0,(projPH->Integral()-(projPH->GetBinContent(1)*projPH->GetNbinsX())) * 0.00028); // Scaling
5082 fPH->SetParameter(1,fPhd[0] - 0.1); // Start
5083 fPH->SetParameter(2,fPhd[1] - fPhd[0]); // AR
5084 fPH->SetParameter(3,fPhd[2] - fPhd[1]); // DR
5085 fPH->SetParameter(4,0.225); // QA/QD
5086 fPH->SetParameter(5,(Float_t) projPH->GetBinContent(1));
5088 if (fDebugLevel != 1) {
5089 projPH->Fit(fPH,"0M","",0.0,upedge);
5092 TCanvas *cpente = new TCanvas("cpente","cpente",50,50,600,800);
5094 projPH->Fit(fPH,"M+","",0.0,upedge);
5096 line->SetLineColor(4);
5097 line->DrawLine(fPH->GetParameter(1)
5099 ,fPH->GetParameter(1)
5100 ,projPH->GetMaximum());
5101 line->DrawLine(fPH->GetParameter(1)+fPH->GetParameter(2)
5103 ,fPH->GetParameter(1)+fPH->GetParameter(2)
5104 ,projPH->GetMaximum());
5105 line->DrawLine(fPH->GetParameter(1)+fPH->GetParameter(2)+fPH->GetParameter(3)
5107 ,fPH->GetParameter(1)+fPH->GetParameter(2)+fPH->GetParameter(3)
5108 ,projPH->GetMaximum());
5111 if (fPH->GetParameter(3) != 0) {
5112 fNumberFitSuccess++;
5113 fCurrentCoef[0] = kDrWidth / (fPH->GetParameter(3));
5114 fCurrentCoefE = (fPH->GetParError(3)/fPH->GetParameter(3))*fCurrentCoef[0];
5115 fCurrentCoef2[0] = fPH->GetParameter(1);
5116 fCurrentCoefE2 = fPH->GetParError(1);
5119 fCurrentCoef[0] = -TMath::Abs(fCurrentCoef[1]);
5120 fCurrentCoef2[0] = fCurrentCoef2[1] + 100.0;
5126 // Put the default value
5127 fCurrentCoef[0] = -TMath::Abs(fCurrentCoef[1]);
5128 fCurrentCoef2[0] = fCurrentCoef2[1] + 100.0;
5131 if (fDebugLevel != 1) {
5136 //_____________________________________________________________________________
5137 Bool_t AliTRDCalibraFit::FitPRFGausMI(Double_t *arraye, Double_t *arraym, Double_t *arrayme, Int_t nBins, Float_t xMin, Float_t xMax)
5140 // Fit methode for the sigma of the pad response function
5145 fCurrentCoef[0] = 0.0;
5146 fCurrentCoefE = 0.0;
5148 Double_t ret = FitGausMI(arraye, arraym, arrayme, nBins, xMin, xMax,¶m);
5150 if(TMath::Abs(ret+4) <= 0.000000001){
5151 fCurrentCoef[0] = -fCurrentCoef[1];
5155 fNumberFitSuccess++;
5156 fCurrentCoef[0] = param[2];
5157 fCurrentCoefE = ret;
5161 //_____________________________________________________________________________
5162 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)
5165 // Fit methode for the sigma of the pad response function
5168 //We should have at least 3 points
5169 if(nBins <=3) return -4.0;
5171 TLinearFitter fitter(3,"pol2");
5172 fitter.StoreData(kFALSE);
5173 fitter.ClearPoints();
5175 Float_t binWidth = (xMax-xMin)/(Float_t)nBins;
5176 Float_t entries = 0;
5177 Int_t nbbinwithentries = 0;
5178 for (Int_t i=0; i<nBins; i++){
5180 if(arraye[i] > 15) nbbinwithentries++;
5181 //printf("entries for i %d: %f\n",i,arraye[i]);
5183 if ((entries<700) || (nbbinwithentries < ((Int_t)(nBins/2)))) return -4;
5184 //printf("entries %f\n",entries);
5185 //printf("nbbinwithentries %d\n",nbbinwithentries);
5188 Float_t errorm = 0.0;
5189 Float_t errorn = 0.0;
5190 Float_t error = 0.0;
5193 for (Int_t ibin=0;ibin<nBins; ibin++){
5194 Float_t entriesI = arraye[ibin];
5195 Float_t valueI = arraym[ibin];
5196 Double_t xcenter = 0.0;
5198 if ((entriesI>15) && (valueI>0.0)){
5199 xcenter = xMin+(ibin+0.5)*binWidth;
5204 if((valueI + 0.01) > 0.0) errorm = TMath::Log((valueI + 0.01)/valueI);
5205 if((valueI - 0.01) > 0.0) errorn = TMath::Log((valueI - 0.01)/valueI);
5206 error = TMath::Max(TMath::Abs(errorm),TMath::Abs(errorn));
5209 if((valueI + arrayme[ibin]) > 0.0) errorm = TMath::Log((valueI + arrayme[ibin])/valueI);
5210 if((valueI - arrayme[ibin]) > 0.0) errorn = TMath::Log((valueI - arrayme[ibin])/valueI);
5211 error = TMath::Max(TMath::Abs(errorm),TMath::Abs(errorn));
5213 if(TMath::Abs(error) < 0.000000001) continue;
5214 val = TMath::Log(Float_t(valueI));
5215 fitter.AddPoint(&xcenter,val,error);
5219 if(fDebugLevel > 1){
5221 if ( !fDebugStreamer ) {
5223 TDirectory *backup = gDirectory;
5224 fDebugStreamer = new TTreeSRedirector("TRDDebugFitPRF.root");
5225 if ( backup ) backup->cd(); //we don't want to be cd'd to the debug streamer
5228 Int_t detector = fCountDet;
5229 Int_t layer = GetLayer(fCountDet);
5232 (* fDebugStreamer) << "FitGausMIFill"<<
5233 "detector="<<detector<<
5237 "entriesI="<<entriesI<<
5240 "xcenter="<<xcenter<<
5250 if(npoints <=3) return -4.0;
5255 fitter.GetParameters(par);
5256 chi2 = fitter.GetChisquare()/Float_t(npoints);
5259 if (!param) param = new TVectorD(3);
5260 if(TMath::Abs(par[2]) <= 0.000000001) return -4.0;
5261 Double_t x = TMath::Sqrt(TMath::Abs(-2*par[2]));
5262 Double_t deltax = (fitter.GetParError(2))/x;
5263 Double_t errorparam2 = TMath::Abs(deltax)/(x*x);
5266 (*param)[1] = par[1]/(-2.*par[2]);
5267 (*param)[2] = 1./TMath::Sqrt(TMath::Abs(-2.*par[2]));
5268 Double_t lnparam0 = par[0]+ par[1]* (*param)[1] + par[2]*(*param)[1]*(*param)[1];
5269 if ( lnparam0>307 ) return -4;
5270 (*param)[0] = TMath::Exp(lnparam0);
5272 if(fDebugLevel > 1){
5274 if ( !fDebugStreamer ) {
5276 TDirectory *backup = gDirectory;
5277 fDebugStreamer = new TTreeSRedirector("TRDDebugFitPRF.root");
5278 if ( backup ) backup->cd(); //we don't want to be cd'd to the debug streamer
5281 Int_t detector = fCountDet;
5282 Int_t layer = GetLayer(fCountDet);
5285 (* fDebugStreamer) << "FitGausMIFit"<<
5286 "detector="<<detector<<
5289 "errorsigma="<<chi2<<
5290 "mean="<<(*param)[1]<<
5291 "sigma="<<(*param)[2]<<
5292 "constant="<<(*param)[0]<<
5297 if((chi2/(*param)[2]) > 0.1){
5299 chi2 = FitGausMI(arraye,arraym,arrayme,nBins,xMin,xMax,param,kFALSE);
5304 if(fDebugLevel == 1){
5305 TString name("PRF");
5306 name += (Int_t)xMin;
5307 name += (Int_t)xMax;
5308 TCanvas *c1 = new TCanvas((const char *)name,(const char *)name,50,50,600,800);
5311 TH1F *histo = new TH1F((const char *)name,(const char *)name,nBins,xMin,xMax);
5312 for(Int_t k = 0; k < nBins; k++){
5313 histo->SetBinContent(k+1,arraym[k]);
5314 histo->SetBinError(k+1,arrayme[k]);
5317 name += "functionf";
5318 TF1 *f1= new TF1((const char*)name,"[0]*exp(-(x-[1])^2/(2*[2]*[2]))",xMin,xMax);
5319 f1->SetParameter(0, (*param)[0]);
5320 f1->SetParameter(1, (*param)[1]);
5321 f1->SetParameter(2, (*param)[2]);
5329 //_____________________________________________________________________________
5330 void AliTRDCalibraFit::FitPRF(TH1 *projPRF)
5333 // Fit methode for the sigma of the pad response function
5336 fCurrentCoef[0] = 0.0;
5337 fCurrentCoefE = 0.0;
5339 if (fDebugLevel != 1) {
5340 projPRF->Fit("gaus","0M","",-fRangeFitPRF,fRangeFitPRF);
5343 TCanvas *cfit = new TCanvas("cfit","cfit",50,50,600,800);
5345 projPRF->Fit("gaus","M+","",-fRangeFitPRF,fRangeFitPRF);
5348 fCurrentCoef[0] = projPRF->GetFunction("gaus")->GetParameter(2);
5349 fCurrentCoefE = projPRF->GetFunction("gaus")->GetParError(2);
5350 if(fCurrentCoef[0] <= 0.0) fCurrentCoef[0] = -fCurrentCoef[1];
5352 fNumberFitSuccess++;
5355 //_____________________________________________________________________________
5356 void AliTRDCalibraFit::RmsPRF(TH1 *projPRF)
5359 // Fit methode for the sigma of the pad response function
5361 fCurrentCoef[0] = 0.0;
5362 fCurrentCoefE = 0.0;
5363 if (fDebugLevel == 1) {
5364 TCanvas *cfit = new TCanvas("cfit","cfit",50,50,600,800);
5368 fCurrentCoef[0] = projPRF->GetRMS();
5369 if(fCurrentCoef[0] <= 0.0) fCurrentCoef[0] = -fCurrentCoef[1];
5371 fNumberFitSuccess++;
5374 //_____________________________________________________________________________
5375 void AliTRDCalibraFit::FitTnpRange(Double_t *arraye, Double_t *arraym, Double_t *arrayme, Int_t nbg, Int_t nybins)
5378 // Fit methode for the sigma of the pad response function with 2*nbg tan bins
5381 TLinearFitter linearfitter = TLinearFitter(3,"pol2");
5384 Int_t nbins = (Int_t)(nybins/(2*nbg));
5385 Float_t lowedge = -3.0*nbg;
5386 Float_t upedge = lowedge + 3.0;
5389 Double_t xvalues = -0.2*nbg+0.1;
5391 Int_t total = 2*nbg;
5394 for(Int_t k = 0; k < total; k++){
5395 if(FitPRFGausMI(arraye+offset, arraym+offset, arrayme+offset, nbins, lowedge, upedge)){
5397 y = fCurrentCoef[0]*fCurrentCoef[0];
5398 linearfitter.AddPoint(&xvalues,y,2*fCurrentCoefE*fCurrentCoef[0]);
5401 if(fDebugLevel > 1){
5403 if ( !fDebugStreamer ) {
5405 TDirectory *backup = gDirectory;
5406 fDebugStreamer = new TTreeSRedirector("TRDDebugFitPRF.root");
5407 if ( backup ) backup->cd(); //we don't want to be cd'd to the debug streamer
5410 Int_t detector = fCountDet;
5411 Int_t layer = GetLayer(fCountDet);
5412 Int_t nbtotal = total;
5414 Float_t low = lowedge;
5415 Float_t up = upedge;
5416 Float_t tnp = xvalues;
5417 Float_t wid = fCurrentCoef[0];
5418 Float_t widfE = fCurrentCoefE;
5420 (* fDebugStreamer) << "FitTnpRange0"<<
5421 "detector="<<detector<<
5423 "nbtotal="<<nbtotal<<
5441 fCurrentCoefE = 0.0;
5442 fCurrentCoef[0] = 0.0;
5444 //printf("npoints\n",npoints);
5447 fCurrentCoef[0] = -TMath::Abs(fCurrentCoef[1]);
5452 linearfitter.Eval();
5453 linearfitter.GetParameters(pars0);
5454 Double_t pointError0 = TMath::Sqrt(linearfitter.GetChisquare()/npoints);
5455 Double_t errorsx0 = linearfitter.GetParError(2)*pointError0;
5456 Double_t min0 = 0.0;
5457 Double_t ermin0 = 0.0;
5458 //Double_t prfe0 = 0.0;
5459 Double_t prf0 = 0.0;
5460 if((pars0[2] > 0.000000000001) && (TMath::Abs(pars0[1]) >= 0.000000000001)) {
5461 min0 = -pars0[1]/(2*pars0[2]);
5462 ermin0 = TMath::Abs(min0*(errorsx0/pars0[2]+linearfitter.GetParError(1)*pointError0/pars0[1]));
5463 prf0 = pars0[0]+pars0[1]*min0+pars0[2]*min0*min0;
5466 prfe0 = linearfitter->GetParError(0)*pointError0
5467 +(linearfitter->GetParError(1)*pointError0/pars0[1]+ermin0/min0)*pars0[1]*min0
5468 +(linearfitter->GetParError(2)*pointError0/pars0[2]+2*ermin0/min0)*pars0[2]*min0*min0;
5469 prfe0 = prfe0/(2*TMath::Sqrt(prf0));
5470 fCurrentCoefE = (Float_t) prfe0;
5472 fCurrentCoef[0] = (Float_t) TMath::Sqrt(TMath::Abs(prf0));
5475 fCurrentCoef[0] = -TMath::Abs(fCurrentCoef[1]);
5479 fCurrentCoef[0] = -TMath::Abs(fCurrentCoef[1]);
5482 if(fDebugLevel > 1){
5484 if ( !fDebugStreamer ) {
5486 TDirectory *backup = gDirectory;
5487 fDebugStreamer = new TTreeSRedirector("TRDDebugFitPRF.root");
5488 if ( backup ) backup->cd(); //we don't want to be cd'd to the debug streamer
5491 Int_t detector = fCountDet;
5492 Int_t layer = GetLayer(fCountDet);
5493 Int_t nbtotal = total;
5494 Double_t colsize[6] = {0.635,0.665,0.695,0.725,0.755,0.785};
5495 Double_t sigmax = TMath::Sqrt(TMath::Abs(pars0[2]))*10000*colsize[layer];
5497 (* fDebugStreamer) << "FitTnpRange1"<<
5498 "detector="<<detector<<
5500 "nbtotal="<<nbtotal<<
5504 "npoints="<<npoints<<
5507 "sigmaprf="<<fCurrentCoef[0]<<
5508 "sigprf="<<fCurrentCoef[1]<<
5515 //_____________________________________________________________________________
5516 void AliTRDCalibraFit::FitMean(TH1 *projch, Double_t nentries, Double_t mean)
5519 // Only mean methode for the gain factor
5522 fCurrentCoef[0] = mean;
5523 fCurrentCoefE = 0.0;
5524 if(nentries > 0) fCurrentCoefE = projch->GetRMS()/TMath::Sqrt(nentries);
5525 if (fDebugLevel == 1) {
5526 TCanvas *cpmean = new TCanvas("cpmean","cpmean",50,50,600,800);
5530 CalculChargeCoefMean(kTRUE);
5531 fNumberFitSuccess++;
5533 //_____________________________________________________________________________
5534 void AliTRDCalibraFit::FitMeanW(TH1 *projch, Double_t nentries)
5537 // mean w methode for the gain factor
5541 Int_t nybins = projch->GetNbinsX();
5543 //The weight function
5544 Double_t a = 0.00228515;
5545 Double_t b = -0.00231487;
5546 Double_t c = 0.00044298;
5547 Double_t d = -0.00379239;
5548 Double_t e = 0.00338349;
5558 //A arbitrary error for the moment
5559 fCurrentCoefE = 0.0;
5560 fCurrentCoef[0] = 0.0;
5563 Double_t sumw = 0.0;
5565 Float_t sumAll = (Float_t) nentries;
5566 Int_t sumCurrent = 0;
5567 for(Int_t k = 0; k <nybins; k++){
5568 Double_t fraction = Float_t(sumCurrent)/Float_t(sumAll);
5569 if (fraction>0.95) break;
5570 Double_t weight = a + b*fraction + c*fraction*fraction + d *fraction*fraction*fraction+
5571 e*fraction*fraction*fraction*fraction;
5572 sumw += weight*projch->GetBinContent(k+1)*projch->GetBinCenter(k+1);
5573 sum += weight*projch->GetBinContent(k+1);
5574 sumCurrent += (Int_t) projch->GetBinContent(k+1);
5575 //printf("fraction %f, weight %f, bincontent %f\n",fraction,weight,projch->GetBinContent(k+1));
5577 if(sum > 0.0) fCurrentCoef[0] = (sumw/sum);
5579 if (fDebugLevel == 1) {
5580 TCanvas *cpmeanw = new TCanvas("cpmeanw","cpmeanw",50,50,600,800);
5584 fNumberFitSuccess++;
5585 CalculChargeCoefMean(kTRUE);
5587 //_____________________________________________________________________________
5588 void AliTRDCalibraFit::FitMeanWSm(TH1 *projch, Float_t sumAll)
5591 // mean w methode for the gain factor
5595 Int_t nybins = projch->GetNbinsX();
5597 //The weight function
5598 Double_t a = 0.00228515;
5599 Double_t b = -0.00231487;
5600 Double_t c = 0.00044298;
5601 Double_t d = -0.00379239;
5602 Double_t e = 0.00338349;
5612 //A arbitrary error for the moment
5613 fCurrentCoefE = 0.0;
5614 fCurrentCoef[0] = 0.0;
5617 Double_t sumw = 0.0;
5619 Int_t sumCurrent = 0;
5620 for(Int_t k = 0; k <nybins; k++){
5621 Double_t fraction = Float_t(sumCurrent)/Float_t(sumAll);
5622 if (fraction>0.95) break;
5623 Double_t weight = a + b*fraction + c*fraction*fraction + d *fraction*fraction*fraction+
5624 e*fraction*fraction*fraction*fraction;
5625 sumw += weight*projch->GetBinContent(k+1)*projch->GetBinCenter(k+1);
5626 sum += weight*projch->GetBinContent(k+1);
5627 sumCurrent += (Int_t) projch->GetBinContent(k+1);
5628 //printf("fraction %f, weight %f, bincontent %f\n",fraction,weight,projch->GetBinContent(k+1));
5630 if(sum > 0.0) fCurrentCoef[0] = (sumw/sum);
5632 if (fDebugLevel == 1) {
5633 TCanvas *cpmeanw = new TCanvas("cpmeanw","cpmeanw",50,50,600,800);
5637 fNumberFitSuccess++;
5639 //_____________________________________________________________________________
5640 void AliTRDCalibraFit::FitCH(TH1 *projch, Double_t mean)
5643 // Fit methode for the gain factor
5646 fCurrentCoef[0] = 0.0;
5647 fCurrentCoefE = 0.0;
5648 Double_t chisqrl = 0.0;
5649 Double_t chisqrg = 0.0;
5650 Double_t chisqr = 0.0;
5651 TF1 *fLandauGaus = new TF1("fLandauGaus",FuncLandauGaus,0,300,5);
5653 projch->Fit("landau","0",""
5654 ,(Double_t) mean/fBeginFitCharge
5655 ,projch->GetBinCenter(projch->GetNbinsX()));
5656 Double_t l3P0 = projch->GetFunction("landau")->GetParameter(0);
5657 Double_t l3P1 = projch->GetFunction("landau")->GetParameter(1);
5658 Double_t l3P2 = projch->GetFunction("landau")->GetParameter(2);
5659 chisqrl = projch->GetFunction("landau")->GetChisquare();
5661 projch->Fit("gaus","0",""
5662 ,(Double_t) mean/fBeginFitCharge
5663 ,projch->GetBinCenter(projch->GetNbinsX()));
5664 Double_t g3P0 = projch->GetFunction("gaus")->GetParameter(0);
5665 Double_t g3P2 = projch->GetFunction("gaus")->GetParameter(2);
5666 chisqrg = projch->GetFunction("gaus")->GetChisquare();
5668 fLandauGaus->SetParameters(l3P0,l3P1,l3P2,g3P0,g3P2);
5669 if (fDebugLevel != 1) {
5670 projch->Fit("fLandauGaus","0",""
5671 ,(Double_t) mean/fBeginFitCharge
5672 ,projch->GetBinCenter(projch->GetNbinsX()));
5673 chisqr = projch->GetFunction("fLandauGaus")->GetChisquare();
5676 TCanvas *cp = new TCanvas("cp","cp",50,50,600,800);
5678 projch->Fit("fLandauGaus","+",""
5679 ,(Double_t) mean/fBeginFitCharge
5680 ,projch->GetBinCenter(projch->GetNbinsX()));
5681 chisqr = projch->GetFunction("fLandauGaus")->GetChisquare();
5683 fLandauGaus->Draw("same");
5686 if ((projch->GetFunction("fLandauGaus")->GetParameter(1) > 0) && (projch->GetFunction("fLandauGaus")->GetParError(1) < (0.05*projch->GetFunction("fLandauGaus")->GetParameter(1))) && (chisqr < chisqrl) && (chisqr < chisqrg)) {
5687 //if ((projch->GetFunction("fLandauGaus")->GetParameter(1) > 0) && (chisqr < chisqrl) && (chisqr < chisqrg)) {
5688 fNumberFitSuccess++;
5689 CalculChargeCoefMean(kTRUE);
5690 fCurrentCoef[0] = projch->GetFunction("fLandauGaus")->GetParameter(1);
5691 fCurrentCoefE = projch->GetFunction("fLandauGaus")->GetParError(1);
5694 CalculChargeCoefMean(kFALSE);
5695 fCurrentCoef[0] = -TMath::Abs(fCurrentCoef[1]);
5698 if (fDebugLevel != 1) {
5703 //_____________________________________________________________________________
5704 void AliTRDCalibraFit::FitBisCH(TH1* projch, Double_t mean)
5707 // Fit methode for the gain factor more time consuming
5711 //Some parameters to initialise
5712 Double_t widthLandau, widthGaus, mPV, integral;
5713 Double_t chisquarel = 0.0;
5714 Double_t chisquareg = 0.0;
5715 projch->Fit("landau","0M+",""
5717 ,projch->GetBinCenter(projch->GetNbinsX()));
5718 widthLandau = projch->GetFunction("landau")->GetParameter(2);
5719 chisquarel = projch->GetFunction("landau")->GetChisquare();
5720 projch->Fit("gaus","0M+",""
5722 ,projch->GetBinCenter(projch->GetNbinsX()));
5723 widthGaus = projch->GetFunction("gaus")->GetParameter(2);
5724 chisquareg = projch->GetFunction("gaus")->GetChisquare();
5726 mPV = (projch->GetFunction("landau")->GetParameter(1))/2;
5727 integral = (projch->GetFunction("gaus")->Integral(0.3*mean,3*mean)+projch->GetFunction("landau")->Integral(0.3*mean,3*mean))/2;
5729 // Setting fit range and start values
5731 //Double_t sv[4] = { l3P2, fChargeCoef[1], projch->Integral("width"), fG3P2 };
5732 //Double_t sv[4] = { fL3P2, fChargeCoef[1], fL3P0, fG3P2 };
5733 Double_t sv[4] = { widthLandau, mPV, integral, widthGaus};
5734 Double_t pllo[4] = { 0.001, 0.001, projch->Integral()/3, 0.001};
5735 Double_t plhi[4] = { 300.0, 300.0, 30*projch->Integral(), 300.0};
5736 Double_t fp[4] = { 1.0, 1.0, 1.0, 1.0 };
5737 Double_t fpe[4] = { 1.0, 1.0, 1.0, 1.0 };
5740 fCurrentCoef[0] = 0.0;
5741 fCurrentCoefE = 0.0;
5745 TF1 *fitsnr = LanGauFit(projch,&fr[0],&sv[0]
5750 Double_t projchPeak;
5751 Double_t projchFWHM;
5752 LanGauPro(fp,projchPeak,projchFWHM);
5754 if ((fp[1] > 0) && ((fpe[1] < (0.05*fp[1])) && (chisqr < chisquarel) && (chisqr < chisquareg))) {
5755 //if ((fp[1] > 0) && ((chisqr < chisquarel) && (chisqr < chisquareg))) {
5756 fNumberFitSuccess++;
5757 CalculChargeCoefMean(kTRUE);
5758 fCurrentCoef[0] = fp[1];
5759 fCurrentCoefE = fpe[1];
5760 //chargeCoefE2 = chisqr;
5763 CalculChargeCoefMean(kFALSE);
5764 fCurrentCoef[0] = -TMath::Abs(fCurrentCoef[1]);
5766 if (fDebugLevel == 1) {
5767 AliInfo(Form("fChargeCoef[0]: %f",(Float_t) fCurrentCoef[0]));
5768 TCanvas *cpy = new TCanvas("cpy","cpy",50,50,600,800);
5771 fitsnr->Draw("same");
5777 //_____________________________________________________________________________
5778 void AliTRDCalibraFit::CalculPolynomeLagrange2(const Double_t *x, const Double_t *y, Double_t &c0, Double_t &c1, Double_t &c2, Double_t &c3, Double_t &c4) const
5781 // Calcul the coefficients of the polynome passant par ces trois points de degre 2
5783 Double_t x0 = y[0]/((x[0]-x[1])*(x[0]-x[2]));
5784 Double_t x1 = y[1]/((x[1]-x[0])*(x[1]-x[2]));
5785 Double_t x2 = y[2]/((x[2]-x[0])*(x[2]-x[1]));
5790 c1 = -(x0*(x[1]+x[2])+x1*(x[0]+x[2])+x2*(x[0]+x[1]));
5791 c0 = x0*x[1]*x[2]+x1*x[0]*x[2]+x2*x[0]*x[1];
5795 //_____________________________________________________________________________
5796 void AliTRDCalibraFit::CalculPolynomeLagrange3(const Double_t *x, const Double_t *y, Double_t &c0, Double_t &c1, Double_t &c2, Double_t &c3, Double_t &c4) const
5799 // Calcul the coefficients of the polynome passant par ces quatre points de degre 3
5801 Double_t x0 = y[0]/((x[0]-x[1])*(x[0]-x[2])*(x[0]-x[3]));
5802 Double_t x1 = y[1]/((x[1]-x[0])*(x[1]-x[2])*(x[1]-x[3]));
5803 Double_t x2 = y[2]/((x[2]-x[0])*(x[2]-x[1])*(x[2]-x[3]));
5804 Double_t x3 = y[3]/((x[3]-x[0])*(x[3]-x[1])*(x[3]-x[2]));
5808 c2 = -(x0*(x[1]+x[2]+x[3])
5809 +x1*(x[0]+x[2]+x[3])
5810 +x2*(x[0]+x[1]+x[3])
5811 +x3*(x[0]+x[1]+x[2]));
5812 c1 = (x0*(x[1]*x[2]+x[1]*x[3]+x[2]*x[3])
5813 +x1*(x[0]*x[2]+x[0]*x[3]+x[2]*x[3])
5814 +x2*(x[0]*x[1]+x[0]*x[3]+x[1]*x[3])
5815 +x3*(x[0]*x[1]+x[0]*x[2]+x[1]*x[2]));
5817 c0 = -(x0*x[1]*x[2]*x[3]
5820 +x3*x[0]*x[1]*x[2]);
5825 //_____________________________________________________________________________
5826 void AliTRDCalibraFit::CalculPolynomeLagrange4(const Double_t *x, const Double_t *y, Double_t &c0, Double_t &c1, Double_t &c2, Double_t &c3, Double_t &c4) const
5829 // Calcul the coefficients of the polynome passant par ces cinqs points de degre 4
5831 Double_t x0 = y[0]/((x[0]-x[1])*(x[0]-x[2])*(x[0]-x[3])*(x[0]-x[4]));
5832 Double_t x1 = y[1]/((x[1]-x[0])*(x[1]-x[2])*(x[1]-x[3])*(x[1]-x[4]));
5833 Double_t x2 = y[2]/((x[2]-x[0])*(x[2]-x[1])*(x[2]-x[3])*(x[2]-x[4]));
5834 Double_t x3 = y[3]/((x[3]-x[0])*(x[3]-x[1])*(x[3]-x[2])*(x[3]-x[4]));
5835 Double_t x4 = y[4]/((x[4]-x[0])*(x[4]-x[1])*(x[4]-x[2])*(x[4]-x[3]));
5838 c4 = x0+x1+x2+x3+x4;
5839 c3 = -(x0*(x[1]+x[2]+x[3]+x[4])
5840 +x1*(x[0]+x[2]+x[3]+x[4])
5841 +x2*(x[0]+x[1]+x[3]+x[4])
5842 +x3*(x[0]+x[1]+x[2]+x[4])
5843 +x4*(x[0]+x[1]+x[2]+x[3]));
5844 c2 = (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])
5845 +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])
5846 +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])
5847 +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])
5848 +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]));
5850 c1 = -(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])
5851 +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])
5852 +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])
5853 +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])
5854 +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]));
5856 c0 = (x0*x[1]*x[2]*x[3]*x[4]
5857 +x1*x[0]*x[2]*x[3]*x[4]
5858 +x2*x[0]*x[1]*x[3]*x[4]
5859 +x3*x[0]*x[1]*x[2]*x[4]
5860 +x4*x[0]*x[1]*x[2]*x[3]);
5863 //_____________________________________________________________________________
5864 void AliTRDCalibraFit::NormierungCharge()
5867 // Normalisation of the gain factor resulting for the fits
5870 // Calcul of the mean of choosen method by fFitChargeNDB
5872 //printf("total number of entries %d\n",fVectorFitCH->GetEntriesFast());
5873 for (Int_t k = 0; k < (Int_t) fVectorFit.GetEntriesFast(); k++) {
5875 Int_t detector = ((AliTRDFitInfo *) fVectorFit.At(k))->GetDetector();
5876 Float_t *coef = ((AliTRDFitInfo *) fVectorFit.At(k))->GetCoef();
5877 //printf("detector %d coef[0] %f\n",detector,coef[0]);
5878 if (GetStack(detector) == 2) {
5881 if (GetStack(detector) != 2) {
5884 for (Int_t j = 0; j < total; j++) {
5892 fScaleFitFactor = fScaleFitFactor / sum;
5895 fScaleFitFactor = 1.0;
5898 //methode de boeuf mais bon...
5899 Double_t scalefactor = fScaleFitFactor;
5901 if(fDebugLevel > 1){
5903 if ( !fDebugStreamer ) {
5905 TDirectory *backup = gDirectory;
5906 fDebugStreamer = new TTreeSRedirector("TRDDebugFitCH.root");
5907 if ( backup ) backup->cd(); //we don't want to be cd'd to the debug streamer
5909 (* fDebugStreamer) << "NormierungCharge"<<
5910 "scalefactor="<<scalefactor<<
5914 //_____________________________________________________________________________
5915 TH1I *AliTRDCalibraFit::ReBin(const TH1I *hist) const
5918 // Rebin of the 1D histo for the gain calibration if needed.
5919 // you have to choose fRebin, divider of fNumberBinCharge
5922 TAxis *xhist = hist->GetXaxis();
5923 TH1I *rehist = new TH1I("projrebin","",(Int_t) xhist->GetNbins()/fRebin
5924 ,xhist->GetBinLowEdge(1)
5925 ,xhist->GetBinUpEdge(xhist->GetNbins()));
5927 AliInfo(Form("fRebin: %d",fRebin));
5929 for (Int_t k = 1; k <= (Int_t) xhist->GetNbins()/fRebin; k++) {
5931 for (Int_t ji = i; ji < i+fRebin; ji++) {
5932 sum += hist->GetBinContent(ji);
5935 rehist->SetBinContent(k,sum);
5943 //_____________________________________________________________________________
5944 TH1F *AliTRDCalibraFit::ReBin(const TH1F *hist) const
5947 // Rebin of the 1D histo for the gain calibration if needed
5948 // you have to choose fRebin divider of fNumberBinCharge
5951 TAxis *xhist = hist->GetXaxis();
5952 TH1F *rehist = new TH1F("projrebin","",(Int_t) xhist->GetNbins()/fRebin
5953 ,xhist->GetBinLowEdge(1)
5954 ,xhist->GetBinUpEdge(xhist->GetNbins()));
5956 AliInfo(Form("fRebin: %d",fRebin));
5958 for (Int_t k = 1; k <= (Int_t) xhist->GetNbins()/fRebin; k++) {
5960 for (Int_t ji = i; ji < i+fRebin; ji++) {
5961 sum += hist->GetBinContent(ji);
5964 rehist->SetBinContent(k,sum);
5972 //____________Some basic geometry function_____________________________________
5975 //_____________________________________________________________________________
5976 Int_t AliTRDCalibraFit::GetLayer(Int_t d) const
5979 // Reconstruct the plane number from the detector number
5982 return ((Int_t) (d % 6));
5986 //_____________________________________________________________________________
5987 Int_t AliTRDCalibraFit::GetStack(Int_t d) const
5990 // Reconstruct the stack number from the detector number
5992 const Int_t kNlayer = 6;
5994 return ((Int_t) (d % 30) / kNlayer);
5998 //_____________________________________________________________________________
5999 Int_t AliTRDCalibraFit::GetSector(Int_t d) const
6002 // Reconstruct the sector number from the detector number
6006 return ((Int_t) (d / fg));
6011 //____________Fill and Init tree Gain, PRF, Vdrift and T0______________________
6013 //_______________________________________________________________________________
6014 void AliTRDCalibraFit::ResetVectorFit()
6017 // Reset the VectorFits
6020 fVectorFit.SetOwner();
6022 fVectorFit2.SetOwner();
6023 fVectorFit2.Clear();
6027 //____________Private Functions________________________________________________
6030 //_____________________________________________________________________________
6031 Double_t AliTRDCalibraFit::PH(const Double_t *x, const Double_t *par)
6034 // Function for the fit
6037 //TF1 *fAsymmGauss = new TF1("fAsymmGauss",AsymmGauss,0,4,6);
6039 //PARAMETERS FOR FIT PH
6041 //fAsymmGauss->SetParameter(0,0.113755);
6042 //fAsymmGauss->SetParameter(1,0.350706);
6043 //fAsymmGauss->SetParameter(2,0.0604244);
6044 //fAsymmGauss->SetParameter(3,7.65596);
6045 //fAsymmGauss->SetParameter(4,1.00124);
6046 //fAsymmGauss->SetParameter(5,0.870597); // No tail cancelation
6054 Double_t dx = 0.005;
6055 Double_t xs = par[1];
6057 Double_t paras[2] = { 0.0, 0.0 };
6060 if ((xs >= par[1]) &&
6061 (xs < (par[1]+par[2]))) {
6062 //fAsymmGauss->SetParameter(0,par[0]);
6063 //fAsymmGauss->SetParameter(1,xs);
6064 //ss += fAsymmGauss->Eval(xx);
6067 ss += AsymmGauss(&xx,paras);
6069 if ((xs >= (par[1]+par[2])) &&
6070 (xs < (par[1]+par[2]+par[3]))) {
6071 //fAsymmGauss->SetParameter(0,par[0]*par[4]);
6072 //fAsymmGauss->SetParameter(1,xs);
6073 //ss += fAsymmGauss->Eval(xx);
6074 paras[0] = par[0]*par[4];
6076 ss += AsymmGauss(&xx,paras);
6085 //_____________________________________________________________________________
6086 Double_t AliTRDCalibraFit::AsymmGauss(const Double_t *x, const Double_t *par)
6089 // Function for the fit
6092 //par[0] = normalization
6100 Double_t par1save = par[1];
6101 //Double_t par2save = par[2];
6102 Double_t par2save = 0.0604244;
6103 //Double_t par3save = par[3];
6104 Double_t par3save = 7.65596;
6105 //Double_t par5save = par[5];
6106 Double_t par5save = 0.870597;
6107 Double_t dx = x[0] - par1save;
6109 Double_t sigma2 = par2save*par2save;
6110 Double_t sqrt2 = TMath::Sqrt(2.0);
6111 Double_t exp1 = par3save * TMath::Exp(-par3save * (dx - 0.5 * par3save * sigma2))
6112 * (1.0 - AliMathBase::ErfFast((par3save * sigma2 - dx) / (sqrt2 * par2save)));
6113 Double_t exp2 = par5save * TMath::Exp(-par5save * (dx - 0.5 * par5save * sigma2))
6114 * (1.0 - AliMathBase::ErfFast((par5save * sigma2 - dx) / (sqrt2 * par2save)));
6116 //return par[0]*(exp1+par[4]*exp2);
6117 return par[0] * (exp1 + 1.00124 * exp2);
6121 //_____________________________________________________________________________
6122 Double_t AliTRDCalibraFit::FuncLandauGaus(const Double_t *x, const Double_t *par)
6125 // Sum Landau + Gaus with identical mean
6128 Double_t valLandau = par[0] * TMath::Landau(x[0],par[1],par[2]);
6129 //Double_t valGaus = par[3] * TMath::Gaus(x[0],par[4],par[5]);
6130 Double_t valGaus = par[3] * TMath::Gaus(x[0],par[1],par[4]);
6131 Double_t val = valLandau + valGaus;
6137 //_____________________________________________________________________________
6138 Double_t AliTRDCalibraFit::LanGauFun(const Double_t *x, const Double_t *par)
6141 // Function for the fit
6144 // par[0]=Width (scale) parameter of Landau density
6145 // par[1]=Most Probable (MP, location) parameter of Landau density
6146 // par[2]=Total area (integral -inf to inf, normalization constant)
6147 // par[3]=Width (sigma) of convoluted Gaussian function
6149 // In the Landau distribution (represented by the CERNLIB approximation),
6150 // the maximum is located at x=-0.22278298 with the location parameter=0.
6151 // This shift is corrected within this function, so that the actual
6152 // maximum is identical to the MP parameter.
6155 // Numeric constants
6156 Double_t invsq2pi = 0.3989422804014; // (2 pi)^(-1/2)
6157 Double_t mpshift = -0.22278298; // Landau maximum location
6159 // Control constants
6160 Double_t np = 100.0; // Number of convolution steps
6161 Double_t sc = 5.0; // Convolution extends to +-sc Gaussian sigmas
6173 // MP shift correction
6174 mpc = par[1] - mpshift * par[0];
6176 // Range of convolution integral
6177 xlow = x[0] - sc * par[3];
6178 xupp = x[0] + sc * par[3];
6180 step = (xupp - xlow) / np;
6182 // Convolution integral of Landau and Gaussian by sum
6183 for (i = 1.0; i <= np/2; i++) {
6185 xx = xlow + (i-.5) * step;
6186 fland = TMath::Landau(xx,mpc,par[0]) / par[0];
6187 sum += fland * TMath::Gaus(x[0],xx,par[3]);
6189 xx = xupp - (i-.5) * step;
6190 fland = TMath::Landau(xx,mpc,par[0]) / par[0];
6191 sum += fland * TMath::Gaus(x[0],xx,par[3]);
6195 return (par[2] * step * sum * invsq2pi / par[3]);
6198 //_____________________________________________________________________________
6199 TF1 *AliTRDCalibraFit::LanGauFit(TH1 *his, const Double_t *fitrange, const Double_t *startvalues
6200 , const Double_t *parlimitslo, const Double_t *parlimitshi
6201 , Double_t *fitparams, Double_t *fiterrors
6202 , Double_t *chiSqr, Int_t *ndf) const
6205 // Function for the fit
6209 Char_t funname[100];
6211 TF1 *ffitold = (TF1 *) gROOT->GetListOfFunctions()->FindObject(funname);
6216 TF1 *ffit = new TF1(funname,LanGauFun,fitrange[0],fitrange[1],4);
6217 ffit->SetParameters(startvalues);
6218 ffit->SetParNames("Width","MP","Area","GSigma");
6220 for (i = 0; i < 4; i++) {
6221 ffit->SetParLimits(i,parlimitslo[i],parlimitshi[i]);
6224 his->Fit(funname,"RB0"); // Fit within specified range, use ParLimits, do not plot
6226 ffit->GetParameters(fitparams); // Obtain fit parameters
6227 for (i = 0; i < 4; i++) {
6228 fiterrors[i] = ffit->GetParError(i); // Obtain fit parameter errors
6230 chiSqr[0] = ffit->GetChisquare(); // Obtain chi^2
6231 ndf[0] = ffit->GetNDF(); // Obtain ndf
6233 return (ffit); // Return fit function
6237 //_____________________________________________________________________________
6238 Int_t AliTRDCalibraFit::LanGauPro(const Double_t *params, Double_t &maxx, Double_t &fwhm)
6241 // Function for the fit
6254 Int_t maxcalls = 10000;
6256 // Search for maximum
6257 p = params[1] - 0.1 * params[0];
6258 step = 0.05 * params[0];
6262 while ((l != lold) && (i < maxcalls)) {
6266 l = LanGauFun(&x,params);
6268 step = -step / 10.0;
6273 if (i == maxcalls) {
6279 // Search for right x location of fy
6280 p = maxx + params[0];
6286 while ( (l != lold) && (i < maxcalls) ) {
6291 l = TMath::Abs(LanGauFun(&x,params) - fy);
6305 // Search for left x location of fy
6307 p = maxx - 0.5 * params[0];
6313 while ((l != lold) && (i < maxcalls)) {
6317 l = TMath::Abs(LanGauFun(&x,params) - fy);
6319 step = -step / 10.0;
6324 if (i == maxcalls) {
6333 //_____________________________________________________________________________
6334 Double_t AliTRDCalibraFit::GausConstant(const Double_t *x, const Double_t *par)
6337 // Gaus with identical mean
6340 Double_t gauss = par[0] * TMath::Gaus(x[0],0.0,par[1])+par[2];