1 /**************************************************************************
2 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
4 * Author: The ALICE Off-line Project. *
5 * Contributors are mentioned in the code where appropriate. *
7 * Permission to use, copy, modify and distribute this software and its *
8 * documentation strictly for non-commercial purposes is hereby granted *
9 * without fee, provided that the above copyright notice appears in all *
10 * copies and that both the copyright notice and this permission notice *
11 * appear in the supporting documentation. The authors make no claims *
12 * about the suitability of this software for any purpose. It is *
13 * provided "as is" without express or implied warranty. *
14 **************************************************************************/
18 /////////////////////////////////////////////////////////////////////////////////
22 // This class is for the TRD calibration of the relative gain factor, the drift velocity,
23 // the time 0 and the pad response function. It fits the histos.
24 // The 2D histograms or vectors (first converted in 1D histos) will be fitted
25 // if they have enough entries, otherwise the (default) value of the choosen database
26 // will be put. For the relative gain calibration the resulted factors will be globally
27 // normalized to the gain factors of the choosen database. It unables to precise
28 // previous calibration procedure.
29 // The function SetDebug enables the user to see:
30 // _fDebug = 0: nothing, only the values are written in the tree if wanted
31 // _fDebug = 1: a comparaison of the coefficients found and the default values
32 // in the choosen database.
33 // fCoef , histogram of the coefs as function of the calibration group number
34 // fDelta , histogram of the relative difference of the coef with the default
35 // value in the database as function of the calibration group number
36 // fError , dirstribution of this relative difference
37 // _fDebug = 2: only the fit of the choosen calibration group fFitVoir (SetFitVoir)
38 // _fDebug = 3: The coefficients in the choosen detector fDet (SetDet) as function of the
39 // pad row and col number
40 // _fDebug = 4; The coeffcicients in the choosen detector fDet (SetDet) like in the 3 but with
41 // also the comparaison histograms of the 1 for this detector
45 // R. Bailhache (R.Bailhache@gsi.de)
47 //////////////////////////////////////////////////////////////////////////////////////
53 #include <TProfile2D.h>
56 #include <TGraphErrors.h>
58 #include <TObjArray.h>
64 #include <TStopwatch.h>
67 #include <TDirectory.h>
69 #include <TTreeStream.h>
70 #include <TLinearFitter.h>
75 #include "AliCDBManager.h"
76 #include "AliMathBase.h"
78 #include "AliTRDCalibraFit.h"
79 #include "AliTRDCalibraMode.h"
80 #include "AliTRDCalibraVector.h"
81 #include "AliTRDCalibraVdriftLinearFit.h"
82 #include "AliTRDcalibDB.h"
83 #include "AliTRDgeometry.h"
84 #include "AliTRDpadPlane.h"
85 #include "AliTRDgeometry.h"
86 #include "./Cal/AliTRDCalROC.h"
87 #include "./Cal/AliTRDCalPad.h"
88 #include "./Cal/AliTRDCalDet.h"
91 ClassImp(AliTRDCalibraFit)
93 AliTRDCalibraFit* AliTRDCalibraFit::fgInstance = 0;
94 Bool_t AliTRDCalibraFit::fgTerminated = kFALSE;
96 //_____________singleton implementation_________________________________________________
97 AliTRDCalibraFit *AliTRDCalibraFit::Instance()
100 // Singleton implementation
103 if (fgTerminated != kFALSE) {
107 if (fgInstance == 0) {
108 fgInstance = new AliTRDCalibraFit();
115 //______________________________________________________________________________________
116 void AliTRDCalibraFit::Terminate()
119 // Singleton implementation
120 // Deletes the instance of this class
123 fgTerminated = kTRUE;
125 if (fgInstance != 0) {
132 //______________________________________________________________________________________
133 AliTRDCalibraFit::AliTRDCalibraFit()
136 ,fNumberOfBinsExpected(0)
138 ,fBeginFitCharge(3.5)
140 ,fTakeTheMaxPH(kFALSE)
147 ,fNumberFitSuccess(0)
154 ,fCalibraMode(new AliTRDCalibraMode())
159 ,fScaleFitFactor(0.0)
167 ,fCurrentCoefDetector(0x0)
168 ,fCurrentCoefDetector2(0x0)
173 // Default constructor
176 fGeo = new AliTRDgeometry();
178 // Current variables initialised
179 for (Int_t k = 0; k < 2; k++) {
180 fCurrentCoef[k] = 0.0;
181 fCurrentCoef2[k] = 0.0;
183 for (Int_t i = 0; i < 3; i++) {
189 //______________________________________________________________________________________
190 AliTRDCalibraFit::AliTRDCalibraFit(const AliTRDCalibraFit &c)
193 ,fNumberOfBinsExpected(c.fNumberOfBinsExpected)
195 ,fBeginFitCharge(c.fBeginFitCharge)
196 ,fFitPHPeriode(c.fFitPHPeriode)
197 ,fTakeTheMaxPH(c.fTakeTheMaxPH)
198 ,fT0Shift(c.fT0Shift)
199 ,fRangeFitPRF(c.fRangeFitPRF)
201 ,fMinEntries(c.fMinEntries)
203 ,fNumberFit(c.fNumberFit)
204 ,fNumberFitSuccess(c.fNumberFitSuccess)
205 ,fNumberEnt(c.fNumberEnt)
206 ,fStatisticMean(c.fStatisticMean)
208 ,fDebugLevel(c.fDebugLevel)
209 ,fFitVoir(c.fFitVoir)
210 ,fMagneticField(c.fMagneticField)
212 ,fCurrentCoefE(c.fCurrentCoefE)
213 ,fCurrentCoefE2(c.fCurrentCoefE2)
216 ,fScaleFitFactor(c.fScaleFitFactor)
217 ,fEntriesCurrent(c.fEntriesCurrent)
218 ,fCountDet(c.fCountDet)
224 ,fCurrentCoefDetector(0x0)
225 ,fCurrentCoefDetector2(0x0)
233 if(c.fCalibraMode) fCalibraMode = new AliTRDCalibraMode(*c.fCalibraMode);
235 //Current variables initialised
236 for (Int_t k = 0; k < 2; k++) {
237 fCurrentCoef[k] = 0.0;
238 fCurrentCoef2[k] = 0.0;
240 for (Int_t i = 0; i < 3; i++) {
244 if(c.fCalDet) fCalDet = new AliTRDCalDet(*c.fCalDet);
245 if(c.fCalDet2) fCalDet2 = new AliTRDCalDet(*c.fCalDet2);
247 if(c.fCalROC) fCalROC = new AliTRDCalROC(*c.fCalROC);
248 if(c.fCalROC2) fCalROC = new AliTRDCalROC(*c.fCalROC2);
250 fVectorFit.SetName(c.fVectorFit.GetName());
251 for(Int_t k = 0; k < c.fVectorFit.GetEntriesFast(); k++){
252 AliTRDFitInfo *fitInfo = new AliTRDFitInfo();
253 Int_t detector = ((AliTRDFitInfo *)c.fVectorFit.UncheckedAt(k))->GetDetector();
255 if (GetChamber(detector) == 2) {
261 Float_t *coef = new Float_t[ntotal];
262 for (Int_t i = 0; i < ntotal; i++) {
263 coef[i] = ((AliTRDFitInfo *)c.fVectorFit.UncheckedAt(k))->GetCoef()[i];
265 fitInfo->SetCoef(coef);
266 fitInfo->SetDetector(detector);
267 fVectorFit.Add((TObject *) fitInfo);
269 fVectorFit.SetName(c.fVectorFit.GetName());
270 for(Int_t k = 0; k < c.fVectorFit2.GetEntriesFast(); k++){
271 AliTRDFitInfo *fitInfo = new AliTRDFitInfo();
272 Int_t detector = ((AliTRDFitInfo *)c.fVectorFit2.UncheckedAt(k))->GetDetector();
274 if (GetChamber(detector) == 2) {
280 Float_t *coef = new Float_t[ntotal];
281 for (Int_t i = 0; i < ntotal; i++) {
282 coef[i] = ((AliTRDFitInfo *)c.fVectorFit2.UncheckedAt(k))->GetCoef()[i];
284 fitInfo->SetCoef(coef);
285 fitInfo->SetDetector(detector);
286 fVectorFit2.Add((TObject *) fitInfo);
291 fGeo = new AliTRDgeometry();
294 //____________________________________________________________________________________
295 AliTRDCalibraFit::~AliTRDCalibraFit()
298 // AliTRDCalibraFit destructor
300 if ( fDebugStreamer ) delete fDebugStreamer;
301 if ( fCalDet ) delete fCalDet;
302 if ( fCalDet2 ) delete fCalDet2;
303 if ( fCalROC ) delete fCalROC;
304 if ( fCalROC2 ) delete fCalROC2;
306 fVectorFit2.Delete();
312 //_____________________________________________________________________________
313 void AliTRDCalibraFit::Destroy()
325 //____________Functions fit Online CH2d________________________________________
326 Bool_t AliTRDCalibraFit::AnalyseCH(TH2I *ch)
329 // Fit the 1D histos, projections of the 2D ch on the Xaxis, for each
330 // calibration group normalized the resulted coefficients (to 1 normally)
333 // Set the calibration mode
334 const char *name = ch->GetTitle();
335 SetModeCalibration(name,0);
337 // Number of Ybins (detectors or groups of pads)
338 Int_t nbins = ch->GetNbinsX();// charge
339 Int_t nybins = ch->GetNbinsY();// groups number
340 if (!InitFit(nybins,0)) {
346 fStatisticMean = 0.0;
348 fNumberFitSuccess = 0;
350 // Init fCountDet and fCount
351 InitfCountDetAndfCount(0);
352 // Beginning of the loop betwwen dect1 and dect2
353 for (Int_t idect = fDect1; idect < fDect2; idect++) {
354 // Determination of fNnZ, fNnRphi, fNfragZ and fNfragRphi...
355 UpdatefCountDetAndfCount(idect,0);
356 ReconstructFitRowMinRowMax(idect, 0);
358 TH1I *projch = (TH1I *) ch->ProjectionX("projch",idect+1,idect+1,(Option_t *)"e");
359 projch->SetDirectory(0);
360 // Number of entries for this calibration group
361 Double_t nentries = 0.0;
363 for (Int_t k = 0; k < nbins; k++) {
364 Int_t binnb = (nbins+2)*(idect+1)+(k+1);
365 nentries += ch->GetBinContent(binnb);
366 mean += projch->GetBinCenter(k+1)*projch->GetBinContent(k+1);
367 projch->SetBinError(k+1,TMath::Sqrt(projch->GetBinContent(k+1)));
369 projch->SetEntries(nentries);
370 //printf("The number of entries for the group %d is %f\n",idect,nentries);
375 // Rebin and statistic stuff
377 projch = ReBin((TH1I *) projch);
379 // This detector has not enough statistics or was off
380 if (nentries <= fMinEntries) {
381 NotEnoughStatisticCH(idect);
382 if (fDebugLevel != 1) {
387 // Statistics of the group fitted
388 fStatisticMean += nentries;
393 case 0: FitMeanW((TH1 *) projch, nentries); break;
394 case 1: FitMean((TH1 *) projch, nentries, mean); break;
395 case 2: FitCH((TH1 *) projch, mean); break;
396 case 3: FitBisCH((TH1 *) projch, mean); break;
397 default: return kFALSE;
400 FillInfosFitCH(idect);
402 if (fDebugLevel != 1) {
407 if (fDebugLevel != 1) {
411 if (fNumberFit > 0) {
412 AliInfo(Form("There are %d with at least one entries. %d fits have been proceeded (sucessfully or not...). There is a mean statistic of: %d over these fitted histograms and %d successfulled fits",fNumberEnt, fNumberFit, (Int_t) fStatisticMean/fNumberFit, fNumberFitSuccess));
413 fStatisticMean = fStatisticMean / fNumberFit;
416 AliInfo(Form("There are %d with at least one entries. There is no fit!",fNumberEnt));
418 delete fDebugStreamer;
419 fDebugStreamer = 0x0;
423 //____________Functions fit Online CH2d________________________________________
424 Bool_t AliTRDCalibraFit::AnalyseCH(AliTRDCalibraVector *calvect)
427 // Reconstruct a 1D histo from the vectorCH for each calibration group,
428 // fit the histo, normalized the resulted coefficients (to 1 normally)
431 // Set the calibraMode
432 const char *name = calvect->GetNameCH();
433 SetModeCalibration(name,0);
435 // Number of Xbins (detectors or groups of pads)
436 if (!InitFit((432*calvect->GetDetCha0(0)+108*calvect->GetDetCha2(0)),0)) {
442 fStatisticMean = 0.0;
444 fNumberFitSuccess = 0;
446 // Init fCountDet and fCount
447 InitfCountDetAndfCount(0);
448 // Beginning of the loop between dect1 and dect2
449 for (Int_t idect = fDect1; idect < fDect2; idect++) {
450 // Determination of fNnZ, fNnRphi, fNfragZ and fNfragRphi...........
451 UpdatefCountDetAndfCount(idect,0);
452 ReconstructFitRowMinRowMax(idect,0);
454 Double_t nentries = 0.0;
457 Bool_t something = kTRUE;
458 if(!calvect->GetCHEntries(fCountDet)) something = kFALSE;
462 projch = calvect->ConvertVectorCHHisto(fCountDet,(idect-(fCount-(fCalibraMode->GetNfragZ(0)*fCalibraMode->GetNfragRphi(0)))),(const char *) name);
463 projch->SetDirectory(0);
464 for (Int_t k = 0; k < calvect->GetNumberBinCharge(); k++) {
465 nentries += projch->GetBinContent(k+1);
466 mean += projch->GetBinCenter(k+1)*projch->GetBinContent(k+1);
472 //printf("The number of entries for the group %d is %f\n",idect,nentries);
475 projch = ReBin((TH1F *) projch);
478 // This detector has not enough statistics or was not found in VectorCH
479 if (nentries <= fMinEntries) {
480 NotEnoughStatisticCH(idect);
481 if (fDebugLevel != 1) {
482 if(projch) delete projch;
486 // Statistic of the histos fitted
487 fStatisticMean += nentries;
492 case 0: FitMeanW((TH1 *) projch, nentries); break;
493 case 1: FitMean((TH1 *) projch, nentries, mean); break;
494 case 2: FitCH((TH1 *) projch, mean); break;
495 case 3: FitBisCH((TH1 *) projch, mean); break;
496 default: return kFALSE;
499 FillInfosFitCH(idect);
501 if (fDebugLevel != 1) {
506 if (fDebugLevel != 1) {
510 if (fNumberFit > 0) {
511 AliInfo(Form("There are %d with at least one entries. %d fits have been proceeded (sucessfully or not...). There is a mean statistic of: %d over these fitted histograms and %d successfulled fits",fNumberEnt, fNumberFit, (Int_t) fStatisticMean/fNumberFit, fNumberFitSuccess));
512 fStatisticMean = fStatisticMean / fNumberFit;
515 AliInfo(Form("There are %d with at least one entries. There is no fit!",fNumberEnt));
517 delete fDebugStreamer;
518 fDebugStreamer = 0x0;
521 //________________functions fit Online PH2d____________________________________
522 Bool_t AliTRDCalibraFit::AnalysePH(TProfile2D *ph)
525 // Take the 1D profiles (average pulse height), projections of the 2D PH
526 // on the Xaxis, for each calibration group
527 // Reconstruct a drift velocity
528 // A first calibration of T0 is also made using the same method
531 // Set the calibration mode
532 const char *name = ph->GetTitle();
533 SetModeCalibration(name,1);
535 // Number of Xbins (detectors or groups of pads)
536 Int_t nbins = ph->GetNbinsX();// time
537 Int_t nybins = ph->GetNbinsY();// calibration group
538 if (!InitFit(nybins,1)) {
544 fStatisticMean = 0.0;
546 fNumberFitSuccess = 0;
548 // Init fCountDet and fCount
549 InitfCountDetAndfCount(1);
550 // Beginning of the loop
551 for (Int_t idect = fDect1; idect < fDect2; idect++) {
552 // Determination of fNnZ, fNnRphi, fNfragZ and fNfragRphi.......
553 UpdatefCountDetAndfCount(idect,1);
554 ReconstructFitRowMinRowMax(idect,1);
556 TH1D *projph = (TH1D *) ph->ProjectionX("projph",idect+1,idect+1,(Option_t *) "e");
557 projph->SetDirectory(0);
558 // Number of entries for this calibration group
559 Double_t nentries = 0;
560 for (Int_t k = 0; k < nbins; k++) {
561 Int_t binnb = (nbins+2)*(idect+1)+(k+1);
562 nentries += ph->GetBinEntries(binnb);
567 //printf("The number of entries for the group %d is %f\n",idect,nentries);
568 // This detector has not enough statistics or was off
569 if (nentries <= fMinEntries) {
570 //printf("Not enough statistic!\n");
571 NotEnoughStatisticPH(idect);
572 if (fDebugLevel != 1) {
577 // Statistics of the histos fitted
579 fStatisticMean += nentries;
580 // Calcul of "real" coef
581 CalculVdriftCoefMean();
586 case 0: FitLagrangePoly((TH1 *) projph); break;
587 case 1: FitPente((TH1 *) projph); break;
588 case 2: FitPH((TH1 *) projph,(Int_t) (idect - fDect1)); break;
589 default: return kFALSE;
591 // Fill the tree if end of a detector or only the pointer to the branch!!!
592 FillInfosFitPH(idect);
594 if (fDebugLevel != 1) {
599 if (fNumberFit > 0) {
600 AliInfo(Form("There are %d with at least one entries. %d fits have been proceeded (sucessfully or not...). There is a mean statistic of: %d over these fitted histograms and %d successfulled fits",fNumberEnt, fNumberFit, (Int_t) fStatisticMean/fNumberFit,fNumberFitSuccess));
601 fStatisticMean = fStatisticMean / fNumberFit;
604 AliInfo(Form("There are %d with at least one entries. There is no fit!",fNumberEnt));
606 delete fDebugStreamer;
607 fDebugStreamer = 0x0;
610 //____________Functions fit Online PH2d________________________________________
611 Bool_t AliTRDCalibraFit::AnalysePH(AliTRDCalibraVector *calvect)
614 // Reconstruct the average pulse height from the vectorPH for each
616 // Reconstruct a drift velocity
617 // A first calibration of T0 is also made using the same method (slope method)
620 // Set the calibration mode
621 const char *name = calvect->GetNamePH();
622 SetModeCalibration(name,1);
624 // Number of Xbins (detectors or groups of pads)
625 if (!InitFit((432*calvect->GetDetCha0(1)+108*calvect->GetDetCha2(1)),1)) {
631 fStatisticMean = 0.0;
633 fNumberFitSuccess = 0;
635 // Init fCountDet and fCount
636 InitfCountDetAndfCount(1);
637 // Beginning of the loop
638 for (Int_t idect = fDect1; idect < fDect2; idect++) {
639 // Determination of fNnZ, fNnRphi, fNfragZ and fNfragRphi...........
640 UpdatefCountDetAndfCount(idect,1);
641 ReconstructFitRowMinRowMax(idect,1);
645 Bool_t something = kTRUE;
646 if(!calvect->GetPHEntries(fCountDet)) something = kFALSE;
650 projph = CorrectTheError((TGraphErrors *) (calvect->ConvertVectorPHTGraphErrors(fCountDet,(idect-(fCount-(fCalibraMode->GetNfragZ(1)*fCalibraMode->GetNfragRphi(1)))),(const char *) name)));
651 projph->SetDirectory(0);
653 //printf("The number of entries for the group %d is %d\n",idect,fEntriesCurrent);
654 // This detector has not enough statistics or was off
655 if (fEntriesCurrent <= fMinEntries) {
656 //printf("Not enough stat!\n");
657 NotEnoughStatisticPH(idect);
658 if (fDebugLevel != 1) {
659 if(projph) delete projph;
663 // Statistic of the histos fitted
665 fStatisticMean += fEntriesCurrent;
666 // Calcul of "real" coef
667 CalculVdriftCoefMean();
672 case 0: FitLagrangePoly((TH1 *) projph); break;
673 case 1: FitPente((TH1 *) projph); break;
674 case 2: FitPH((TH1 *) projph,(Int_t) (idect - fDect1)); break;
675 default: return kFALSE;
677 // Fill the tree if end of a detector or only the pointer to the branch!!!
678 FillInfosFitPH(idect);
680 if (fDebugLevel != 1) {
686 if (fNumberFit > 0) {
687 AliInfo(Form("There are %d with at least one entries. %d fits have been proceeded (sucessfully or not...). There is a mean statistic of: %d over these fitted histograms and %d successfulled fits",fNumberEnt, fNumberFit, (Int_t) fStatisticMean/fNumberFit,fNumberFitSuccess));
688 fStatisticMean = fStatisticMean / fNumberFit;
691 AliInfo(Form("There are %d with at least one entries. There is no fit!",fNumberEnt));
693 delete fDebugStreamer;
694 fDebugStreamer = 0x0;
697 //____________Functions fit Online PRF2d_______________________________________
698 Bool_t AliTRDCalibraFit::AnalysePRF(TProfile2D *prf)
701 // Take the 1D profiles (pad response function), projections of the 2D PRF
702 // on the Xaxis, for each calibration group
703 // Fit with a gaussian to reconstruct the sigma of the pad response function
706 // Set the calibration mode
707 const char *name = prf->GetTitle();
708 SetModeCalibration(name,2);
710 // Number of Ybins (detectors or groups of pads)
711 Int_t nybins = prf->GetNbinsY();// calibration groups
712 Int_t nbins = prf->GetNbinsX();// bins
713 Int_t nbg = GetNumberOfGroupsPRF((const char *)prf->GetTitle());
714 if((nbg > 0) || (nbg == -1)) return kFALSE;
715 if (!InitFit(nybins,2)) {
721 fStatisticMean = 0.0;
723 fNumberFitSuccess = 0;
725 // Init fCountDet and fCount
726 InitfCountDetAndfCount(2);
727 // Beginning of the loop
728 for (Int_t idect = fDect1; idect < fDect2; idect++) {
729 // Determination of fNnZ, fNnRphi, fNfragZ and fNfragRphi......
730 UpdatefCountDetAndfCount(idect,2);
731 ReconstructFitRowMinRowMax(idect,2);
733 TH1D *projprf = (TH1D *) prf->ProjectionX("projprf",idect+1,idect+1,(Option_t *) "e");
734 projprf->SetDirectory(0);
735 // Number of entries for this calibration group
736 Double_t nentries = 0;
737 for (Int_t k = 0; k < nbins; k++) {
738 Int_t binnb = (nbins+2)*(idect+1)+(k+1);
739 nentries += prf->GetBinEntries(binnb);
741 if(nentries > 0) fNumberEnt++;
742 // This detector has not enough statistics or was off
743 if (nentries <= fMinEntries) {
744 NotEnoughStatisticPRF(idect);
745 if (fDebugLevel != 1) {
750 // Statistics of the histos fitted
752 fStatisticMean += nentries;
753 // Calcul of "real" coef
758 case 0: FitPRF((TH1 *) projprf); break;
759 case 1: RmsPRF((TH1 *) projprf); break;
760 default: return kFALSE;
762 // Fill the tree if end of a detector or only the pointer to the branch!!!
763 FillInfosFitPRF(idect);
765 if (fDebugLevel != 1) {
770 if (fNumberFit > 0) {
771 AliInfo(Form("There are %d with at least one entries.",fNumberEnt));
772 AliInfo(Form("%d fits have been proceeded (sucessfully or not...).",fNumberFit));
773 AliInfo(Form("There is a mean statistic of: %d over these fitted histograms and %d successfulled fits"
774 ,(Int_t) fStatisticMean/fNumberFit,fNumberFitSuccess));
775 fStatisticMean = fStatisticMean / fNumberFit;
778 AliInfo(Form("There are %d with at least one entries. There is no fit!",fNumberEnt));
780 delete fDebugStreamer;
781 fDebugStreamer = 0x0;
784 //____________Functions fit Online PRF2d_______________________________________
785 Bool_t AliTRDCalibraFit::AnalysePRFMarianFit(TProfile2D *prf)
788 // Take the 1D profiles (pad response function), projections of the 2D PRF
789 // on the Xaxis, for each calibration group
790 // Fit with a gaussian to reconstruct the sigma of the pad response function
793 // Set the calibration mode
794 const char *name = prf->GetTitle();
795 SetModeCalibration(name,2);
797 // Number of Ybins (detectors or groups of pads)
798 TAxis *xprf = prf->GetXaxis();
799 TAxis *yprf = prf->GetYaxis();
800 Int_t nybins = yprf->GetNbins();// calibration groups
801 Int_t nbins = xprf->GetNbins();// bins
802 Float_t lowedge = (Float_t) xprf->GetBinLowEdge(1);//lowedge in bins
803 Float_t upedge = (Float_t) xprf->GetBinUpEdge(nbins);//upedge in bins
804 Int_t nbg = GetNumberOfGroupsPRF((const char *)name);
805 if(nbg == -1) return kFALSE;
806 if(nbg > 0) fMethod = 1;
808 if (!InitFit(nybins,2)) {
814 fStatisticMean = 0.0;
816 fNumberFitSuccess = 0;
818 // Init fCountDet and fCount
819 InitfCountDetAndfCount(2);
820 // Beginning of the loop
821 for (Int_t idect = fDect1; idect < fDect2; idect++) {
822 // Determination of fNnZ, fNnRphi, fNfragZ and fNfragRphi.......
823 UpdatefCountDetAndfCount(idect,2);
824 ReconstructFitRowMinRowMax(idect,2);
825 // Build the array of entries and sum
826 TArrayD arraye = TArrayD(nbins);
827 TArrayD arraym = TArrayD(nbins);
828 TArrayD arrayme = TArrayD(nbins);
829 Double_t nentries = 0;
830 //printf("nbins %d\n",nbins);
831 for (Int_t k = 0; k < nbins; k++) {
832 Int_t binnb = (nbins+2)*(idect+1)+(k+1);
833 Double_t entries = (Double_t)prf->GetBinEntries(binnb);
834 Double_t mean = (Double_t)prf->GetBinContent(binnb);
835 Double_t error = (Double_t)prf->GetBinError(binnb);
836 //printf("for %d we have %f\n",k,entries);
838 arraye.AddAt(entries,k);
839 arraym.AddAt(mean,k);
840 arrayme.AddAt(error,k);
842 if(nentries > 0) fNumberEnt++;
843 //printf("The number of entries for the group %d is %f\n",idect,nentries);
844 // This detector has not enough statistics or was off
845 if (nentries <= fMinEntries) {
846 NotEnoughStatisticPRF(idect);
849 // Statistics of the histos fitted
851 fStatisticMean += nentries;
852 // Calcul of "real" coef
857 case 0: FitPRFGausMI( arraye.GetArray(), arraym.GetArray(), arrayme.GetArray(), nbins, lowedge, upedge); break;
858 case 1: FitTnpRange( arraye.GetArray(), arraym.GetArray(), arrayme.GetArray(), nbg, nbins); break;
859 default: return kFALSE;
861 // Fill the tree if end of a detector or only the pointer to the branch!!!
862 FillInfosFitPRF(idect);
865 if (fNumberFit > 0) {
866 AliInfo(Form("There are %d with at least one entries.",fNumberEnt));
867 AliInfo(Form("%d fits have been proceeded (sucessfully or not...).",fNumberFit));
868 AliInfo(Form("There is a mean statistic of: %d over these fitted histograms and %d successfulled fits"
869 ,(Int_t) fStatisticMean/fNumberFit,fNumberFitSuccess));
870 fStatisticMean = fStatisticMean / fNumberFit;
873 AliInfo(Form("There are %d with at least one entries. There is no fit!",fNumberEnt));
875 delete fDebugStreamer;
876 fDebugStreamer = 0x0;
879 //____________Functions fit Online PRF2d_______________________________________
880 Bool_t AliTRDCalibraFit::AnalysePRF(AliTRDCalibraVector *calvect)
883 // Reconstruct the 1D histo (pad response function) from the vectorPRD for
884 // each calibration group
885 // Fit with a gaussian to reconstruct the sigma of the pad response function
888 // Set the calibra mode
889 const char *name = calvect->GetNamePRF();
890 SetModeCalibration(name,2);
891 //printf("test0 %s\n",name);
893 // Number of Xbins (detectors or groups of pads)
894 if (!InitFit((432*calvect->GetDetCha0(2)+108*calvect->GetDetCha2(2)),2)) {
899 ///printf("test2\n");
902 fStatisticMean = 0.0;
904 fNumberFitSuccess = 0;
906 // Init fCountDet and fCount
907 InitfCountDetAndfCount(2);
908 // Beginning of the loop
909 for (Int_t idect = fDect1; idect < fDect2; idect++) {
910 // Determination of fNnZ, fNnRphi, fNfragZ and fNfragRphi........
911 UpdatefCountDetAndfCount(idect,2);
912 ReconstructFitRowMinRowMax(idect,2);
916 Bool_t something = kTRUE;
917 if(!calvect->GetPRFEntries(fCountDet)) something = kFALSE;
921 projprf = CorrectTheError((TGraphErrors *) (calvect->ConvertVectorPRFTGraphErrors(fCountDet,(idect-(fCount-(fCalibraMode->GetNfragZ(1)*fCalibraMode->GetNfragRphi(1)))),(const char *) name)));
922 projprf->SetDirectory(0);
924 // This detector has not enough statistics or was off
925 if (fEntriesCurrent <= fMinEntries) {
926 NotEnoughStatisticPRF(idect);
927 if (fDebugLevel != 1) {
928 if(projprf) delete projprf;
932 // Statistic of the histos fitted
934 fStatisticMean += fEntriesCurrent;
935 // Calcul of "real" coef
940 case 1: FitPRF((TH1 *) projprf); break;
941 case 2: RmsPRF((TH1 *) projprf); break;
942 default: return kFALSE;
944 // Fill the tree if end of a detector or only the pointer to the branch!!!
945 FillInfosFitPRF(idect);
947 if (fDebugLevel != 1) {
952 if (fNumberFit > 0) {
953 AliInfo(Form("There are %d with at least one entries. %d fits have been proceeded (sucessfully or not...). There is a mean statistic of: %d over these fitted histograms and %d successfulled fits",fNumberEnt, fNumberFit, (Int_t) fStatisticMean/fNumberFit,fNumberFitSuccess));
956 AliInfo(Form("There are %d with at least one entries. There is no fit!",fNumberEnt));
958 delete fDebugStreamer;
959 fDebugStreamer = 0x0;
962 //____________Functions fit Online PRF2d_______________________________________
963 Bool_t AliTRDCalibraFit::AnalysePRFMarianFit(AliTRDCalibraVector *calvect)
966 // Reconstruct the 1D histo (pad response function) from the vectorPRD for
967 // each calibration group
968 // Fit with a gaussian to reconstruct the sigma of the pad response function
971 // Set the calibra mode
972 const char *name = calvect->GetNamePRF();
973 SetModeCalibration(name,2);
974 //printf("test0 %s\n",name);
975 Int_t nbg = GetNumberOfGroupsPRF((const char *)name);
976 printf("test1 %d\n",nbg);
977 if(nbg == -1) return kFALSE;
978 if(nbg > 0) fMethod = 1;
980 // Number of Xbins (detectors or groups of pads)
981 if (!InitFit((432*calvect->GetDetCha0(2)+108*calvect->GetDetCha2(2)),2)) {
989 fStatisticMean = 0.0;
991 fNumberFitSuccess = 0;
995 Double_t *arrayx = 0;
996 Double_t *arraye = 0;
997 Double_t *arraym = 0;
998 Double_t *arrayme = 0;
999 Float_t lowedge = 0.0;
1000 Float_t upedge = 0.0;
1001 // Init fCountDet and fCount
1002 InitfCountDetAndfCount(2);
1003 // Beginning of the loop
1004 for (Int_t idect = fDect1; idect < fDect2; idect++) {
1005 // Determination of fNnZ, fNnRphi, fNfragZ and fNfragRphi......
1006 UpdatefCountDetAndfCount(idect,2);
1007 ReconstructFitRowMinRowMax(idect,2);
1009 TGraphErrors *projprftree = 0x0;
1010 fEntriesCurrent = 0;
1011 Bool_t something = kTRUE;
1012 if(!calvect->GetPRFEntries(fCountDet)) something = kFALSE;
1014 TString name("PRF");
1016 projprftree = calvect->ConvertVectorPRFTGraphErrors(fCountDet,(idect-(fCount-(fCalibraMode->GetNfragZ(1)*fCalibraMode->GetNfragRphi(1)))),(const char *) name);
1017 nbins = projprftree->GetN();
1018 arrayx = (Double_t *)projprftree->GetX();
1019 arraye = (Double_t *)projprftree->GetEX();
1020 arraym = (Double_t *)projprftree->GetY();
1021 arrayme = (Double_t *)projprftree->GetEY();
1022 Float_t step = arrayx[1]-arrayx[0];
1023 lowedge = arrayx[0] - step/2.0;
1024 upedge = arrayx[(nbins-1)] + step/2.0;
1025 //printf("nbins est %d\n",nbins);
1026 for(Int_t k = 0; k < nbins; k++){
1027 fEntriesCurrent += (Int_t)arraye[k];
1028 //printf("for %d we have %f, %f\n",k,arraye[k],((projprftree->GetEX())[k]));
1029 if(arraye[k]>0.0) arrayme[k] = TMath::Sqrt(TMath::Abs(arrayme[k]-arraym[k]*arraym[k])/arraye[k]);
1031 if(fEntriesCurrent > 0) fNumberEnt++;
1033 //printf("The number of entries for the group %d is %d\n",idect,fEntriesCurrent);
1034 // This detector has not enough statistics or was off
1035 if (fEntriesCurrent <= fMinEntries) {
1036 NotEnoughStatisticPRF(idect);
1037 if(projprftree) delete projprftree;
1040 // Statistic of the histos fitted
1042 fStatisticMean += fEntriesCurrent;
1043 // Calcul of "real" coef
1044 CalculPRFCoefMean();
1048 case 0: FitPRFGausMI(arraye,arraym,arrayme,nbins,lowedge,upedge); break;
1049 case 1: FitTnpRange(arraye,arraym,arrayme,nbg,nbins); break;
1050 default: return kFALSE;
1052 // Fill the tree if end of a detector or only the pointer to the branch!!!
1053 FillInfosFitPRF(idect);
1055 if (fDebugLevel != 1) {
1060 if (fNumberFit > 0) {
1061 AliInfo(Form("There are %d with at least one entries. %d fits have been proceeded (sucessfully or not...). There is a mean statistic of: %d over these fitted histograms and %d successfulled fits",fNumberEnt, fNumberFit, (Int_t) fStatisticMean/fNumberFit,fNumberFitSuccess));
1064 AliInfo(Form("There are %d with at least one entries. There is no fit!",fNumberEnt));
1066 delete fDebugStreamer;
1067 fDebugStreamer = 0x0;
1070 //____________Functions fit Online CH2d________________________________________
1071 Bool_t AliTRDCalibraFit::AnalyseLinearFitters(AliTRDCalibraVdriftLinearFit *calivdli)
1074 // The linear method
1077 fStatisticMean = 0.0;
1079 fNumberFitSuccess = 0;
1081 if(!InitFitLinearFitter()) return kFALSE;
1084 for(Int_t idet = 0; idet < 540; idet++){
1087 //printf("detector number %d\n",idet);
1092 fEntriesCurrent = 0;
1094 Bool_t here = calivdli->GetParam(idet,¶m);
1095 Bool_t heree = calivdli->GetError(idet,&error);
1096 //printf("here %d and heree %d\n",here, heree);
1098 fEntriesCurrent = (Int_t) error[2];
1101 //printf("Number of entries %d\n",fEntriesCurrent);
1102 // Nothing found or not enough statistic
1103 if((!heree) || (!here) || (fEntriesCurrent <= fMinEntries)) {
1104 NotEnoughStatisticLinearFitter();
1111 fStatisticMean += fEntriesCurrent;
1114 if((-(param[1])) <= 0.0) {
1115 NotEnoughStatisticLinearFitter();
1119 // CalculDatabaseVdriftandTan
1120 CalculVdriftLorentzCoef();
1123 fNumberFitSuccess ++;
1125 // Put the fCurrentCoef
1126 fCurrentCoef[0] = -param[1];
1127 // here the database must be the one of the reconstruction for the lorentz angle....
1128 fCurrentCoef2[0] = (param[0]+fCurrentCoef[1]*fCurrentCoef2[1])/fCurrentCoef[0];
1129 fCurrentCoefE = error[1];
1130 fCurrentCoefE2 = error[0];
1131 if((fCurrentCoef2[0] != 0.0) && (param[0] != 0.0)){
1132 fCurrentCoefE2 = (fCurrentCoefE2/param[0]+fCurrentCoefE/fCurrentCoef[0])*fCurrentCoef2[0];
1136 FillInfosFitLinearFitter();
1141 if (fNumberFit > 0) {
1142 AliInfo(Form("There are %d with at least one entries. %d fits have been proceeded (sucessfully or not...). There is a mean statistic of: %d over these fitted histograms and %d successfulled fits",fNumberEnt, fNumberFit, (Int_t) fStatisticMean/fNumberFit,fNumberFitSuccess));
1145 AliInfo(Form("There are %d with at least one entries. There is no fit!",fNumberEnt));
1147 delete fDebugStreamer;
1148 fDebugStreamer = 0x0;
1152 //____________Functions for seeing if the pad is really okey___________________
1153 //_____________________________________________________________________________
1154 Int_t AliTRDCalibraFit::GetNumberOfGroupsPRF(const char* nametitle)
1157 // Get numberofgroupsprf
1161 const Char_t *pattern0 = "Ngp0";
1162 const Char_t *pattern1 = "Ngp1";
1163 const Char_t *pattern2 = "Ngp2";
1164 const Char_t *pattern3 = "Ngp3";
1165 const Char_t *pattern4 = "Ngp4";
1166 const Char_t *pattern5 = "Ngp5";
1167 const Char_t *pattern6 = "Ngp6";
1170 if (strstr(nametitle,pattern0)) {
1173 if (strstr(nametitle,pattern1)) {
1176 if (strstr(nametitle,pattern2)) {
1179 if (strstr(nametitle,pattern3)) {
1182 if (strstr(nametitle,pattern4)) {
1185 if (strstr(nametitle,pattern5)) {
1188 if (strstr(nametitle,pattern6)){
1195 //_____________________________________________________________________________
1196 Bool_t AliTRDCalibraFit::SetModeCalibration(const char *name, Int_t i)
1199 // Set fNz[i] and fNrphi[i] of the AliTRDCalibraFit::Instance()
1200 // corresponding to the given name
1203 if(!SetNzFromTObject(name,i)) return kFALSE;
1204 if(!SetNrphiFromTObject(name,i)) return kFALSE;
1209 //_____________________________________________________________________________
1210 Bool_t AliTRDCalibraFit::SetNrphiFromTObject(const char *name, Int_t i)
1213 // Set fNrphi[i] of the AliTRDCalibraFit::Instance()
1214 // corresponding to the given TObject
1218 const Char_t *patternrphi0 = "Nrphi0";
1219 const Char_t *patternrphi1 = "Nrphi1";
1220 const Char_t *patternrphi2 = "Nrphi2";
1221 const Char_t *patternrphi3 = "Nrphi3";
1222 const Char_t *patternrphi4 = "Nrphi4";
1223 const Char_t *patternrphi5 = "Nrphi5";
1224 const Char_t *patternrphi6 = "Nrphi6";
1227 if (strstr(name,patternrphi0)) {
1228 fCalibraMode->SetNrphi(i ,0);
1231 if (strstr(name,patternrphi1)) {
1232 fCalibraMode->SetNrphi(i, 1);
1235 if (strstr(name,patternrphi2)) {
1236 fCalibraMode->SetNrphi(i, 2);
1239 if (strstr(name,patternrphi3)) {
1240 fCalibraMode->SetNrphi(i, 3);
1243 if (strstr(name,patternrphi4)) {
1244 fCalibraMode->SetNrphi(i, 4);
1247 if (strstr(name,patternrphi5)) {
1248 fCalibraMode->SetNrphi(i, 5);
1251 if (strstr(name,patternrphi6)) {
1252 fCalibraMode->SetNrphi(i, 6);
1256 fCalibraMode->SetNrphi(i ,0);
1260 //_____________________________________________________________________________
1261 Bool_t AliTRDCalibraFit::SetNzFromTObject(const char *name, Int_t i)
1264 // Set fNz[i] of the AliTRDCalibraFit::Instance()
1265 // corresponding to the given TObject
1269 const Char_t *patternz0 = "Nz0";
1270 const Char_t *patternz1 = "Nz1";
1271 const Char_t *patternz2 = "Nz2";
1272 const Char_t *patternz3 = "Nz3";
1273 const Char_t *patternz4 = "Nz4";
1275 if (strstr(name,patternz0)) {
1276 fCalibraMode->SetNz(i, 0);
1279 if (strstr(name,patternz1)) {
1280 fCalibraMode->SetNz(i ,1);
1283 if (strstr(name,patternz2)) {
1284 fCalibraMode->SetNz(i ,2);
1287 if (strstr(name,patternz3)) {
1288 fCalibraMode->SetNz(i ,3);
1291 if (strstr(name,patternz4)) {
1292 fCalibraMode->SetNz(i ,4);
1296 fCalibraMode->SetNz(i ,0);
1299 //_____________________________________________________________________________
1300 AliTRDCalDet *AliTRDCalibraFit::CreateDetObjectVdrift(TObjArray *vectorFit, Bool_t perdetector)
1303 // It creates the AliTRDCalDet object from the AliTRDFitInfo
1304 // It takes the mean value of the coefficients per detector
1305 // This object has to be written in the database
1308 // Create the DetObject
1309 AliTRDCalDet *object = new AliTRDCalDet("ChamberVdrift","TRD drift velocities (detector value)");
1311 Int_t loop = (Int_t) vectorFit->GetEntriesFast();
1312 if(loop != 540) AliInfo("The Vector Fit is not complete!");
1313 Int_t detector = -1;
1314 Float_t value = 0.0;
1316 for (Int_t k = 0; k < loop; k++) {
1317 detector = ((AliTRDFitInfo *) vectorFit->At(k))->GetDetector();
1320 mean = TMath::Abs(((AliTRDFitInfo *) vectorFit->At(k))->GetCoef()[0]);
1324 Int_t rowMax = fGeo->GetRowMax(GetPlane(detector),GetChamber(detector),GetSector(detector));
1325 Int_t colMax = fGeo->GetColMax(GetPlane(detector));
1326 for (Int_t row = 0; row < rowMax; row++) {
1327 for (Int_t col = 0; col < colMax; col++) {
1328 value = ((AliTRDFitInfo *) vectorFit->At(k))->GetCoef()[(Int_t)(col*rowMax+row)];
1329 mean += TMath::Abs(value);
1333 if(count > 0) mean = mean/count;
1335 object->SetValue(detector,mean);
1340 //_____________________________________________________________________________
1341 AliTRDCalDet *AliTRDCalibraFit::CreateDetObjectGain(TObjArray *vectorFit, Double_t scaleFitFactor, Bool_t perdetector)
1344 // It creates the AliTRDCalDet object from the AliTRDFitInfo
1345 // It takes the mean value of the coefficients per detector
1346 // This object has to be written in the database
1349 // Create the DetObject
1350 AliTRDCalDet *object = new AliTRDCalDet("ChamberGainFactor","GainFactor (detector value)");
1353 Int_t loop = (Int_t) vectorFit->GetEntriesFast();
1354 if(loop != 540) AliInfo("The Vector Fit is not complete!");
1355 Int_t detector = -1;
1356 Float_t value = 0.0;
1358 for (Int_t k = 0; k < loop; k++) {
1359 detector = ((AliTRDFitInfo *) vectorFit->At(k))->GetDetector();
1362 value = ((AliTRDFitInfo *) vectorFit->At(k))->GetCoef()[0];
1363 if(value > 0) value = value*scaleFitFactor;
1364 mean = TMath::Abs(value);
1368 Int_t rowMax = fGeo->GetRowMax(GetPlane(detector),GetChamber(detector),GetSector(detector));
1369 Int_t colMax = fGeo->GetColMax(GetPlane(detector));
1370 for (Int_t row = 0; row < rowMax; row++) {
1371 for (Int_t col = 0; col < colMax; col++) {
1372 value = ((AliTRDFitInfo *) vectorFit->At(k))->GetCoef()[(Int_t)(col*rowMax+row)];
1373 if(value > 0) value = value*scaleFitFactor;
1374 mean += TMath::Abs(value);
1378 if(count > 0) mean = mean/count;
1380 object->SetValue(detector,mean);
1385 //_____________________________________________________________________________
1386 AliTRDCalDet *AliTRDCalibraFit::CreateDetObjectT0(TObjArray *vectorFit, Bool_t perdetector)
1389 // It creates the AliTRDCalDet object from the AliTRDFitInfo2
1390 // It takes the min value of the coefficients per detector
1391 // This object has to be written in the database
1394 // Create the DetObject
1395 AliTRDCalDet *object = new AliTRDCalDet("ChamberT0","T0 (detector value)");
1397 Int_t loop = (Int_t) vectorFit->GetEntriesFast();
1398 if(loop != 540) AliInfo("The Vector Fit is not complete!");
1399 Int_t detector = -1;
1400 Float_t value = 0.0;
1402 for (Int_t k = 0; k < loop; k++) {
1403 detector = ((AliTRDFitInfo *) vectorFit->At(k))->GetDetector();
1404 Float_t min = 100.0;
1406 min = ((AliTRDFitInfo *) vectorFit->At(k))->GetCoef()[0];
1409 Int_t rowMax = fGeo->GetRowMax(GetPlane(detector),GetChamber(detector),GetSector(detector));
1410 Int_t colMax = fGeo->GetColMax(GetPlane(detector));
1411 for (Int_t row = 0; row < rowMax; row++) {
1412 for (Int_t col = 0; col < colMax; col++) {
1413 value = ((AliTRDFitInfo *) vectorFit->At(k))->GetCoef()[(Int_t)(col*rowMax+row)];
1414 if(min > value) min = value;
1418 object->SetValue(detector,min);
1424 //_____________________________________________________________________________
1425 AliTRDCalDet *AliTRDCalibraFit::CreateDetObjectLorentzAngle(TObjArray *vectorFit)
1428 // It creates the AliTRDCalDet object from the AliTRDFitInfo2
1429 // It takes the min value of the coefficients per detector
1430 // This object has to be written in the database
1433 // Create the DetObject
1434 AliTRDCalDet *object = new AliTRDCalDet("tan(lorentzangle)","tan(lorentzangle) (detector value)");
1437 Int_t loop = (Int_t) vectorFit->GetEntriesFast();
1438 if(loop != 540) AliInfo("The Vector Fit is not complete!");
1439 Int_t detector = -1;
1440 Float_t value = 0.0;
1442 for (Int_t k = 0; k < loop; k++) {
1443 detector = ((AliTRDFitInfo *) vectorFit->At(k))->GetDetector();
1445 Int_t rowMax = fGeo->GetRowMax(GetPlane(detector),GetChamber(detector),GetSector(detector));
1446 Int_t colMax = fGeo->GetColMax(GetPlane(detector));
1447 Float_t min = 100.0;
1448 for (Int_t row = 0; row < rowMax; row++) {
1449 for (Int_t col = 0; col < colMax; col++) {
1450 value = ((AliTRDFitInfo *) fVectorFit2.At(k))->GetCoef()[(Int_t)(col*rowMax+row)];
1451 mean += -TMath::Abs(value);
1455 if(count > 0) mean = mean/count;
1457 value = ((AliTRDFitInfo *) vectorFit->At(k))->GetCoef()[0];
1458 object->SetValue(detector,-TMath::Abs(value));
1464 //_____________________________________________________________________________
1465 TObject *AliTRDCalibraFit::CreatePadObjectGain(TObjArray *vectorFit, Double_t scaleFitFactor, AliTRDCalDet *detobject)
1468 // It Creates the AliTRDCalPad object from AliTRDFitInfo
1469 // You need first to create the object for the detectors,
1470 // where the mean value is put.
1471 // This object has to be written in the database
1474 // Create the DetObject
1475 AliTRDCalPad *object = new AliTRDCalPad("GainFactor","GainFactor (local variations)");
1478 for(Int_t k = 0; k < 540; k++){
1479 AliTRDCalROC *calROC = object->GetCalROC(k);
1480 Int_t nchannels = calROC->GetNchannels();
1481 for(Int_t ch = 0; ch < nchannels; ch++){
1482 calROC->SetValue(ch,1.0);
1488 Int_t loop = (Int_t) vectorFit->GetEntriesFast();
1489 if(loop != 540) AliInfo("The Vector Fit is not complete!");
1490 Int_t detector = -1;
1491 Float_t value = 0.0;
1493 for (Int_t k = 0; k < loop; k++) {
1494 detector = ((AliTRDFitInfo *) vectorFit->At(k))->GetDetector();
1495 AliTRDCalROC *calROC = object->GetCalROC(detector);
1496 Float_t mean = detobject->GetValue(detector);
1497 if(mean == 0) continue;
1498 Int_t rowMax = calROC->GetNrows();
1499 Int_t colMax = calROC->GetNcols();
1500 for (Int_t row = 0; row < rowMax; row++) {
1501 for (Int_t col = 0; col < colMax; col++) {
1502 value = ((AliTRDFitInfo *) vectorFit->At(k))->GetCoef()[(Int_t)(col*rowMax+row)];
1503 if(value > 0) value = value*scaleFitFactor;
1504 calROC->SetValue(col,row,TMath::Abs(value)/mean);
1512 //_____________________________________________________________________________
1513 TObject *AliTRDCalibraFit::CreatePadObjectVdrift(TObjArray *vectorFit, AliTRDCalDet *detobject)
1516 // It Creates the AliTRDCalPad object from AliTRDFitInfo
1517 // You need first to create the object for the detectors,
1518 // where the mean value is put.
1519 // This object has to be written in the database
1522 // Create the DetObject
1523 AliTRDCalPad *object = new AliTRDCalPad("LocalVdrift","TRD drift velocities (local variations)");
1526 for(Int_t k = 0; k < 540; k++){
1527 AliTRDCalROC *calROC = object->GetCalROC(k);
1528 Int_t nchannels = calROC->GetNchannels();
1529 for(Int_t ch = 0; ch < nchannels; ch++){
1530 calROC->SetValue(ch,1.0);
1536 Int_t loop = (Int_t) vectorFit->GetEntriesFast();
1537 if(loop != 540) AliInfo("The Vector Fit is not complete!");
1538 Int_t detector = -1;
1539 Float_t value = 0.0;
1541 for (Int_t k = 0; k < loop; k++) {
1542 detector = ((AliTRDFitInfo *) vectorFit->At(k))->GetDetector();
1543 AliTRDCalROC *calROC = object->GetCalROC(detector);
1544 Float_t mean = detobject->GetValue(detector);
1545 if(mean == 0) continue;
1546 Int_t rowMax = calROC->GetNrows();
1547 Int_t colMax = calROC->GetNcols();
1548 for (Int_t row = 0; row < rowMax; row++) {
1549 for (Int_t col = 0; col < colMax; col++) {
1550 value = ((AliTRDFitInfo *) vectorFit->At(k))->GetCoef()[(Int_t)(col*rowMax+row)];
1551 calROC->SetValue(col,row,TMath::Abs(value)/mean);
1559 //_____________________________________________________________________________
1560 TObject *AliTRDCalibraFit::CreatePadObjectT0(TObjArray *vectorFit, AliTRDCalDet *detobject)
1563 // It Creates the AliTRDCalPad object from AliTRDFitInfo2
1564 // You need first to create the object for the detectors,
1565 // where the mean value is put.
1566 // This object has to be written in the database
1569 // Create the DetObject
1570 AliTRDCalPad *object = new AliTRDCalPad("LocalT0","T0 (local variations)");
1573 for(Int_t k = 0; k < 540; k++){
1574 AliTRDCalROC *calROC = object->GetCalROC(k);
1575 Int_t nchannels = calROC->GetNchannels();
1576 for(Int_t ch = 0; ch < nchannels; ch++){
1577 calROC->SetValue(ch,0.0);
1583 Int_t loop = (Int_t) vectorFit->GetEntriesFast();
1584 if(loop != 540) AliInfo("The Vector Fit is not complete!");
1585 Int_t detector = -1;
1586 Float_t value = 0.0;
1588 for (Int_t k = 0; k < loop; k++) {
1589 detector = ((AliTRDFitInfo *) vectorFit->At(k))->GetDetector();
1590 AliTRDCalROC *calROC = object->GetCalROC(detector);
1591 Float_t min = detobject->GetValue(detector);
1592 Int_t rowMax = calROC->GetNrows();
1593 Int_t colMax = calROC->GetNcols();
1594 for (Int_t row = 0; row < rowMax; row++) {
1595 for (Int_t col = 0; col < colMax; col++) {
1596 value = ((AliTRDFitInfo *) vectorFit->At(k))->GetCoef()[(Int_t)(col*rowMax+row)];
1597 calROC->SetValue(col,row,value-min);
1605 //_____________________________________________________________________________
1606 TObject *AliTRDCalibraFit::CreatePadObjectPRF(TObjArray *vectorFit)
1609 // It Creates the AliTRDCalPad object from AliTRDFitInfo
1610 // This object has to be written in the database
1613 // Create the DetObject
1614 AliTRDCalPad *object = new AliTRDCalPad("PRFWidth","PRFWidth");
1616 Int_t loop = (Int_t) vectorFit->GetEntriesFast();
1617 if(loop != 540) AliInfo("The Vector Fit is not complete!");
1618 Int_t detector = -1;
1619 Float_t value = 0.0;
1621 for (Int_t k = 0; k < loop; k++) {
1622 detector = ((AliTRDFitInfo *) vectorFit->At(k))->GetDetector();
1623 AliTRDCalROC *calROC = object->GetCalROC(detector);
1624 Int_t rowMax = calROC->GetNrows();
1625 Int_t colMax = calROC->GetNcols();
1626 for (Int_t row = 0; row < rowMax; row++) {
1627 for (Int_t col = 0; col < colMax; col++) {
1628 value = ((AliTRDFitInfo *) vectorFit->At(k))->GetCoef()[(Int_t)(col*rowMax+row)];
1629 calROC->SetValue(col,row,TMath::Abs(value));
1637 //_____________________________________________________________________________
1638 AliTRDCalDet *AliTRDCalibraFit::MakeOutliersStatDet(TObjArray *vectorFit, const char *name, Double_t &mean)
1641 // It Creates the AliTRDCalDet object from AliTRDFitInfo
1642 // 0 successful fit 1 not successful fit
1643 // mean is the mean value over the successful fit
1644 // do not use it for t0: no meaning
1647 // Create the CalObject
1648 AliTRDCalDet *object = new AliTRDCalDet(name,name);
1652 Int_t loop = (Int_t) vectorFit->GetEntriesFast();
1654 AliInfo("The Vector Fit is not complete! We initialise all outliers");
1655 for(Int_t k = 0; k < 540; k++){
1656 object->SetValue(k,1.0);
1659 Int_t detector = -1;
1660 Float_t value = 0.0;
1662 for (Int_t k = 0; k < loop; k++) {
1663 detector = ((AliTRDFitInfo *) vectorFit->At(k))->GetDetector();
1664 value = ((AliTRDFitInfo *) vectorFit->At(k))->GetCoef()[0];
1665 if(value <= 0) object->SetValue(detector,1.0);
1667 object->SetValue(detector,0.0);
1672 if(count > 0) mean /= count;
1675 //_____________________________________________________________________________
1676 TObject *AliTRDCalibraFit::MakeOutliersStatPad(TObjArray *vectorFit, const char *name, Double_t &mean)
1679 // It Creates the AliTRDCalPad object from AliTRDFitInfo
1680 // 0 not successful fit 1 successful fit
1681 // mean mean value over the successful fit
1684 // Create the CalObject
1685 AliTRDCalPad *object = new AliTRDCalPad(name,name);
1689 Int_t loop = (Int_t) vectorFit->GetEntriesFast();
1691 AliInfo("The Vector Fit is not complete! We initialise all outliers");
1692 for(Int_t k = 0; k < 540; k++){
1693 AliTRDCalROC *calROC = object->GetCalROC(k);
1694 Int_t nchannels = calROC->GetNchannels();
1695 for(Int_t ch = 0; ch < nchannels; ch++){
1696 calROC->SetValue(ch,1.0);
1700 Int_t detector = -1;
1701 Float_t value = 0.0;
1703 for (Int_t k = 0; k < loop; k++) {
1704 detector = ((AliTRDFitInfo *) vectorFit->At(k))->GetDetector();
1705 AliTRDCalROC *calROC = object->GetCalROC(detector);
1706 Int_t nchannels = calROC->GetNchannels();
1707 for (Int_t ch = 0; ch < nchannels; ch++) {
1708 value = ((AliTRDFitInfo *) vectorFit->At(k))->GetCoef()[ch];
1709 if(value <= 0) calROC->SetValue(ch,1.0);
1711 calROC->SetValue(ch,0.0);
1717 if(count > 0) mean /= count;
1720 //_____________________________________________________________________________
1721 void AliTRDCalibraFit::SetPeriodeFitPH(Int_t periodeFitPH)
1724 // Set FitPH if 1 then each detector will be fitted
1727 if (periodeFitPH > 0) {
1728 fFitPHPeriode = periodeFitPH;
1731 AliInfo("periodeFitPH must be higher than 0!");
1735 //_____________________________________________________________________________
1736 void AliTRDCalibraFit::SetBeginFitCharge(Float_t beginFitCharge)
1739 // The fit of the deposited charge distribution begins at
1740 // histo->Mean()/beginFitCharge
1741 // You can here set beginFitCharge
1744 if (beginFitCharge > 0) {
1745 fBeginFitCharge = beginFitCharge;
1748 AliInfo("beginFitCharge must be strict positif!");
1753 //_____________________________________________________________________________
1754 void AliTRDCalibraFit::SetT0Shift(Float_t t0Shift)
1757 // The t0 calculated with the maximum positif slope is shift from t0Shift
1758 // You can here set t0Shift
1765 AliInfo("t0Shift must be strict positif!");
1770 //_____________________________________________________________________________
1771 void AliTRDCalibraFit::SetRangeFitPRF(Float_t rangeFitPRF)
1774 // The fit of the PRF is from -rangeFitPRF to rangeFitPRF
1775 // You can here set rangeFitPRF
1778 if ((rangeFitPRF > 0) &&
1779 (rangeFitPRF <= 1.5)) {
1780 fRangeFitPRF = rangeFitPRF;
1783 AliInfo("rangeFitPRF must be between 0 and 1.0");
1788 //_____________________________________________________________________________
1789 void AliTRDCalibraFit::SetMinEntries(Int_t minEntries)
1792 // Minimum entries for fitting
1795 if (minEntries > 0) {
1796 fMinEntries = minEntries;
1799 AliInfo("fMinEntries must be >= 0.");
1804 //_____________________________________________________________________________
1805 void AliTRDCalibraFit::SetRebin(Short_t rebin)
1808 // Rebin with rebin time less bins the Ch histo
1809 // You can set here rebin that should divide the number of bins of CH histo
1814 AliInfo("You have to be sure that fRebin divides fNumberBinCharge used!");
1817 AliInfo("You have to choose a positiv value!");
1821 //_____________________________________________________________________________
1822 Bool_t AliTRDCalibraFit::FillVectorFit()
1825 // For the Fit functions fill the vector Fit
1828 AliTRDFitInfo *fitInfo = new AliTRDFitInfo();
1831 if (GetChamber(fCountDet) == 2) {
1838 //printf("For the detector %d , ntotal %d and fCoefCH[0] %f\n",countdet,ntotal,fCoefCH[0]);
1839 Float_t *coef = new Float_t[ntotal];
1840 for (Int_t i = 0; i < ntotal; i++) {
1841 coef[i] = fCurrentCoefDetector[i];
1844 Int_t detector = fCountDet;
1846 fitInfo->SetCoef(coef);
1847 fitInfo->SetDetector(detector);
1848 fVectorFit.Add((TObject *) fitInfo);
1853 //_____________________________________________________________________________
1854 Bool_t AliTRDCalibraFit::FillVectorFit2()
1857 // For the Fit functions fill the vector Fit
1860 AliTRDFitInfo *fitInfo = new AliTRDFitInfo();
1863 if (GetChamber(fCountDet) == 2) {
1870 //printf("For the detector %d , ntotal %d and fCoefCH[0] %f\n",countdet,ntotal,fCoefCH[0]);
1871 Float_t *coef = new Float_t[ntotal];
1872 for (Int_t i = 0; i < ntotal; i++) {
1873 coef[i] = fCurrentCoefDetector2[i];
1876 Int_t detector = fCountDet;
1878 fitInfo->SetCoef(coef);
1879 fitInfo->SetDetector(detector);
1880 fVectorFit2.Add((TObject *) fitInfo);
1885 //____________Functions for initialising the AliTRDCalibraFit in the code_________
1886 Bool_t AliTRDCalibraFit::InitFit(Int_t nbins, Int_t i)
1889 // Init the number of expected bins and fDect1[i] fDect2[i]
1892 gStyle->SetPalette(1);
1893 gStyle->SetOptStat(1111);
1894 gStyle->SetPadBorderMode(0);
1895 gStyle->SetCanvasColor(10);
1896 gStyle->SetPadLeftMargin(0.13);
1897 gStyle->SetPadRightMargin(0.01);
1899 // Mode groups of pads: the total number of bins!
1900 CalculNumberOfBinsExpected(i);
1902 // Quick verification that we have the good pad calibration mode!
1903 if (fNumberOfBinsExpected != nbins) {
1904 AliInfo("It doesn't correspond to the mode of pad group calibration!");
1908 // Security for fDebug 3 and 4
1909 if ((fDebugLevel >= 3) &&
1913 AliInfo("This detector doesn't exit!");
1917 // Determine fDet1 and fDet2 and set the fNfragZ and fNfragRphi for debug 3 and 4
1918 CalculDect1Dect2(i);
1923 //____________Functions for initialising the AliTRDCalibraFit in the code_________
1924 Bool_t AliTRDCalibraFit::InitFitCH()
1927 // Init the fVectorFitCH for normalisation
1928 // Init the histo for debugging
1933 fScaleFitFactor = 0.0;
1934 fCurrentCoefDetector = new Float_t[2304];
1935 for (Int_t k = 0; k < 2304; k++) {
1936 fCurrentCoefDetector[k] = 0.0;
1938 fVectorFit.SetName("gainfactorscoefficients");
1940 // fDebug == 0 nothing
1941 // fDebug == 1 and fFitVoir no histo
1942 if (fDebugLevel == 1) {
1943 if(!CheckFitVoir()) return kFALSE;
1945 //Get the CalDet object
1947 AliTRDcalibDB *cal = AliTRDcalibDB::Instance();
1949 AliInfo("Could not get calibDB");
1952 if(fCalDet) delete fCalDet;
1953 fCalDet = new AliTRDCalDet(*(cal->GetGainFactorDet()));
1956 Float_t devalue = 1.0;
1957 if(fCalDet) delete fCalDet;
1958 fCalDet = new AliTRDCalDet("ChamberGainFactor","GainFactor (detector value)");
1959 for(Int_t k = 0; k < 540; k++){
1960 fCalDet->SetValue(k,devalue);
1966 //____________Functions for initialising the AliTRDCalibraFit in the code_________
1967 Bool_t AliTRDCalibraFit::InitFitPH()
1970 // Init the arrays of results
1971 // Init the histos for debugging
1975 fVectorFit.SetName("driftvelocitycoefficients");
1976 fVectorFit2.SetName("t0coefficients");
1978 fCurrentCoefDetector = new Float_t[2304];
1979 for (Int_t k = 0; k < 2304; k++) {
1980 fCurrentCoefDetector[k] = 0.0;
1983 fCurrentCoefDetector2 = new Float_t[2304];
1984 for (Int_t k = 0; k < 2304; k++) {
1985 fCurrentCoefDetector2[k] = 0.0;
1988 //fDebug == 0 nothing
1989 // fDebug == 1 and fFitVoir no histo
1990 if (fDebugLevel == 1) {
1991 if(!CheckFitVoir()) return kFALSE;
1993 //Get the CalDet object
1995 AliTRDcalibDB *cal = AliTRDcalibDB::Instance();
1997 AliInfo("Could not get calibDB");
2000 if(fCalDet) delete fCalDet;
2001 if(fCalDet2) delete fCalDet2;
2002 fCalDet = new AliTRDCalDet(*(cal->GetVdriftDet()));
2003 fCalDet2 = new AliTRDCalDet(*(cal->GetT0Det()));
2006 Float_t devalue = 1.5;
2007 Float_t devalue2 = 0.0;
2008 if(fCalDet) delete fCalDet;
2009 if(fCalDet2) delete fCalDet2;
2010 fCalDet = new AliTRDCalDet("ChamberVdrift","TRD drift velocities (detector value)");
2011 fCalDet2 = new AliTRDCalDet("ChamberT0","T0 (detector value)");
2012 for(Int_t k = 0; k < 540; k++){
2013 fCalDet->SetValue(k,devalue);
2014 fCalDet2->SetValue(k,devalue2);
2019 //____________Functions for initialising the AliTRDCalibraFit in the code_________
2020 Bool_t AliTRDCalibraFit::InitFitPRF()
2023 // Init the calibration mode (Nz, Nrphi), the histograms for
2024 // debugging the fit methods if fDebug > 0,
2028 fVectorFit.SetName("prfwidthcoefficients");
2030 fCurrentCoefDetector = new Float_t[2304];
2031 for (Int_t k = 0; k < 2304; k++) {
2032 fCurrentCoefDetector[k] = 0.0;
2035 // fDebug == 0 nothing
2036 // fDebug == 1 and fFitVoir no histo
2037 if (fDebugLevel == 1) {
2038 if(!CheckFitVoir()) return kFALSE;
2042 //____________Functions for initialising the AliTRDCalibraFit in the code_________
2043 Bool_t AliTRDCalibraFit::InitFitLinearFitter()
2046 // Init the fCalDet, fVectorFit fCurrentCoefDetector
2051 fCurrentCoefDetector = new Float_t[2304];
2052 fCurrentCoefDetector2 = new Float_t[2304];
2053 for (Int_t k = 0; k < 2304; k++) {
2054 fCurrentCoefDetector[k] = 0.0;
2055 fCurrentCoefDetector2[k] = 0.0;
2058 //printf("test0\n");
2060 AliTRDcalibDB *cal = AliTRDcalibDB::Instance();
2062 AliInfo("Could not get calibDB");
2066 //Get the CalDet object
2068 if(fCalDet) delete fCalDet;
2069 if(fCalDet2) delete fCalDet2;
2070 fCalDet = new AliTRDCalDet(*(cal->GetVdriftDet()));
2071 //printf("test1\n");
2072 fCalDet2 = new AliTRDCalDet("lorentz angle tan","lorentz angle tan (detector value)");
2073 //printf("test2\n");
2074 for(Int_t k = 0; k < 540; k++){
2075 fCalDet2->SetValue(k,cal->GetOmegaTau(fCalDet->GetValue(k),-fMagneticField));
2077 //printf("test3\n");
2080 Float_t devalue = 1.5;
2081 Float_t devalue2 = cal->GetOmegaTau(1.5,-fMagneticField);
2082 if(fCalDet) delete fCalDet;
2083 if(fCalDet2) delete fCalDet2;
2084 //printf("test1\n");
2085 fCalDet = new AliTRDCalDet("ChamberVdrift","TRD drift velocities (detector value)");
2086 fCalDet2 = new AliTRDCalDet("lorentz angle tan","lorentz angle tan (detector value)");
2087 //printf("test2\n");
2088 for(Int_t k = 0; k < 540; k++){
2089 fCalDet->SetValue(k,devalue);
2090 fCalDet2->SetValue(k,devalue2);
2092 //printf("test3\n");
2097 //____________Functions for initialising the AliTRDCalibraFit in the code_________
2098 void AliTRDCalibraFit::InitfCountDetAndfCount(Int_t i)
2101 // Init the current detector where we are fCountDet and the
2102 // next fCount for the functions Fit...
2105 // Loop on the Xbins of ch!!
2106 fCountDet = -1; // Current detector
2107 fCount = 0; // To find the next detector
2110 if (fDebugLevel >= 3) {
2111 // Set countdet to the detector
2112 fCountDet = AliTRDgeometry::GetDetector(fDet[0],fDet[1],fDet[2]);
2113 // Set counter to write at the end of the detector
2115 // Get the right calib objects
2118 if(fDebugLevel == 1) {
2120 fCalibraMode->CalculXBins(fCountDet,i);
2121 while(fCalibraMode->GetXbins(i) <=fFitVoir){
2123 fCalibraMode->CalculXBins(fCountDet,i);
2125 fCount = fCalibraMode->GetXbins(i);
2127 // Determination of fNnZ, fNnRphi, fNfragZ and fNfragRphi
2128 fCalibraMode->ModePadCalibration((Int_t) GetChamber(fCountDet),i);
2129 fCalibraMode->ModePadFragmentation((Int_t) GetPlane(fCountDet)
2130 ,(Int_t) GetChamber(fCountDet)
2131 ,(Int_t) GetSector(fCountDet),i);
2134 //_______________________________________________________________________________
2135 void AliTRDCalibraFit::CalculNumberOfBinsExpected(Int_t i)
2138 // Calculate the number of bins expected (calibration groups)
2141 fNumberOfBinsExpected = 0;
2142 fCalibraMode->ModePadCalibration(2,i);
2143 fCalibraMode->ModePadFragmentation(0,2,0,i);
2144 fCalibraMode->SetDetChamb2(i);
2145 if (fDebugLevel > 1) {
2146 AliInfo(Form("For the chamber 2: %d",fCalibraMode->GetDetChamb2(i)));
2148 fNumberOfBinsExpected += 6 * 18 * fCalibraMode->GetDetChamb2(i);
2149 fCalibraMode->ModePadCalibration(0,i);
2150 fCalibraMode->ModePadFragmentation(0,0,0,i);
2151 fCalibraMode->SetDetChamb0(i);
2152 if (fDebugLevel > 1) {
2153 AliInfo(Form("For the other chamber 0: %d",fCalibraMode->GetDetChamb0(i)));
2155 fNumberOfBinsExpected += 6 * 4 * 18 * fCalibraMode->GetDetChamb0(i);
2158 //_______________________________________________________________________________
2159 void AliTRDCalibraFit::CalculDect1Dect2(Int_t i)
2162 // Calculate the range of fits
2167 if (fDebugLevel == 1) {
2171 if ((fDebugLevel == 2) || (fDebugLevel == 0)) {
2173 fDect2 = fNumberOfBinsExpected;
2175 if (fDebugLevel >= 3) {
2176 fCountDet = AliTRDgeometry::GetDetector(fDet[0],fDet[1],fDet[2]);
2177 fCalibraMode->CalculXBins(fCountDet,i);
2178 fDect1 = fCalibraMode->GetXbins(i);
2179 // Determination of fNnZ, fNnRphi, fNfragZ and fNfragRphi
2180 fCalibraMode->ModePadCalibration((Int_t) GetChamber(fCountDet),i);
2181 fCalibraMode->ModePadFragmentation((Int_t) GetPlane(fCountDet)
2182 ,(Int_t) GetChamber(fCountDet)
2183 ,(Int_t) GetSector(fCountDet),i);
2184 // Set for the next detector
2185 fDect2 = fDect1 + fCalibraMode->GetNfragZ(i)*fCalibraMode->GetNfragRphi(i);
2188 //_______________________________________________________________________________
2189 Bool_t AliTRDCalibraFit::CheckFitVoir()
2192 // Check if fFitVoir is in the range
2195 if (fFitVoir < fNumberOfBinsExpected) {
2196 AliInfo(Form("We will see the fit of the object %d",fFitVoir));
2199 AliInfo("fFitVoir is out of range of the histo!");
2204 //____________Functions for initialising the AliTRDCalibraFit in the code_________
2205 void AliTRDCalibraFit::UpdatefCountDetAndfCount(Int_t idect, Int_t i)
2208 // See if we are in a new detector and update the
2209 // variables fNfragZ and fNfragRphi if yes
2210 // Will never happen for only one detector (3 and 4)
2211 // Doesn't matter for 2
2213 if (fCount == idect) {
2214 // On en est au detector
2216 // Determination of fNnZ, fNnRphi, fNfragZ and fNfragRphi
2217 fCalibraMode->ModePadCalibration((Int_t) GetChamber(fCountDet),i);
2218 fCalibraMode->ModePadFragmentation((Int_t) GetPlane(fCountDet)
2219 ,(Int_t) GetChamber(fCountDet)
2220 ,(Int_t) GetSector(fCountDet),i);
2221 // Set for the next detector
2222 fCount += fCalibraMode->GetNfragZ(i)*fCalibraMode->GetNfragRphi(i);
2227 //____________Functions for initialising the AliTRDCalibraFit in the code_________
2228 void AliTRDCalibraFit::ReconstructFitRowMinRowMax(Int_t idect, Int_t i)
2231 // Reconstruct the min pad row, max pad row, min pad col and
2232 // max pad col of the calibration group for the Fit functions
2234 if (fDebugLevel != 1) {
2235 fCalibraMode->ReconstructionRowPadGroup((Int_t) (idect-(fCount-(fCalibraMode->GetNfragZ(i)*fCalibraMode->GetNfragRphi(i)))),i);
2238 //____________Functions for initialising the AliTRDCalibraFit in the code_________
2239 Bool_t AliTRDCalibraFit::NotEnoughStatisticCH(Int_t idect)
2242 // For the case where there are not enough entries in the histograms
2243 // of the calibration group, the value present in the choosen database
2244 // will be put. A negativ sign enables to know that a fit was not possible.
2247 if (fDebugLevel == 1) {
2248 AliInfo("The element has not enough statistic to be fitted");
2253 AliInfo(Form("The element %d in this detector %d has not enough statistic to be fitted"
2254 ,idect-(fCount-(fCalibraMode->GetNfragZ(0)*fCalibraMode->GetNfragRphi(0))),fCountDet));
2256 // Calcul the coef from the database choosen
2257 CalculChargeCoefMean(kFALSE);
2259 //chamber 2, not chamber 2
2261 if(GetChamber(fCountDet) == 2) factor = 12;
2264 // Fill the fCurrentCoefDetector with negative value to say: not fitted
2265 for (Int_t k = fCalibraMode->GetRowMin(0); k < fCalibraMode->GetRowMax(0); k++) {
2266 for (Int_t j = fCalibraMode->GetColMin(0); j < fCalibraMode->GetColMax(0); j++) {
2267 fCurrentCoefDetector[(Int_t)(j*factor+k)] = -TMath::Abs(fCurrentCoef[1]);
2271 //Put default value negative
2272 fCurrentCoef[0] = -TMath::Abs(fCurrentCoef[1]);
2273 fCurrentCoefE = 0.0;
2282 //____________Functions for initialising the AliTRDCalibraFit in the code_________
2283 Bool_t AliTRDCalibraFit::NotEnoughStatisticPH(Int_t idect)
2286 // For the case where there are not enough entries in the histograms
2287 // of the calibration group, the value present in the choosen database
2288 // will be put. A negativ sign enables to know that a fit was not possible.
2290 if (fDebugLevel == 1) {
2291 AliInfo("The element has not enough statistic to be fitted");
2295 AliInfo(Form("The element %d in this detector %d has not enough statistic to be fitted"
2296 ,idect-(fCount-(fCalibraMode->GetNfragZ(1)*fCalibraMode->GetNfragRphi(1))),fCountDet));
2298 CalculVdriftCoefMean();
2301 //chamber 2 and not chamber 2
2303 if(GetChamber(fCountDet) == 2) factor = 12;
2307 // Fill the fCurrentCoefDetector 2
2308 for (Int_t k = fCalibraMode->GetRowMin(1); k < fCalibraMode->GetRowMax(1); k++) {
2309 for (Int_t j = fCalibraMode->GetColMin(1); j < fCalibraMode->GetColMax(1); j++) {
2310 fCurrentCoefDetector[(Int_t)(j*factor+k)] = -TMath::Abs(fCurrentCoef[1]);
2311 fCurrentCoefDetector2[(Int_t)(j*factor+k)] = fCurrentCoef2[1];
2315 // Put the default value
2316 fCurrentCoef[0] = -TMath::Abs(fCurrentCoef[1]);
2317 fCurrentCoefE = 0.0;
2318 fCurrentCoef2[0] = fCurrentCoef2[1];
2319 fCurrentCoefE2 = 0.0;
2330 //____________Functions for initialising the AliTRDCalibraFit in the code_________
2331 Bool_t AliTRDCalibraFit::NotEnoughStatisticPRF(Int_t idect)
2334 // For the case where there are not enough entries in the histograms
2335 // of the calibration group, the value present in the choosen database
2336 // will be put. A negativ sign enables to know that a fit was not possible.
2339 if (fDebugLevel == 1) {
2340 AliInfo("The element has not enough statistic to be fitted");
2344 AliInfo(Form("The element %d in this detector %d has not enough statistic to be fitted"
2345 ,idect-(fCount-(fCalibraMode->GetNfragZ(2)*fCalibraMode->GetNfragRphi(2))),fCountDet));
2347 CalculPRFCoefMean();
2349 // chamber 2 and not chamber 2
2351 if(GetChamber(fCountDet) == 2) factor = 12;
2355 // Fill the fCurrentCoefDetector
2356 for (Int_t k = fCalibraMode->GetRowMin(2); k < fCalibraMode->GetRowMax(2); k++) {
2357 for (Int_t j = fCalibraMode->GetColMin(2); j < fCalibraMode->GetColMax(2); j++) {
2358 fCurrentCoefDetector[(Int_t)(j*factor+k)] = -fCurrentCoef[1];
2362 // Put the default value
2363 fCurrentCoef[0] = -fCurrentCoef[1];
2364 fCurrentCoefE = 0.0;
2372 //____________Functions for initialising the AliTRDCalibraFit in the code_________
2373 Bool_t AliTRDCalibraFit::NotEnoughStatisticLinearFitter()
2376 // For the case where there are not enough entries in the histograms
2377 // of the calibration group, the value present in the choosen database
2378 // will be put. A negativ sign enables to know that a fit was not possible.
2381 // Calcul the coef from the database choosen
2382 CalculVdriftLorentzCoef();
2385 if(GetChamber(fCountDet) == 2) factor = 1728;
2389 // Fill the fCurrentCoefDetector
2390 for (Int_t k = 0; k < factor; k++) {
2391 fCurrentCoefDetector[k] = -TMath::Abs(fCurrentCoef[1]);
2392 // should be negative
2393 fCurrentCoefDetector2[k] = +TMath::Abs(fCurrentCoef2[1]);
2397 //Put default opposite sign
2398 fCurrentCoef[0] = -TMath::Abs(fCurrentCoef[1]);
2399 fCurrentCoefE = 0.0;
2400 fCurrentCoef2[0] = +TMath::Abs(fCurrentCoef2[1]);
2401 fCurrentCoefE2 = 0.0;
2403 FillFillLinearFitter();
2408 //____________Functions for initialising the AliTRDCalibraFit in the code_________
2409 Bool_t AliTRDCalibraFit::FillInfosFitCH(Int_t idect)
2412 // Fill the coefficients found with the fits or other
2413 // methods from the Fit functions
2416 if (fDebugLevel != 1) {
2419 if(GetChamber(fCountDet) == 2) factor = 12;
2422 for (Int_t k = fCalibraMode->GetRowMin(0); k < fCalibraMode->GetRowMax(0); k++) {
2423 for (Int_t j = fCalibraMode->GetColMin(0); j < fCalibraMode->GetColMax(0); j++) {
2424 fCurrentCoefDetector[(Int_t)(j*factor+k)] = fCurrentCoef[0];
2435 //____________Functions for initialising the AliTRDCalibraFit in the code_________
2436 Bool_t AliTRDCalibraFit::FillInfosFitPH(Int_t idect)
2439 // Fill the coefficients found with the fits or other
2440 // methods from the Fit functions
2443 if (fDebugLevel != 1) {
2446 if(GetChamber(fCountDet) == 2) factor = 12;
2449 for (Int_t k = fCalibraMode->GetRowMin(1); k < fCalibraMode->GetRowMax(1); k++) {
2450 for (Int_t j = fCalibraMode->GetColMin(1); j < fCalibraMode->GetColMax(1); j++) {
2451 fCurrentCoefDetector[(Int_t)(j*factor+k)] = fCurrentCoef[0];
2452 fCurrentCoefDetector2[(Int_t)(j*factor+k)] = fCurrentCoef2[0];
2459 //____________Functions for initialising the AliTRDCalibraFit in the code_________
2460 Bool_t AliTRDCalibraFit::FillInfosFitPRF(Int_t idect)
2463 // Fill the coefficients found with the fits or other
2464 // methods from the Fit functions
2467 if (fDebugLevel != 1) {
2470 if(GetChamber(fCountDet) == 2) factor = 12;
2473 // Pointer to the branch
2474 for (Int_t k = fCalibraMode->GetRowMin(2); k < fCalibraMode->GetRowMax(2); k++) {
2475 for (Int_t j = fCalibraMode->GetColMin(2); j < fCalibraMode->GetColMax(2); j++) {
2476 fCurrentCoefDetector[(Int_t)(j*factor+k)] = fCurrentCoef[0];
2485 //____________Functions for initialising the AliTRDCalibraFit in the code_________
2486 Bool_t AliTRDCalibraFit::FillInfosFitLinearFitter()
2489 // Fill the coefficients found with the fits or other
2490 // methods from the Fit functions
2494 if(GetChamber(fCountDet) == 2) factor = 1728;
2497 // Pointer to the branch
2498 for (Int_t k = 0; k < factor; k++) {
2499 fCurrentCoefDetector[k] = fCurrentCoef[0];
2500 fCurrentCoefDetector2[k] = fCurrentCoef2[0];
2503 FillFillLinearFitter();
2508 //________________________________________________________________________________
2509 void AliTRDCalibraFit::FillFillCH(Int_t idect)
2512 // DebugStream and fVectorFit
2515 // End of one detector
2516 if ((idect == (fCount-1))) {
2519 for (Int_t k = 0; k < 2304; k++) {
2520 fCurrentCoefDetector[k] = 0.0;
2524 if(fDebugLevel > 1){
2526 if ( !fDebugStreamer ) {
2528 TDirectory *backup = gDirectory;
2529 fDebugStreamer = new TTreeSRedirector("TRDDebugFit.root");
2530 if ( backup ) backup->cd(); //we don't want to be cd'd to the debug streamer
2533 Int_t detector = fCountDet;
2534 Int_t caligroup = idect;
2535 Short_t rowmin = fCalibraMode->GetRowMin(0);
2536 Short_t rowmax = fCalibraMode->GetRowMax(0);
2537 Short_t colmin = fCalibraMode->GetColMin(0);
2538 Short_t colmax = fCalibraMode->GetColMax(0);
2539 Float_t gf = fCurrentCoef[0];
2540 Float_t gfs = fCurrentCoef[1];
2541 Float_t gfE = fCurrentCoefE;
2543 (*fDebugStreamer) << "GAIN" <<
2544 "detector=" << detector <<
2545 "caligroup=" << caligroup <<
2546 "rowmin=" << rowmin <<
2547 "rowmax=" << rowmax <<
2548 "colmin=" << colmin <<
2549 "colmax=" << colmax <<
2557 //________________________________________________________________________________
2558 void AliTRDCalibraFit::FillFillPH(Int_t idect)
2561 // DebugStream and fVectorFit and fVectorFit2
2564 // End of one detector
2565 if ((idect == (fCount-1))) {
2569 for (Int_t k = 0; k < 2304; k++) {
2570 fCurrentCoefDetector[k] = 0.0;
2571 fCurrentCoefDetector2[k] = 0.0;
2575 if(fDebugLevel > 1){
2577 if ( !fDebugStreamer ) {
2579 TDirectory *backup = gDirectory;
2580 fDebugStreamer = new TTreeSRedirector("TRDDebugFit.root");
2581 if ( backup ) backup->cd(); //we don't want to be cd'd to the debug streamer
2585 Int_t detector = fCountDet;
2586 Int_t caligroup = idect;
2587 Short_t rowmin = fCalibraMode->GetRowMin(1);
2588 Short_t rowmax = fCalibraMode->GetRowMax(1);
2589 Short_t colmin = fCalibraMode->GetColMin(1);
2590 Short_t colmax = fCalibraMode->GetColMax(1);
2591 Float_t vf = fCurrentCoef[0];
2592 Float_t vs = fCurrentCoef[1];
2593 Float_t vfE = fCurrentCoefE;
2594 Float_t t0f = fCurrentCoef2[0];
2595 Float_t t0s = fCurrentCoef2[1];
2596 Float_t t0E = fCurrentCoefE2;
2600 (* fDebugStreamer) << "PH"<<
2601 "detector="<<detector<<
2602 "caligroup="<<caligroup<<
2617 //________________________________________________________________________________
2618 void AliTRDCalibraFit::FillFillPRF(Int_t idect)
2621 // DebugStream and fVectorFit
2624 // End of one detector
2625 if ((idect == (fCount-1))) {
2628 for (Int_t k = 0; k < 2304; k++) {
2629 fCurrentCoefDetector[k] = 0.0;
2634 if(fDebugLevel > 1){
2636 if ( !fDebugStreamer ) {
2638 TDirectory *backup = gDirectory;
2639 fDebugStreamer = new TTreeSRedirector("TRDDebugFit.root");
2640 if ( backup ) backup->cd(); //we don't want to be cd'd to the debug streamer
2643 Int_t detector = fCountDet;
2644 Int_t plane = GetPlane(fCountDet);
2645 Int_t caligroup = idect;
2646 Short_t rowmin = fCalibraMode->GetRowMin(2);
2647 Short_t rowmax = fCalibraMode->GetRowMax(2);
2648 Short_t colmin = fCalibraMode->GetColMin(2);
2649 Short_t colmax = fCalibraMode->GetColMax(2);
2650 Float_t widf = fCurrentCoef[0];
2651 Float_t wids = fCurrentCoef[1];
2652 Float_t widfE = fCurrentCoefE;
2654 (* fDebugStreamer) << "PRF"<<
2655 "detector="<<detector<<
2657 "caligroup="<<caligroup<<
2669 //________________________________________________________________________________
2670 void AliTRDCalibraFit::FillFillLinearFitter()
2673 // DebugStream and fVectorFit
2676 // End of one detector
2682 for (Int_t k = 0; k < 2304; k++) {
2683 fCurrentCoefDetector[k] = 0.0;
2684 fCurrentCoefDetector2[k] = 0.0;
2688 if(fDebugLevel > 1){
2690 if ( !fDebugStreamer ) {
2692 TDirectory *backup = gDirectory;
2693 fDebugStreamer = new TTreeSRedirector("TRDDebugFit.root");
2694 if ( backup ) backup->cd(); //we don't want to be cd'd to the debug streamer
2697 //Debug: comparaison of the different methods (okey for first time but not for iterative procedure)
2698 AliTRDpadPlane *padplane = fGeo->GetPadPlane(GetPlane(fCountDet),GetChamber(fCountDet));
2699 Float_t rowmd = (padplane->GetRow0()+padplane->GetRowEnd())/2.;
2700 Float_t r = AliTRDgeometry::GetTime0(GetPlane(fCountDet));
2701 Float_t tiltangle = padplane->GetTiltingAngle();
2702 Int_t detector = fCountDet;
2703 Int_t chamber = GetChamber(fCountDet);
2704 Int_t plane = GetChamber(fCountDet);
2705 Float_t vf = fCurrentCoef[0];
2706 Float_t vs = fCurrentCoef[1];
2707 Float_t vfE = fCurrentCoefE;
2708 Float_t lorentzangler = fCurrentCoef2[0];
2709 Float_t Elorentzangler = fCurrentCoefE2;
2710 Float_t lorentzangles = fCurrentCoef2[1];
2712 (* fDebugStreamer) << "LinearFitter"<<
2713 "detector="<<detector<<
2714 "chamber="<<chamber<<
2718 "tiltangle="<<tiltangle<<
2722 "lorentzangler="<<lorentzangler<<
2723 "Elorentzangler="<<Elorentzangler<<
2724 "lorentzangles="<<lorentzangles<<
2730 //____________Calcul Coef Mean_________________________________________________
2732 //_____________________________________________________________________________
2733 Bool_t AliTRDCalibraFit::CalculT0CoefMean()
2736 // For the detector Dect calcul the mean time 0
2737 // for the calibration group idect from the choosen database
2740 fCurrentCoef2[1] = 0.0;
2741 if(fDebugLevel != 1){
2742 if((fCalibraMode->GetNz(1) > 0) ||
2743 (fCalibraMode->GetNrphi(1) > 0)) {
2744 for (Int_t row = fCalibraMode->GetRowMin(1); row < fCalibraMode->GetRowMax(1); row++) {
2745 for (Int_t col = fCalibraMode->GetColMin(1); col < fCalibraMode->GetColMax(1); col++) {
2746 fCurrentCoef2[1] += (Float_t) (fCalROC2->GetValue(col,row)+fCalDet2->GetValue(fCountDet));
2749 fCurrentCoef2[1] = fCurrentCoef2[1] / ((fCalibraMode->GetColMax(1)-fCalibraMode->GetColMin(1))*(fCalibraMode->GetRowMax(1)-fCalibraMode->GetRowMin(1)));
2753 fCurrentCoef2[1] = fCalDet2->GetValue(fCountDet);
2756 for(Int_t row = 0; row < fGeo->GetRowMax(GetPlane(fCountDet),GetChamber(fCountDet),GetSector(fCountDet)); row++){
2757 for(Int_t col = 0; col < fGeo->GetColMax(GetPlane(fCountDet)); col++){
2758 fCurrentCoef2[1] += (Float_t) (fCalROC2->GetValue(col,row)+fCalDet2->GetValue(fCountDet));
2761 fCurrentCoef2[1] = fCurrentCoef2[1] / ((fGeo->GetRowMax(GetPlane(fCountDet),GetChamber(fCountDet),GetSector(fCountDet)))*(fGeo->GetColMax(GetPlane(fCountDet))));
2768 //_____________________________________________________________________________
2769 Bool_t AliTRDCalibraFit::CalculChargeCoefMean(Bool_t vrai)
2772 // For the detector Dect calcul the mean gain factor
2773 // for the calibration group idect from the choosen database
2776 fCurrentCoef[1] = 0.0;
2777 if(fDebugLevel != 1){
2778 if ((fCalibraMode->GetNz(0) > 0) ||
2779 (fCalibraMode->GetNrphi(0) > 0)) {
2780 for (Int_t row = fCalibraMode->GetRowMin(0); row < fCalibraMode->GetRowMax(0); row++) {
2781 for (Int_t col = fCalibraMode->GetColMin(0); col < fCalibraMode->GetColMax(0); col++) {
2782 fCurrentCoef[1] += (Float_t) (fCalROC->GetValue(col,row)*fCalDet->GetValue(fCountDet));
2783 if (vrai) fScaleFitFactor += (Float_t) (fCalROC->GetValue(col,row)*fCalDet->GetValue(fCountDet));
2786 fCurrentCoef[1] = fCurrentCoef[1] / ((fCalibraMode->GetColMax(0)-fCalibraMode->GetColMin(0))*(fCalibraMode->GetRowMax(0)-fCalibraMode->GetRowMin(0)));
2790 fCurrentCoef[1] = (Float_t) fCalDet->GetValue(fCountDet);
2791 if (vrai) fScaleFitFactor += ((Float_t) fCalDet->GetValue(fCountDet))*(fCalibraMode->GetColMax(0)-fCalibraMode->GetColMin(0))*(fCalibraMode->GetRowMax(0)-fCalibraMode->GetRowMin(0));
2796 //_____________________________________________________________________________
2797 Bool_t AliTRDCalibraFit::CalculPRFCoefMean()
2800 // For the detector Dect calcul the mean sigma of pad response
2801 // function for the calibration group idect from the choosen database
2804 fCurrentCoef[1] = 0.0;
2805 if(fDebugLevel != 1){
2806 for (Int_t row = fCalibraMode->GetRowMin(2); row < fCalibraMode->GetRowMax(2); row++) {
2807 for (Int_t col = fCalibraMode->GetColMin(2); col < fCalibraMode->GetColMax(2); col++) {
2808 fCurrentCoef[1] += (Float_t) fCalROC->GetValue(col,row);
2811 fCurrentCoef[1] = fCurrentCoef[1] / ((fCalibraMode->GetColMax(2)-fCalibraMode->GetColMin(2))*(fCalibraMode->GetRowMax(2)-fCalibraMode->GetRowMin(2)));
2815 //_____________________________________________________________________________
2816 Bool_t AliTRDCalibraFit::CalculVdriftCoefMean()
2819 // For the detector dect calcul the mean drift velocity for the
2820 // calibration group idect from the choosen database
2823 fCurrentCoef[1] = 0.0;
2824 if(fDebugLevel != 1){
2825 if ((fCalibraMode->GetNz(1) > 0) ||
2826 (fCalibraMode->GetNrphi(1) > 0)) {
2827 for (Int_t row = fCalibraMode->GetRowMin(1); row < fCalibraMode->GetRowMax(1); row++) {
2828 for (Int_t col = fCalibraMode->GetColMin(1); col < fCalibraMode->GetColMax(1); col++) {
2829 fCurrentCoef[1] += (Float_t) (fCalROC->GetValue(col,row)*fCalDet->GetValue(fCountDet));
2832 fCurrentCoef[1] = fCurrentCoef[1] / ((fCalibraMode->GetColMax(1)-fCalibraMode->GetColMin(1))*(fCalibraMode->GetRowMax(1)-fCalibraMode->GetRowMin(1)));
2836 fCurrentCoef[1] = (Float_t) fCalDet->GetValue(fCountDet);
2841 //_____________________________________________________________________________
2842 Bool_t AliTRDCalibraFit::CalculVdriftLorentzCoef()
2845 // For the detector fCountDet, mean drift velocity and tan lorentzangle
2848 fCurrentCoef[1] = fCalDet->GetValue(fCountDet);
2849 fCurrentCoef2[1] = fCalDet2->GetValue(fCountDet);
2853 //_____________________________________________________________________________
2854 Float_t AliTRDCalibraFit::GetPRFDefault(Int_t plane) const
2857 // Default width of the PRF if there is no database as reference
2862 //case 0: return 0.515;
2863 //case 1: return 0.502;
2864 //case 2: return 0.491;
2865 //case 3: return 0.481;
2866 //case 4: return 0.471;
2867 //case 5: return 0.463;
2868 //default: return 0.0;
2871 case 0: return 0.538429;
2872 case 1: return 0.524302;
2873 case 2: return 0.511591;
2874 case 3: return 0.500140;
2875 case 4: return 0.489821;
2876 case 5: return 0.480524;
2877 default: return 0.0;
2880 //________________________________________________________________________________
2881 void AliTRDCalibraFit::SetCalROC(Int_t i)
2884 // Set the calib object for fCountDet
2887 Float_t value = 0.0;
2889 //Get the CalDet object
2891 AliTRDcalibDB *cal = AliTRDcalibDB::Instance();
2893 AliInfo("Could not get calibDB");
2899 if(fCalROC) delete fCalROC;
2900 fCalROC = new AliTRDCalROC(*(cal->GetGainFactorROC(fCountDet)));
2903 if(fCalROC) delete fCalROC;
2904 if(fCalROC2) delete fCalROC2;
2905 fCalROC = new AliTRDCalROC(*(cal->GetVdriftROC(fCountDet)));
2906 fCalROC2 = new AliTRDCalROC(*(cal->GetT0ROC(fCountDet)));
2909 if(fCalROC) delete fCalROC;
2910 fCalROC = new AliTRDCalROC(*(cal->GetPRFROC(fCountDet))); break;
2918 if(fCalROC) delete fCalROC;
2919 fCalROC = new AliTRDCalROC(GetPlane(fCountDet),GetChamber(fCountDet));
2920 for(Int_t k = 0; k < fCalROC->GetNchannels(); k++){
2921 fCalROC->SetValue(k,1.0);
2925 if(fCalROC) delete fCalROC;
2926 if(fCalROC2) delete fCalROC2;
2927 fCalROC = new AliTRDCalROC(GetPlane(fCountDet),GetChamber(fCountDet));
2928 fCalROC2 = new AliTRDCalROC(GetPlane(fCountDet),GetChamber(fCountDet));
2929 for(Int_t k = 0; k < fCalROC->GetNchannels(); k++){
2930 fCalROC->SetValue(k,1.0);
2931 fCalROC2->SetValue(k,0.0);
2935 if(fCalROC) delete fCalROC;
2936 value = GetPRFDefault(GetPlane(fCountDet));
2937 fCalROC = new AliTRDCalROC(GetPlane(fCountDet),GetChamber(fCountDet));
2938 for(Int_t k = 0; k < fCalROC->GetNchannels(); k++){
2939 fCalROC->SetValue(k,value);
2947 //____________Fit Methods______________________________________________________
2949 //_____________________________________________________________________________
2950 void AliTRDCalibraFit::FitPente(TH1* projPH)
2953 // Slope methode for the drift velocity
2957 const Float_t kDrWidth = AliTRDgeometry::DrThick();
2964 fCurrentCoefE = 0.0;
2965 fCurrentCoefE2 = 0.0;
2966 fCurrentCoef[0] = 0.0;
2967 fCurrentCoef2[0] = 0.0;
2968 TLine *line = new TLine();
2971 TAxis *xpph = projPH->GetXaxis();
2972 Int_t nbins = xpph->GetNbins();
2973 Double_t lowedge = xpph->GetBinLowEdge(1);
2974 Double_t upedge = xpph->GetBinUpEdge(xpph->GetNbins());
2975 Double_t widbins = (upedge - lowedge) / nbins;
2976 Double_t limit = upedge + 0.5 * widbins;
2979 // Beginning of the signal
2980 TH1D *pentea = new TH1D("pentea","pentea",projPH->GetNbinsX(),0,(Float_t) limit);
2981 for (Int_t k = 1; k < projPH->GetNbinsX(); k++) {
2982 pentea->SetBinContent(k,(Double_t) (projPH->GetBinContent(k+1) - projPH->GetBinContent(k)));
2984 binmax = (Int_t) pentea->GetMaximumBin();
2987 AliInfo("Put the binmax from 1 to 2 to enable the fit");
2989 if (binmax >= nbins) {
2992 AliInfo("Put the binmax from nbins-1 to nbins-2 to enable the fit");
2994 pentea->Fit("pol2","0MR","",TMath::Max(pentea->GetBinCenter(binmax-1),0.0),pentea->GetBinCenter(binmax+1));
2995 Float_t l3P1am = pentea->GetFunction("pol2")->GetParameter(1);
2996 Float_t l3P2am = pentea->GetFunction("pol2")->GetParameter(2);
2997 Float_t l3P1amE = pentea->GetFunction("pol2")->GetParError(1);
2998 Float_t l3P2amE = pentea->GetFunction("pol2")->GetParError(2);
3000 fPhd[0] = -(l3P1am / (2 * l3P2am));
3003 if((l3P1am != 0.0) && (l3P2am != 0.0)){
3004 fCurrentCoefE2 = (l3P1amE/l3P1am + l3P2amE/l3P2am)*fPhd[0];
3007 // Amplification region
3010 for (Int_t kbin = 1; kbin < projPH->GetNbinsX(); kbin ++) {
3011 if (((projPH->GetBinContent(kbin+1) - projPH->GetBinContent(kbin)) <= 0.0) && (ju == 0) && (kbin > (fPhd[0]/widbins))) {
3018 AliInfo("Put the binmax from 1 to 2 to enable the fit");
3020 if (binmax >= nbins) {
3023 AliInfo("Put the binmax from nbins-1 to nbins-2 to enable the fit");
3025 projPH->Fit("pol2","0MR","",TMath::Max(projPH->GetBinCenter(binmax-1),0.0),projPH->GetBinCenter(binmax+1));
3026 Float_t l3P1amf = projPH->GetFunction("pol2")->GetParameter(1);
3027 Float_t l3P2amf = projPH->GetFunction("pol2")->GetParameter(2);
3028 Float_t l3P1amfE = projPH->GetFunction("pol2")->GetParError(1);
3029 Float_t l3P2amfE = projPH->GetFunction("pol2")->GetParError(2);
3031 fPhd[1] = -(l3P1amf / (2 * l3P2amf));
3033 if((l3P1amf != 0.0) && (l3P2amf != 0.0)){
3034 fCurrentCoefE = (l3P1amfE/l3P1amf + l3P2amfE/l3P2amf)*fPhd[1];
3037 fCurrentCoefE2 = fCurrentCoefE;
3040 TH1D *pente = new TH1D("pente","pente",projPH->GetNbinsX(),0,(Float_t) limit);
3041 for (Int_t k = TMath::Min(binmax+4,projPH->GetNbinsX()); k < projPH->GetNbinsX(); k++) {
3042 pente->SetBinContent(k,(Double_t) (projPH->GetBinContent(k+1) - projPH->GetBinContent(k)));
3045 if(pente->GetEntries() > 0) binmin = (Int_t) pente->GetMinimumBin();
3048 AliInfo("Put the binmax from 1 to 2 to enable the fit");
3050 if (binmin >= nbins) {
3053 AliInfo("Put the binmax from nbins-1 to nbins-2 to enable the fit");
3058 ,TMath::Max(pente->GetBinCenter(binmin-1), 0.0)
3059 ,TMath::Min(pente->GetBinCenter(binmin+1),(Double_t) limit));
3060 Float_t l3P1dr = pente->GetFunction("pol2")->GetParameter(1);
3061 Float_t l3P2dr = pente->GetFunction("pol2")->GetParameter(2);
3062 Float_t l3P1drE = pente->GetFunction("pol2")->GetParError(1);
3063 Float_t l3P2drE = pente->GetFunction("pol2")->GetParError(2);
3065 fPhd[2] = -(l3P1dr / (2 * l3P2dr));
3067 if((l3P1dr != 0.0) && (l3P2dr != 0.0)){
3068 fCurrentCoefE += (l3P1drE/l3P1dr + l3P2drE/l3P2dr)*fPhd[2];
3070 Float_t fPhdt0 = 0.0;
3071 if(fTakeTheMaxPH) fPhdt0 = fPhd[1];
3072 else fPhdt0 = fPhd[0];
3074 if ((fPhd[2] > fPhd[0]) &&
3075 (fPhd[2] > fPhd[1]) &&
3076 (fPhd[1] > fPhd[0]) &&
3078 fCurrentCoef[0] = (kDrWidth) / (fPhd[2]-fPhd[1]);
3079 fNumberFitSuccess++;
3081 if (fPhdt0 >= 0.0) {
3082 fCurrentCoef2[0] = (fPhdt0 - fT0Shift) / widbins;
3083 if (fCurrentCoef2[0] < -1.0) {
3084 fCurrentCoef2[0] = fCurrentCoef2[1];
3088 fCurrentCoef2[0] = fCurrentCoef2[1];
3093 fCurrentCoef[0] = -TMath::Abs(fCurrentCoef[1]);
3094 fCurrentCoef2[0] = fCurrentCoef2[1];
3097 if (fDebugLevel == 1) {
3098 TCanvas *cpentei = new TCanvas("cpentei","cpentei",50,50,600,800);
3101 line->SetLineColor(2);
3102 line->DrawLine(fPhd[0],0,fPhd[0],projPH->GetMaximum());
3103 line->DrawLine(fPhd[1],0,fPhd[1],projPH->GetMaximum());
3104 line->DrawLine(fPhd[2],0,fPhd[2],projPH->GetMaximum());
3105 AliInfo(Form("fPhd[0] (beginning of the signal): %f" ,(Float_t) fPhd[0]));
3106 AliInfo(Form("fPhd[1] (end of the amplification region): %f" ,(Float_t) fPhd[1]));
3107 AliInfo(Form("fPhd[2] (end of the drift region): %f" ,(Float_t) fPhd[2]));
3108 AliInfo(Form("fVriftCoef[1] (with only the drift region(default)): %f",(Float_t) fCurrentCoef[0]));
3109 TCanvas *cpentei2 = new TCanvas("cpentei2","cpentei2",50,50,600,800);
3112 TCanvas *cpentei3 = new TCanvas("cpentei3","cpentei3",50,50,600,800);
3121 //_____________________________________________________________________________
3122 void AliTRDCalibraFit::FitLagrangePoly(TH1* projPH)
3125 // Slope methode but with polynomes de Lagrange
3129 const Float_t kDrWidth = AliTRDgeometry::DrThick();
3132 Double_t *x = new Double_t[5];
3133 Double_t *y = new Double_t[5];
3148 fCurrentCoefE = 0.0;
3149 fCurrentCoefE2 = 1.0;
3150 fCurrentCoef[0] = 0.0;
3151 fCurrentCoef2[0] = 0.0;
3152 TLine *line = new TLine();
3153 TF1 * polynome = 0x0;
3154 TF1 * polynomea = 0x0;
3155 TF1 * polynomeb = 0x0;
3159 TAxis *xpph = projPH->GetXaxis();
3160 Int_t nbins = xpph->GetNbins();
3161 Double_t lowedge = xpph->GetBinLowEdge(1);
3162 Double_t upedge = xpph->GetBinUpEdge(xpph->GetNbins());
3163 Double_t widbins = (upedge - lowedge) / nbins;
3164 Double_t limit = upedge + 0.5 * widbins;
3169 // Beginning of the signal
3170 TH1D *pentea = new TH1D("pentea","pentea",projPH->GetNbinsX(),0,(Float_t) limit);
3171 for (Int_t k = 1; k < projPH->GetNbinsX(); k++) {
3172 pentea->SetBinContent(k,(Double_t) (projPH->GetBinContent(k+1) - projPH->GetBinContent(k)));
3175 binmax = (Int_t) pentea->GetMaximumBin();
3177 Double_t minnn = 0.0;
3178 Double_t maxxx = 0.0;
3180 Int_t kase = nbins-binmax;
3188 minnn = pentea->GetBinCenter(binmax-2);
3189 maxxx = pentea->GetBinCenter(binmax);
3190 x[0] = pentea->GetBinCenter(binmax-2);
3191 x[1] = pentea->GetBinCenter(binmax-1);
3192 x[2] = pentea->GetBinCenter(binmax);
3193 y[0] = pentea->GetBinContent(binmax-2);
3194 y[1] = pentea->GetBinContent(binmax-1);
3195 y[2] = pentea->GetBinContent(binmax);
3196 c = CalculPolynomeLagrange2(x,y);
3197 AliInfo("At the limit for beginning!");
3200 minnn = pentea->GetBinCenter(binmax-2);
3201 maxxx = pentea->GetBinCenter(binmax+1);
3202 x[0] = pentea->GetBinCenter(binmax-2);
3203 x[1] = pentea->GetBinCenter(binmax-1);
3204 x[2] = pentea->GetBinCenter(binmax);
3205 x[3] = pentea->GetBinCenter(binmax+1);
3206 y[0] = pentea->GetBinContent(binmax-2);
3207 y[1] = pentea->GetBinContent(binmax-1);
3208 y[2] = pentea->GetBinContent(binmax);
3209 y[3] = pentea->GetBinContent(binmax+1);
3210 c = CalculPolynomeLagrange3(x,y);
3218 minnn = pentea->GetBinCenter(binmax);
3219 maxxx = pentea->GetBinCenter(binmax+2);
3220 x[0] = pentea->GetBinCenter(binmax);
3221 x[1] = pentea->GetBinCenter(binmax+1);
3222 x[2] = pentea->GetBinCenter(binmax+2);
3223 y[0] = pentea->GetBinContent(binmax);
3224 y[1] = pentea->GetBinContent(binmax+1);
3225 y[2] = pentea->GetBinContent(binmax+2);
3226 c = CalculPolynomeLagrange2(x,y);
3229 minnn = pentea->GetBinCenter(binmax-1);
3230 maxxx = pentea->GetBinCenter(binmax+2);
3231 x[0] = pentea->GetBinCenter(binmax-1);
3232 x[1] = pentea->GetBinCenter(binmax);
3233 x[2] = pentea->GetBinCenter(binmax+1);
3234 x[3] = pentea->GetBinCenter(binmax+2);
3235 y[0] = pentea->GetBinContent(binmax-1);
3236 y[1] = pentea->GetBinContent(binmax);
3237 y[2] = pentea->GetBinContent(binmax+1);
3238 y[3] = pentea->GetBinContent(binmax+2);
3239 c = CalculPolynomeLagrange3(x,y);
3242 minnn = pentea->GetBinCenter(binmax-2);
3243 maxxx = pentea->GetBinCenter(binmax+2);
3244 x[0] = pentea->GetBinCenter(binmax-2);
3245 x[1] = pentea->GetBinCenter(binmax-1);
3246 x[2] = pentea->GetBinCenter(binmax);
3247 x[3] = pentea->GetBinCenter(binmax+1);
3248 x[4] = pentea->GetBinCenter(binmax+2);
3249 y[0] = pentea->GetBinContent(binmax-2);
3250 y[1] = pentea->GetBinContent(binmax-1);
3251 y[2] = pentea->GetBinContent(binmax);
3252 y[3] = pentea->GetBinContent(binmax+1);
3253 y[4] = pentea->GetBinContent(binmax+2);
3254 c = CalculPolynomeLagrange4(x,y);
3262 polynomeb = new TF1("polb","[0]+[1]*x+[2]*x*x+[3]*x*x*x+[4]*x*x*x*x",minnn,maxxx);
3263 polynomeb->SetParameters(c[0],c[1],c[2],c[3],c[4]);
3265 Double_t step = (maxxx-minnn)/10000;
3267 Double_t maxvalue = 0.0;
3268 Double_t placemaximum = minnn;
3269 for(Int_t o = 0; o < 10000; o++){
3270 if(o == 0) maxvalue = polynomeb->Eval(l);
3271 if(maxvalue < (polynomeb->Eval(l))){
3272 maxvalue = polynomeb->Eval(l);
3277 fPhd[0] = placemaximum;
3280 // Amplification region
3283 for (Int_t kbin = 1; kbin < projPH->GetNbinsX(); kbin ++) {
3284 if (((projPH->GetBinContent(kbin+1) - projPH->GetBinContent(kbin)) <= 0.0) && (ju == 0) && (kbin > (fPhd[0]/widbins))) {
3290 Double_t minn = 0.0;
3291 Double_t maxx = 0.0;
3303 Int_t kase1 = nbins - binmax;
3305 //Determination of minn and maxx
3306 //case binmax = nbins
3311 minn = projPH->GetBinCenter(binmax-2);
3312 maxx = projPH->GetBinCenter(binmax);
3313 x[0] = projPH->GetBinCenter(binmax-2);
3314 x[1] = projPH->GetBinCenter(binmax-1);
3315 x[2] = projPH->GetBinCenter(binmax);
3316 y[0] = projPH->GetBinContent(binmax-2);
3317 y[1] = projPH->GetBinContent(binmax-1);
3318 y[2] = projPH->GetBinContent(binmax);
3319 c = CalculPolynomeLagrange2(x,y);
3320 //AliInfo("At the limit for the drift!");
3323 minn = projPH->GetBinCenter(binmax-2);
3324 maxx = projPH->GetBinCenter(binmax+1);
3325 x[0] = projPH->GetBinCenter(binmax-2);
3326 x[1] = projPH->GetBinCenter(binmax-1);
3327 x[2] = projPH->GetBinCenter(binmax);
3328 x[3] = projPH->GetBinCenter(binmax+1);
3329 y[0] = projPH->GetBinContent(binmax-2);
3330 y[1] = projPH->GetBinContent(binmax-1);
3331 y[2] = projPH->GetBinContent(binmax);
3332 y[3] = projPH->GetBinContent(binmax+1);
3333 c = CalculPolynomeLagrange3(x,y);
3342 minn = projPH->GetBinCenter(binmax);
3343 maxx = projPH->GetBinCenter(binmax+2);
3344 x[0] = projPH->GetBinCenter(binmax);
3345 x[1] = projPH->GetBinCenter(binmax+1);
3346 x[2] = projPH->GetBinCenter(binmax+2);
3347 y[0] = projPH->GetBinContent(binmax);
3348 y[1] = projPH->GetBinContent(binmax+1);
3349 y[2] = projPH->GetBinContent(binmax+2);
3350 c = CalculPolynomeLagrange2(x,y);
3353 minn = projPH->GetBinCenter(binmax-1);
3354 maxx = projPH->GetBinCenter(binmax+2);
3355 x[0] = projPH->GetBinCenter(binmax-1);
3356 x[1] = projPH->GetBinCenter(binmax);
3357 x[2] = projPH->GetBinCenter(binmax+1);
3358 x[3] = projPH->GetBinCenter(binmax+2);
3359 y[0] = projPH->GetBinContent(binmax-1);
3360 y[1] = projPH->GetBinContent(binmax);
3361 y[2] = projPH->GetBinContent(binmax+1);
3362 y[3] = projPH->GetBinContent(binmax+2);
3363 c = CalculPolynomeLagrange3(x,y);
3366 minn = projPH->GetBinCenter(binmax-2);
3367 maxx = projPH->GetBinCenter(binmax+2);
3368 x[0] = projPH->GetBinCenter(binmax-2);
3369 x[1] = projPH->GetBinCenter(binmax-1);
3370 x[2] = projPH->GetBinCenter(binmax);
3371 x[3] = projPH->GetBinCenter(binmax+1);
3372 x[4] = projPH->GetBinCenter(binmax+2);
3373 y[0] = projPH->GetBinContent(binmax-2);
3374 y[1] = projPH->GetBinContent(binmax-1);
3375 y[2] = projPH->GetBinContent(binmax);
3376 y[3] = projPH->GetBinContent(binmax+1);
3377 y[4] = projPH->GetBinContent(binmax+2);
3378 c = CalculPolynomeLagrange4(x,y);
3385 polynomea = new TF1("pola","[0]+[1]*x+[2]*x*x+[3]*x*x*x+[4]*x*x*x*x",minn,maxx);
3386 polynomea->SetParameters(c[0],c[1],c[2],c[3],c[4]);
3388 Double_t step = (maxx-minn)/1000;
3390 Double_t maxvalue = 0.0;
3391 Double_t placemaximum = minn;
3392 for(Int_t o = 0; o < 1000; o++){
3393 if(o == 0) maxvalue = polynomea->Eval(l);
3394 if(maxvalue < (polynomea->Eval(l))){
3395 maxvalue = polynomea->Eval(l);
3400 fPhd[1] = placemaximum;
3404 TH1D *pente = new TH1D("pente","pente", projPH->GetNbinsX(),0,(Float_t) limit);
3405 for (Int_t k = TMath::Min(binmax+4, projPH->GetNbinsX()); k < projPH->GetNbinsX(); k++) {
3406 pente->SetBinContent(k,(Double_t) (projPH->GetBinContent(k+1) - projPH->GetBinContent(k)));
3409 if(pente->GetEntries() > 0) binmin = (Int_t) pente->GetMinimumBin();
3415 AliInfo("Put the binmax from 1 to 2 to enable the fit");
3419 if((projPH->GetBinContent(binmin)-projPH->GetBinError(binmin)) < (projPH->GetBinContent(binmin+1))) put = kFALSE;
3420 if((projPH->GetBinContent(binmin)+projPH->GetBinError(binmin)) > (projPH->GetBinContent(binmin-1))) put = kFALSE;
3434 Bool_t case1 = kFALSE;
3435 Bool_t case2 = kFALSE;
3436 Bool_t case4 = kFALSE;
3438 //Determination of min and max
3439 //case binmin <= nbins-3
3441 if((binmin <= (nbins-3)) && ((binmin-2) >= TMath::Min(binmax+4, projPH->GetNbinsX()))){
3442 min = pente->GetBinCenter(binmin-2);
3443 max = pente->GetBinCenter(binmin+2);
3444 x[0] = pente->GetBinCenter(binmin-2);
3445 x[1] = pente->GetBinCenter(binmin-1);
3446 x[2] = pente->GetBinCenter(binmin);
3447 x[3] = pente->GetBinCenter(binmin+1);
3448 x[4] = pente->GetBinCenter(binmin+2);
3449 y[0] = pente->GetBinContent(binmin-2);
3450 y[1] = pente->GetBinContent(binmin-1);
3451 y[2] = pente->GetBinContent(binmin);
3452 y[3] = pente->GetBinContent(binmin+1);
3453 y[4] = pente->GetBinContent(binmin+2);
3454 //Calcul the polynome de Lagrange
3455 c = CalculPolynomeLagrange4(x,y);
3457 if((pente->GetBinContent(binmin+2) <= pente->GetBinContent(binmin+1)) &&
3458 (pente->GetBinContent(binmin-2) <= pente->GetBinContent(binmin-1))) put = kFALSE;
3459 if(((binmin+3) <= (nbins-1)) &&
3460 (pente->GetBinContent(binmin+3) <= pente->GetBinContent(binmin+2)) &&
3461 ((binmin-3) >= TMath::Min(binmax+4, projPH->GetNbinsX())) &&
3462 (pente->GetBinContent(binmin-3) <= pente->GetBinContent(binmin-2))) put = kFALSE;
3463 if((pente->GetBinContent(binmin+2) <= pente->GetBinContent(binmin+1)) &&
3464 (pente->GetBinContent(binmin-2) > pente->GetBinContent(binmin-1))) case1 = kTRUE;
3465 if((pente->GetBinContent(binmin+2) > pente->GetBinContent(binmin+1)) &&
3466 (pente->GetBinContent(binmin-2) <= pente->GetBinContent(binmin-1))) case4 = kTRUE;
3468 //case binmin = nbins-2
3470 if(((binmin == (nbins-2)) && ((binmin-2) >= TMath::Min(binmax+4, projPH->GetNbinsX()))) ||
3472 min = pente->GetBinCenter(binmin-2);
3473 max = pente->GetBinCenter(binmin+1);
3474 x[0] = pente->GetBinCenter(binmin-2);
3475 x[1] = pente->GetBinCenter(binmin-1);
3476 x[2] = pente->GetBinCenter(binmin);
3477 x[3] = pente->GetBinCenter(binmin+1);
3478 y[0] = pente->GetBinContent(binmin-2);
3479 y[1] = pente->GetBinContent(binmin-1);
3480 y[2] = pente->GetBinContent(binmin);
3481 y[3] = pente->GetBinContent(binmin+1);
3482 //Calcul the polynome de Lagrange
3483 c = CalculPolynomeLagrange3(x,y);
3484 //richtung +: nothing
3486 if((pente->GetBinContent(binmin-2) <= pente->GetBinContent(binmin-1))) case2 = kTRUE;
3489 if(((binmin <= (nbins-3)) && ((binmin-1) == TMath::Min(binmax+4, projPH->GetNbinsX()))) ||
3491 min = pente->GetBinCenter(binmin-1);
3492 max = pente->GetBinCenter(binmin+2);
3493 x[0] = pente->GetBinCenter(binmin-1);
3494 x[1] = pente->GetBinCenter(binmin);
3495 x[2] = pente->GetBinCenter(binmin+1);
3496 x[3] = pente->GetBinCenter(binmin+2);
3497 y[0] = pente->GetBinContent(binmin-1);
3498 y[1] = pente->GetBinContent(binmin);
3499 y[2] = pente->GetBinContent(binmin+1);
3500 y[3] = pente->GetBinContent(binmin+2);
3501 //Calcul the polynome de Lagrange
3502 c = CalculPolynomeLagrange3(x,y);
3504 if((pente->GetBinContent(binmin+2) <= pente->GetBinContent(binmin+1))) case2 = kTRUE;
3507 if((binmin <= (nbins-3)) && (binmin == TMath::Min(binmax+4, projPH->GetNbinsX()))){
3508 min = pente->GetBinCenter(binmin);
3509 max = pente->GetBinCenter(binmin+2);
3510 x[0] = pente->GetBinCenter(binmin);
3511 x[1] = pente->GetBinCenter(binmin+1);
3512 x[2] = pente->GetBinCenter(binmin+2);
3513 y[0] = pente->GetBinContent(binmin);
3514 y[1] = pente->GetBinContent(binmin+1);
3515 y[2] = pente->GetBinContent(binmin+2);
3516 //Calcul the polynome de Lagrange
3517 c = CalculPolynomeLagrange2(x,y);
3519 if((pente->GetBinContent(binmin+2) <= pente->GetBinContent(binmin+1))) put = kFALSE;
3522 if(((binmin == (nbins-2)) && ((binmin-1) == TMath::Min(binmax+4, projPH->GetNbinsX()))) ||
3524 min = pente->GetBinCenter(binmin-1);
3525 max = pente->GetBinCenter(binmin+1);
3526 x[0] = pente->GetBinCenter(binmin-1);
3527 x[1] = pente->GetBinCenter(binmin);
3528 x[2] = pente->GetBinCenter(binmin+1);
3529 y[0] = pente->GetBinContent(binmin-1);
3530 y[1] = pente->GetBinContent(binmin);
3531 y[2] = pente->GetBinContent(binmin+1);
3532 //Calcul the polynome de Lagrange
3533 c = CalculPolynomeLagrange2(x,y);
3534 //richtung +: nothing
3535 //richtung -: nothing
3537 //case binmin = nbins-1
3539 if((binmin == (nbins-1)) && ((binmin-2) >= TMath::Min(binmax+4, projPH->GetNbinsX()))){
3540 min = pente->GetBinCenter(binmin-2);
3541 max = pente->GetBinCenter(binmin);
3542 x[0] = pente->GetBinCenter(binmin-2);
3543 x[1] = pente->GetBinCenter(binmin-1);
3544 x[2] = pente->GetBinCenter(binmin);
3545 y[0] = pente->GetBinContent(binmin-2);
3546 y[1] = pente->GetBinContent(binmin-1);
3547 y[2] = pente->GetBinContent(binmin);
3548 //Calcul the polynome de Lagrange
3549 c = CalculPolynomeLagrange2(x,y);
3550 //AliInfo("At the limit for the drift!");
3551 //fluctuation too big!
3552 //richtung +: nothing
3554 if((pente->GetBinContent(binmin-2) <= pente->GetBinContent(binmin-1))) put = kFALSE;
3556 if((binmin == (nbins-1)) && ((binmin-2) < TMath::Min(binmax+4, projPH->GetNbinsX()))) {
3558 AliInfo("At the limit for the drift and not usable!");
3562 if((binmin == (nbins-2)) && ((binmin-1) < TMath::Min(binmax+4, projPH->GetNbinsX()))){
3564 AliInfo("For the drift...problem!");
3566 //pass but should not happen
3567 if((binmin <= (nbins-3)) && (binmin < TMath::Min(binmax+4, projPH->GetNbinsX()))){
3569 AliInfo("For the drift...problem!");
3573 polynome = new TF1("pol","[0]+[1]*x+[2]*x*x+[3]*x*x*x+[4]*x*x*x*x",min,max);
3574 polynome->SetParameters(c[0],c[1],c[2],c[3],c[4]);
3575 //AliInfo(Form("GetMinimum of the function %f",polynome->GetMinimumX()));
3576 Double_t step = (max-min)/1000;
3578 Double_t minvalue = 0.0;
3579 Double_t placeminimum = min;
3580 for(Int_t o = 0; o < 1000; o++){
3581 if(o == 0) minvalue = polynome->Eval(l);
3582 if(minvalue > (polynome->Eval(l))){
3583 minvalue = polynome->Eval(l);
3588 fPhd[2] = placeminimum;
3591 Float_t fPhdt0 = 0.0;
3592 if(fTakeTheMaxPH) fPhdt0 = fPhd[1];
3593 else fPhdt0 = fPhd[0];
3595 if ((fPhd[2] > fPhd[0]) &&
3596 (fPhd[2] > fPhd[1]) &&
3597 (fPhd[1] > fPhd[0]) &&
3599 fCurrentCoef[0] = (kDrWidth) / (fPhd[2]-fPhd[1]);
3600 fNumberFitSuccess++;
3601 if (fPhdt0 >= 0.0) {
3602 fCurrentCoef2[0] = (fPhdt0 - fT0Shift) / widbins;
3603 if (fCurrentCoef2[0] < -1.0) {
3604 fCurrentCoef2[0] = fCurrentCoef2[1];
3608 fCurrentCoef2[0] = fCurrentCoef2[1];
3612 fCurrentCoef[0] = -TMath::Abs(fCurrentCoef[1]);
3613 fCurrentCoef2[0] = fCurrentCoef2[1];
3614 //printf("Fit failed!\n");
3617 if (fDebugLevel == 1) {
3618 TCanvas *cpentei = new TCanvas("cpentei","cpentei",50,50,600,800);
3621 line->SetLineColor(2);
3622 line->DrawLine(fPhd[0],0,fPhd[0],projPH->GetMaximum());
3623 line->DrawLine(fPhd[1],0,fPhd[1],projPH->GetMaximum());
3624 line->DrawLine(fPhd[2],0,fPhd[2],projPH->GetMaximum());
3625 AliInfo(Form("fPhd[0] (beginning of the signal): %f" ,(Float_t) fPhd[0]));
3626 AliInfo(Form("fPhd[1] (end of the amplification region): %f" ,(Float_t) fPhd[1]));
3627 AliInfo(Form("fPhd[2] (end of the drift region): %f" ,(Float_t) fPhd[2]));
3628 AliInfo(Form("fVriftCoef[3] (with only the drift region(default)): %f",(Float_t) fCurrentCoef[0]));
3629 TCanvas *cpentei2 = new TCanvas("cpentei2","cpentei2",50,50,600,800);
3632 TCanvas *cpentei3 = new TCanvas("cpentei3","cpentei3",50,50,600,800);
3644 projPH->SetDirectory(0);
3648 //_____________________________________________________________________________
3649 void AliTRDCalibraFit::FitPH(TH1* projPH, Int_t idect)
3652 // Fit methode for the drift velocity
3656 const Float_t kDrWidth = AliTRDgeometry::DrThick();
3659 TAxis *xpph = projPH->GetXaxis();
3660 Double_t upedge = xpph->GetBinUpEdge(xpph->GetNbins());
3662 TF1 *fPH = new TF1("fPH",AliTRDCalibraFit::PH,-0.05,3.2,6);
3663 fPH->SetParameter(0,0.469); // Scaling
3664 fPH->SetParameter(1,0.18); // Start
3665 fPH->SetParameter(2,0.0857325); // AR
3666 fPH->SetParameter(3,1.89); // DR
3667 fPH->SetParameter(4,0.08); // QA/QD
3668 fPH->SetParameter(5,0.0); // Baseline
3670 TLine *line = new TLine();
3672 fCurrentCoef[0] = 0.0;
3673 fCurrentCoef2[0] = 0.0;
3674 fCurrentCoefE = 0.0;
3675 fCurrentCoefE2 = 0.0;
3677 if (idect%fFitPHPeriode == 0) {
3679 AliInfo(Form("The detector %d will be fitted",idect));
3680 fPH->SetParameter(0,(projPH->Integral()-(projPH->GetBinContent(1)*projPH->GetNbinsX())) * 0.00028); // Scaling
3681 fPH->SetParameter(1,fPhd[0] - 0.1); // Start
3682 fPH->SetParameter(2,fPhd[1] - fPhd[0]); // AR
3683 fPH->SetParameter(3,fPhd[2] - fPhd[1]); // DR
3684 fPH->SetParameter(4,0.225); // QA/QD
3685 fPH->SetParameter(5,(Float_t) projPH->GetBinContent(1));
3687 if (fDebugLevel != 1) {
3688 projPH->Fit(fPH,"0M","",0.0,upedge);
3691 TCanvas *cpente = new TCanvas("cpente","cpente",50,50,600,800);
3693 projPH->Fit(fPH,"M+","",0.0,upedge);
3695 line->SetLineColor(4);
3696 line->DrawLine(fPH->GetParameter(1)
3698 ,fPH->GetParameter(1)
3699 ,projPH->GetMaximum());
3700 line->DrawLine(fPH->GetParameter(1)+fPH->GetParameter(2)
3702 ,fPH->GetParameter(1)+fPH->GetParameter(2)
3703 ,projPH->GetMaximum());
3704 line->DrawLine(fPH->GetParameter(1)+fPH->GetParameter(2)+fPH->GetParameter(3)
3706 ,fPH->GetParameter(1)+fPH->GetParameter(2)+fPH->GetParameter(3)
3707 ,projPH->GetMaximum());
3710 if (fPH->GetParameter(3) != 0) {
3711 fNumberFitSuccess++;
3712 fCurrentCoef[0] = kDrWidth / (fPH->GetParameter(3));
3713 fCurrentCoefE = (fPH->GetParError(3)/fPH->GetParameter(3))*fCurrentCoef[0];
3714 fCurrentCoef2[0] = fPH->GetParameter(1);
3715 fCurrentCoefE2 = fPH->GetParError(1);
3718 fCurrentCoef[0] = -TMath::Abs(fCurrentCoef[1]);
3719 fCurrentCoef2[0] = fCurrentCoef2[1];
3725 // Put the default value
3726 fCurrentCoef[0] = -TMath::Abs(fCurrentCoef[1]);
3727 fCurrentCoef2[0] = fCurrentCoef2[1];
3730 if (fDebugLevel != 1) {
3735 //_____________________________________________________________________________
3736 Bool_t AliTRDCalibraFit::FitPRFGausMI(Double_t *arraye, Double_t *arraym, Double_t *arrayme, Int_t nBins, Float_t xMin, Float_t xMax)
3739 // Fit methode for the sigma of the pad response function
3744 fCurrentCoef[0] = 0.0;
3745 fCurrentCoefE = 0.0;
3747 Double_t ret = FitGausMI(arraye, arraym, arrayme, nBins, xMin, xMax,¶m);
3750 fCurrentCoef[0] = -fCurrentCoef[1];
3754 fNumberFitSuccess++;
3755 fCurrentCoef[0] = param[2];
3756 fCurrentCoefE = ret;
3760 //_____________________________________________________________________________
3761 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 kError)
3764 // Fit methode for the sigma of the pad response function
3767 //We should have at least 3 points
3768 if(nBins <=3) return -4.0;
3770 TLinearFitter fitter(3,"pol2");
3771 fitter.StoreData(kFALSE);
3772 fitter.ClearPoints();
3774 Float_t binWidth = (xMax-xMin)/(Float_t)nBins;
3775 Float_t entries = 0;
3776 Int_t nbbinwithentries = 0;
3777 for (Int_t i=0; i<nBins; i++){
3779 if(arraye[i] > 15) nbbinwithentries++;
3780 //printf("entries for i %d: %f\n",i,arraye[i]);
3782 if ((entries<700) || (nbbinwithentries < ((Int_t)(nBins/2)))) return -4;
3783 //printf("entries %f\n",entries);
3784 //printf("nbbinwithentries %d\n",nbbinwithentries);
3787 Float_t errorm = 0.0;
3788 Float_t errorn = 0.0;
3789 Float_t error = 0.0;
3792 for (Int_t ibin=0;ibin<nBins; ibin++){
3793 Float_t entriesI = arraye[ibin];
3794 Float_t valueI = arraym[ibin];
3795 Double_t xcenter = 0.0;
3797 if ((entriesI>15) && (valueI>0.0)){
3798 xcenter = xMin+(ibin+0.5)*binWidth;
3803 if((valueI + 0.01) > 0.0) errorm = TMath::Log((valueI + 0.01)/valueI);
3804 if((valueI - 0.01) > 0.0) errorn = TMath::Log((valueI - 0.01)/valueI);
3805 error = TMath::Max(TMath::Abs(errorm),TMath::Abs(errorn));
3808 if((valueI + arrayme[ibin]) > 0.0) errorm = TMath::Log((valueI + arrayme[ibin])/valueI);
3809 if((valueI - arrayme[ibin]) > 0.0) errorn = TMath::Log((valueI - arrayme[ibin])/valueI);
3810 error = TMath::Max(TMath::Abs(errorm),TMath::Abs(errorn));
3812 if(error == 0.0) continue;
3813 val = TMath::Log(Float_t(valueI));
3814 fitter.AddPoint(&xcenter,val,error);
3818 if(fDebugLevel > 1){
3820 if ( !fDebugStreamer ) {
3822 TDirectory *backup = gDirectory;
3823 fDebugStreamer = new TTreeSRedirector("TRDDebugFit.root");
3824 if ( backup ) backup->cd(); //we don't want to be cd'd to the debug streamer
3827 Int_t detector = fCountDet;
3828 Int_t plane = GetPlane(fCountDet);
3831 (* fDebugStreamer) << "FitGausMIFill"<<
3832 "detector="<<detector<<
3836 "entriesI="<<entriesI<<
3839 "xcenter="<<xcenter<<
3849 if(npoints <=3) return -4.0;
3854 fitter.GetParameters(par);
3855 chi2 = fitter.GetChisquare()/Float_t(npoints);
3858 if (!param) param = new TVectorD(3);
3859 if(par[2] == 0.0) return -4.0;
3860 Double_t x = TMath::Sqrt(TMath::Abs(-2*par[2]));
3861 Double_t deltax = (fitter.GetParError(2))/x;
3862 Double_t errorparam2 = TMath::Abs(deltax)/(x*x);
3865 (*param)[1] = par[1]/(-2.*par[2]);
3866 (*param)[2] = 1./TMath::Sqrt(TMath::Abs(-2.*par[2]));
3867 Double_t lnparam0 = par[0]+ par[1]* (*param)[1] + par[2]*(*param)[1]*(*param)[1];
3868 if ( lnparam0>307 ) return -4;
3869 (*param)[0] = TMath::Exp(lnparam0);
3871 if(fDebugLevel > 1){
3873 if ( !fDebugStreamer ) {
3875 TDirectory *backup = gDirectory;
3876 fDebugStreamer = new TTreeSRedirector("TRDDebugFit.root");
3877 if ( backup ) backup->cd(); //we don't want to be cd'd to the debug streamer
3880 Int_t detector = fCountDet;
3881 Int_t plane = GetPlane(fCountDet);
3884 (* fDebugStreamer) << "FitGausMIFit"<<
3885 "detector="<<detector<<
3888 "errorsigma="<<chi2<<
3889 "mean="<<(*param)[1]<<
3890 "sigma="<<(*param)[2]<<
3891 "constant="<<(*param)[0]<<
3896 if((chi2/(*param)[2]) > 0.1){
3898 chi2 = FitGausMI(arraye,arraym,arrayme,nBins,xMin,xMax,param,kFALSE);
3903 if(fDebugLevel == 1){
3904 TString name("PRF");
3905 name += (Int_t)xMin;
3906 name += (Int_t)xMax;
3907 TCanvas *c1 = new TCanvas((const char *)name,(const char *)name,50,50,600,800);
3910 TH1F *histo = new TH1F((const char *)name,(const char *)name,nBins,xMin,xMax);
3911 for(Int_t k = 0; k < nBins; k++){
3912 histo->SetBinContent(k+1,arraym[k]);
3913 histo->SetBinError(k+1,arrayme[k]);
3916 name += "functionf";
3917 TF1 *f1= new TF1((const char*)name,"[0]*exp(-(x-[1])^2/(2*[2]*[2]))",xMin,xMax);
3918 f1->SetParameter(0, (*param)[0]);
3919 f1->SetParameter(1, (*param)[1]);
3920 f1->SetParameter(2, (*param)[2]);
3928 //_____________________________________________________________________________
3929 void AliTRDCalibraFit::FitPRF(TH1 *projPRF)
3932 // Fit methode for the sigma of the pad response function
3935 fCurrentCoef[0] = 0.0;
3936 fCurrentCoefE = 0.0;
3938 if (fDebugLevel != 1) {
3939 projPRF->Fit("gaus","0M","",-fRangeFitPRF,fRangeFitPRF);
3942 TCanvas *cfit = new TCanvas("cfit","cfit",50,50,600,800);
3944 projPRF->Fit("gaus","M+","",-fRangeFitPRF,fRangeFitPRF);
3947 fCurrentCoef[0] = projPRF->GetFunction("gaus")->GetParameter(2);
3948 fCurrentCoefE = projPRF->GetFunction("gaus")->GetParError(2);
3949 if(fCurrentCoef[0] <= 0.0) fCurrentCoef[0] = -fCurrentCoef[1];
3951 fNumberFitSuccess++;
3954 //_____________________________________________________________________________
3955 void AliTRDCalibraFit::RmsPRF(TH1 *projPRF)
3958 // Fit methode for the sigma of the pad response function
3960 fCurrentCoef[0] = 0.0;
3961 fCurrentCoefE = 0.0;
3962 if (fDebugLevel == 1) {
3963 TCanvas *cfit = new TCanvas("cfit","cfit",50,50,600,800);
3967 fCurrentCoef[0] = projPRF->GetRMS();
3968 if(fCurrentCoef[0] <= 0.0) fCurrentCoef[0] = -fCurrentCoef[1];
3970 fNumberFitSuccess++;
3973 //_____________________________________________________________________________
3974 void AliTRDCalibraFit::FitTnpRange(Double_t *arraye, Double_t *arraym, Double_t *arrayme, Int_t nbg, Int_t nybins)
3977 // Fit methode for the sigma of the pad response function with 2*nbg tan bins
3980 TLinearFitter linearfitter = TLinearFitter(3,"pol2");
3983 Int_t nbins = (Int_t)(nybins/(2*nbg));
3984 Float_t lowedge = -3.0*nbg;
3985 Float_t upedge = lowedge + 3.0;
3988 Double_t xvalues = -0.2*nbg+0.1;
3990 Int_t total = 2*nbg;
3993 for(Int_t k = 0; k < total; k++){
3994 if(FitPRFGausMI(arraye+offset, arraym+offset, arrayme+offset, nbins, lowedge, upedge)){
3996 y = fCurrentCoef[0]*fCurrentCoef[0];
3997 linearfitter.AddPoint(&xvalues,y,2*fCurrentCoefE*fCurrentCoef[0]);
4000 if(fDebugLevel > 1){
4002 if ( !fDebugStreamer ) {
4004 TDirectory *backup = gDirectory;
4005 fDebugStreamer = new TTreeSRedirector("TRDDebugFit.root");
4006 if ( backup ) backup->cd(); //we don't want to be cd'd to the debug streamer
4009 Int_t detector = fCountDet;
4010 Int_t plane = GetPlane(fCountDet);
4011 Int_t nbtotal = total;
4013 Float_t low = lowedge;
4014 Float_t up = upedge;
4015 Float_t tnp = xvalues;
4016 Float_t wid = fCurrentCoef[0];
4017 Float_t widfE = fCurrentCoefE;
4019 (* fDebugStreamer) << "FitTnpRangeFill"<<
4020 "detector="<<detector<<
4022 "nbtotal="<<nbtotal<<
4040 fCurrentCoefE = 0.0;
4041 fCurrentCoef[0] = 0.0;
4043 //printf("npoints\n",npoints);
4046 fCurrentCoef[0] = -TMath::Abs(fCurrentCoef[1]);
4051 linearfitter.Eval();
4052 linearfitter.GetParameters(pars0);
4053 Double_t pointError0 = TMath::Sqrt(linearfitter.GetChisquare()/npoints);
4054 Double_t errorsx0 = linearfitter.GetParError(2)*pointError0;
4055 Double_t min0 = 0.0;
4056 Double_t ermin0 = 0.0;
4057 //Double_t prfe0 = 0.0;
4058 Double_t prf0 = 0.0;
4059 if((pars0[2] > 0.0) && (pars0[1] != 0.0)) {
4060 min0 = -pars0[1]/(2*pars0[2]);
4061 ermin0 = TMath::Abs(min0*(errorsx0/pars0[2]+linearfitter.GetParError(1)*pointError0/pars0[1]));
4062 prf0 = pars0[0]+pars0[1]*min0+pars0[2]*min0*min0;
4065 prfe0 = linearfitter->GetParError(0)*pointError0
4066 +(linearfitter->GetParError(1)*pointError0/pars0[1]+ermin0/min0)*pars0[1]*min0
4067 +(linearfitter->GetParError(2)*pointError0/pars0[2]+2*ermin0/min0)*pars0[2]*min0*min0;
4068 prfe0 = prfe0/(2*TMath::Sqrt(prf0));
4069 fCurrentCoefE = (Float_t) prfe0;
4071 fCurrentCoef[0] = (Float_t) TMath::Sqrt(TMath::Abs(prf0));
4074 fCurrentCoef[0] = -TMath::Abs(fCurrentCoef[1]);
4078 fCurrentCoef[0] = -TMath::Abs(fCurrentCoef[1]);
4081 if(fDebugLevel > 1){
4083 if ( !fDebugStreamer ) {
4085 TDirectory *backup = gDirectory;
4086 fDebugStreamer = new TTreeSRedirector("TRDDebugFit.root");
4087 if ( backup ) backup->cd(); //we don't want to be cd'd to the debug streamer
4090 Int_t detector = fCountDet;
4091 Int_t plane = GetPlane(fCountDet);
4092 Int_t nbtotal = total;
4093 Double_t colsize[6] = {0.635,0.665,0.695,0.725,0.755,0.785};
4094 Double_t sigmax = TMath::Sqrt(TMath::Abs(pars0[2]))*10000*colsize[plane];
4096 (* fDebugStreamer) << "FitTnpRangeFit"<<
4097 "detector="<<detector<<
4099 "nbtotal="<<nbtotal<<
4103 "npoints="<<npoints<<
4106 "sigmaprf="<<fCurrentCoef[0]<<
4107 "sigprf="<<fCurrentCoef[1]<<
4114 //_____________________________________________________________________________
4115 void AliTRDCalibraFit::FitMean(TH1 *projch, Double_t nentries, Double_t mean)
4118 // Only mean methode for the gain factor
4121 fCurrentCoef[0] = mean;
4122 fCurrentCoefE = 0.0;
4123 if(nentries > 0) fCurrentCoefE = projch->GetRMS()/TMath::Sqrt(nentries);
4124 if (fDebugLevel == 1) {
4125 TCanvas *cpmean = new TCanvas("cpmean","cpmean",50,50,600,800);
4129 CalculChargeCoefMean(kTRUE);
4130 fNumberFitSuccess++;
4132 //_____________________________________________________________________________
4133 void AliTRDCalibraFit::FitMeanW(TH1 *projch, Double_t nentries)
4136 // mean w methode for the gain factor
4140 Int_t nybins = projch->GetNbinsX();
4142 //The weight function
4143 Double_t a = 0.00228515;
4144 Double_t b = -0.00231487;
4145 Double_t c = 0.00044298;
4146 Double_t d = -0.00379239;
4147 Double_t e = 0.00338349;
4157 //A arbitrary error for the moment
4158 fCurrentCoefE = 0.0;
4159 fCurrentCoef[0] = 0.0;
4162 Double_t sumw = 0.0;
4164 Float_t sumAll = (Float_t) nentries;
4165 Int_t sumCurrent = 0;
4166 for(Int_t k = 0; k <nybins; k++){
4167 Double_t fraction = Float_t(sumCurrent)/Float_t(sumAll);
4168 if (fraction>0.95) break;
4169 Double_t weight = a + b*fraction + c*fraction*fraction + d *fraction*fraction*fraction+
4170 e*fraction*fraction*fraction*fraction;
4171 sumw += weight*projch->GetBinContent(k+1)*projch->GetBinCenter(k+1);
4172 sum += weight*projch->GetBinContent(k+1);
4173 sumCurrent += (Int_t) projch->GetBinContent(k+1);
4174 //printf("fraction %f, weight %f, bincontent %f\n",fraction,weight,projch->GetBinContent(k+1));
4176 if(sum > 0.0) fCurrentCoef[0] = (sumw/sum);
4178 if (fDebugLevel == 1) {
4179 TCanvas *cpmeanw = new TCanvas("cpmeanw","cpmeanw",50,50,600,800);
4183 fNumberFitSuccess++;
4184 CalculChargeCoefMean(kTRUE);
4186 //_____________________________________________________________________________
4187 void AliTRDCalibraFit::FitMeanWSm(TH1 *projch, Float_t sumAll)
4190 // mean w methode for the gain factor
4194 Int_t nybins = projch->GetNbinsX();
4196 //The weight function
4197 Double_t a = 0.00228515;
4198 Double_t b = -0.00231487;
4199 Double_t c = 0.00044298;
4200 Double_t d = -0.00379239;
4201 Double_t e = 0.00338349;
4211 //A arbitrary error for the moment
4212 fCurrentCoefE = 0.0;
4213 fCurrentCoef[0] = 0.0;
4216 Double_t sumw = 0.0;
4218 Int_t sumCurrent = 0;
4219 for(Int_t k = 0; k <nybins; k++){
4220 Double_t fraction = Float_t(sumCurrent)/Float_t(sumAll);
4221 if (fraction>0.95) break;
4222 Double_t weight = a + b*fraction + c*fraction*fraction + d *fraction*fraction*fraction+
4223 e*fraction*fraction*fraction*fraction;
4224 sumw += weight*projch->GetBinContent(k+1)*projch->GetBinCenter(k+1);
4225 sum += weight*projch->GetBinContent(k+1);
4226 sumCurrent += (Int_t) projch->GetBinContent(k+1);
4227 //printf("fraction %f, weight %f, bincontent %f\n",fraction,weight,projch->GetBinContent(k+1));
4229 if(sum > 0.0) fCurrentCoef[0] = (sumw/sum);
4231 if (fDebugLevel == 1) {
4232 TCanvas *cpmeanw = new TCanvas("cpmeanw","cpmeanw",50,50,600,800);
4236 fNumberFitSuccess++;
4238 //_____________________________________________________________________________
4239 void AliTRDCalibraFit::FitCH(TH1 *projch, Double_t mean)
4242 // Fit methode for the gain factor
4245 fCurrentCoef[0] = 0.0;
4246 fCurrentCoefE = 0.0;
4247 Double_t chisqrl = 0.0;
4248 Double_t chisqrg = 0.0;
4249 Double_t chisqr = 0.0;
4250 TF1 *fLandauGaus = new TF1("fLandauGaus",FuncLandauGaus,0,300,5);
4252 projch->Fit("landau","0",""
4253 ,(Double_t) mean/fBeginFitCharge
4254 ,projch->GetBinCenter(projch->GetNbinsX()));
4255 Double_t l3P0 = projch->GetFunction("landau")->GetParameter(0);
4256 Double_t l3P1 = projch->GetFunction("landau")->GetParameter(1);
4257 Double_t l3P2 = projch->GetFunction("landau")->GetParameter(2);
4258 chisqrl = projch->GetFunction("landau")->GetChisquare();
4260 projch->Fit("gaus","0",""
4261 ,(Double_t) mean/fBeginFitCharge
4262 ,projch->GetBinCenter(projch->GetNbinsX()));
4263 Double_t g3P0 = projch->GetFunction("gaus")->GetParameter(0);
4264 Double_t g3P2 = projch->GetFunction("gaus")->GetParameter(2);
4265 chisqrg = projch->GetFunction("gaus")->GetChisquare();
4267 fLandauGaus->SetParameters(l3P0,l3P1,l3P2,g3P0,g3P2);
4268 if (fDebugLevel != 1) {
4269 projch->Fit("fLandauGaus","0",""
4270 ,(Double_t) mean/fBeginFitCharge
4271 ,projch->GetBinCenter(projch->GetNbinsX()));
4272 chisqr = projch->GetFunction("fLandauGaus")->GetChisquare();
4275 TCanvas *cp = new TCanvas("cp","cp",50,50,600,800);
4277 projch->Fit("fLandauGaus","+",""
4278 ,(Double_t) mean/fBeginFitCharge
4279 ,projch->GetBinCenter(projch->GetNbinsX()));
4280 chisqr = projch->GetFunction("fLandauGaus")->GetChisquare();
4282 fLandauGaus->Draw("same");
4285 if ((projch->GetFunction("fLandauGaus")->GetParameter(1) > 0) && (projch->GetFunction("fLandauGaus")->GetParError(1) < (0.05*projch->GetFunction("fLandauGaus")->GetParameter(1))) && (chisqr < chisqrl) && (chisqr < chisqrg)) {
4286 //if ((projch->GetFunction("fLandauGaus")->GetParameter(1) > 0) && (chisqr < chisqrl) && (chisqr < chisqrg)) {
4287 fNumberFitSuccess++;
4288 CalculChargeCoefMean(kTRUE);
4289 fCurrentCoef[0] = projch->GetFunction("fLandauGaus")->GetParameter(1);
4290 fCurrentCoefE = projch->GetFunction("fLandauGaus")->GetParError(1);
4293 CalculChargeCoefMean(kFALSE);
4294 fCurrentCoef[0] = -TMath::Abs(fCurrentCoef[1]);
4297 if (fDebugLevel != 1) {
4302 //_____________________________________________________________________________
4303 void AliTRDCalibraFit::FitBisCH(TH1* projch, Double_t mean)
4306 // Fit methode for the gain factor more time consuming
4310 //Some parameters to initialise
4311 Double_t widthLandau, widthGaus, MPV, Integral;
4312 Double_t chisquarel = 0.0;
4313 Double_t chisquareg = 0.0;
4314 projch->Fit("landau","0M+",""
4316 ,projch->GetBinCenter(projch->GetNbinsX()));
4317 widthLandau = projch->GetFunction("landau")->GetParameter(2);
4318 chisquarel = projch->GetFunction("landau")->GetChisquare();
4319 projch->Fit("gaus","0M+",""
4321 ,projch->GetBinCenter(projch->GetNbinsX()));
4322 widthGaus = projch->GetFunction("gaus")->GetParameter(2);
4323 chisquareg = projch->GetFunction("gaus")->GetChisquare();
4325 MPV = (projch->GetFunction("landau")->GetParameter(1))/2;
4326 Integral = (projch->GetFunction("gaus")->Integral(0.3*mean,3*mean)+projch->GetFunction("landau")->Integral(0.3*mean,3*mean))/2;
4328 // Setting fit range and start values
4330 //Double_t sv[4] = { l3P2, fChargeCoef[1], projch->Integral("width"), fG3P2 };
4331 //Double_t sv[4] = { fL3P2, fChargeCoef[1], fL3P0, fG3P2 };
4332 Double_t sv[4] = { widthLandau, MPV, Integral, widthGaus};
4333 Double_t pllo[4] = { 0.001, 0.001, projch->Integral()/3, 0.001};
4334 Double_t plhi[4] = { 300.0, 300.0, 30*projch->Integral(), 300.0};
4335 Double_t fp[4] = { 1.0, 1.0, 1.0, 1.0 };
4336 Double_t fpe[4] = { 1.0, 1.0, 1.0, 1.0 };
4339 fCurrentCoef[0] = 0.0;
4340 fCurrentCoefE = 0.0;
4344 TF1 *fitsnr = LanGauFit(projch,&fr[0],&sv[0]
4349 Double_t projchPeak;
4350 Double_t projchFWHM;
4351 LanGauPro(fp,projchPeak,projchFWHM);
4353 if ((fp[1] > 0) && ((fpe[1] < (0.05*fp[1])) && (chisqr < chisquarel) && (chisqr < chisquareg))) {
4354 //if ((fp[1] > 0) && ((chisqr < chisquarel) && (chisqr < chisquareg))) {
4355 fNumberFitSuccess++;
4356 CalculChargeCoefMean(kTRUE);
4357 fCurrentCoef[0] = fp[1];
4358 fCurrentCoefE = fpe[1];
4359 //chargeCoefE2 = chisqr;
4362 CalculChargeCoefMean(kFALSE);
4363 fCurrentCoef[0] = -TMath::Abs(fCurrentCoef[1]);
4365 if (fDebugLevel == 1) {
4366 AliInfo(Form("fChargeCoef[0]: %f",(Float_t) fCurrentCoef[0]));
4367 TCanvas *cpy = new TCanvas("cpy","cpy",50,50,600,800);
4370 fitsnr->Draw("same");
4376 //_____________________________________________________________________________
4377 Double_t *AliTRDCalibraFit::CalculPolynomeLagrange2(Double_t *x, Double_t *y)
4380 // Calcul the coefficients of the polynome passant par ces trois points de degre 2
4382 Double_t *c = new Double_t[5];
4383 Double_t x0 = y[0]/((x[0]-x[1])*(x[0]-x[2]));
4384 Double_t x1 = y[1]/((x[1]-x[0])*(x[1]-x[2]));
4385 Double_t x2 = y[2]/((x[2]-x[0])*(x[2]-x[1]));
4390 c[1] = -(x0*(x[1]+x[2])+x1*(x[0]+x[2])+x2*(x[0]+x[1]));
4391 c[0] = x0*x[1]*x[2]+x1*x[0]*x[2]+x2*x[0]*x[1];
4398 //_____________________________________________________________________________
4399 Double_t *AliTRDCalibraFit::CalculPolynomeLagrange3(Double_t *x, Double_t *y)
4402 // Calcul the coefficients of the polynome passant par ces quatre points de degre 3
4404 Double_t *c = new Double_t[5];
4405 Double_t x0 = y[0]/((x[0]-x[1])*(x[0]-x[2])*(x[0]-x[3]));
4406 Double_t x1 = y[1]/((x[1]-x[0])*(x[1]-x[2])*(x[1]-x[3]));
4407 Double_t x2 = y[2]/((x[2]-x[0])*(x[2]-x[1])*(x[2]-x[3]));
4408 Double_t x3 = y[3]/((x[3]-x[0])*(x[3]-x[1])*(x[3]-x[2]));
4412 c[2] = -(x0*(x[1]+x[2]+x[3])
4413 +x1*(x[0]+x[2]+x[3])
4414 +x2*(x[0]+x[1]+x[3])
4415 +x3*(x[0]+x[1]+x[2]));
4416 c[1] = (x0*(x[1]*x[2]+x[1]*x[3]+x[2]*x[3])
4417 +x1*(x[0]*x[2]+x[0]*x[3]+x[2]*x[3])
4418 +x2*(x[0]*x[1]+x[0]*x[3]+x[1]*x[3])
4419 +x3*(x[0]*x[1]+x[0]*x[2]+x[1]*x[2]));
4421 c[0] = -(x0*x[1]*x[2]*x[3]
4424 +x3*x[0]*x[1]*x[2]);
4432 //_____________________________________________________________________________
4433 Double_t *AliTRDCalibraFit::CalculPolynomeLagrange4(Double_t *x, Double_t *y)
4436 // Calcul the coefficients of the polynome passant par ces cinqs points de degre 4
4438 Double_t *c = new Double_t[5];
4439 Double_t x0 = y[0]/((x[0]-x[1])*(x[0]-x[2])*(x[0]-x[3])*(x[0]-x[4]));
4440 Double_t x1 = y[1]/((x[1]-x[0])*(x[1]-x[2])*(x[1]-x[3])*(x[1]-x[4]));
4441 Double_t x2 = y[2]/((x[2]-x[0])*(x[2]-x[1])*(x[2]-x[3])*(x[2]-x[4]));
4442 Double_t x3 = y[3]/((x[3]-x[0])*(x[3]-x[1])*(x[3]-x[2])*(x[3]-x[4]));
4443 Double_t x4 = y[4]/((x[4]-x[0])*(x[4]-x[1])*(x[4]-x[2])*(x[4]-x[3]));
4446 c[4] = x0+x1+x2+x3+x4;
4447 c[3] = -(x0*(x[1]+x[2]+x[3]+x[4])
4448 +x1*(x[0]+x[2]+x[3]+x[4])
4449 +x2*(x[0]+x[1]+x[3]+x[4])
4450 +x3*(x[0]+x[1]+x[2]+x[4])
4451 +x4*(x[0]+x[1]+x[2]+x[3]));
4452 c[2] = (x0*(x[1]*x[2]+x[1]*x[3]+x[1]*x[4]+x[2]*x[3]+x[2]*x[4]+x[3]*x[4])
4453 +x1*(x[0]*x[2]+x[0]*x[3]+x[0]*x[4]+x[2]*x[3]+x[2]*x[4]+x[3]*x[4])
4454 +x2*(x[0]*x[1]+x[0]*x[3]+x[0]*x[4]+x[1]*x[3]+x[1]*x[4]+x[3]*x[4])
4455 +x3*(x[0]*x[1]+x[0]*x[2]+x[0]*x[4]+x[1]*x[2]+x[1]*x[4]+x[2]*x[4])
4456 +x4*(x[0]*x[1]+x[0]*x[2]+x[0]*x[3]+x[1]*x[2]+x[1]*x[3]+x[2]*x[3]));
4458 c[1] = -(x0*(x[1]*x[2]*x[3]+x[1]*x[2]*x[4]+x[1]*x[3]*x[4]+x[2]*x[3]*x[4])
4459 +x1*(x[0]*x[2]*x[3]+x[0]*x[2]*x[4]+x[0]*x[3]*x[4]+x[2]*x[3]*x[4])
4460 +x2*(x[0]*x[1]*x[3]+x[0]*x[1]*x[4]+x[0]*x[3]*x[4]+x[1]*x[3]*x[4])
4461 +x3*(x[0]*x[1]*x[2]+x[0]*x[1]*x[4]+x[0]*x[2]*x[4]+x[1]*x[2]*x[4])
4462 +x4*(x[0]*x[1]*x[2]+x[0]*x[1]*x[3]+x[0]*x[2]*x[3]+x[1]*x[2]*x[3]));
4464 c[0] = (x0*x[1]*x[2]*x[3]*x[4]
4465 +x1*x[0]*x[2]*x[3]*x[4]
4466 +x2*x[0]*x[1]*x[3]*x[4]
4467 +x3*x[0]*x[1]*x[2]*x[4]
4468 +x4*x[0]*x[1]*x[2]*x[3]);
4474 //_____________________________________________________________________________
4475 void AliTRDCalibraFit::NormierungCharge()
4478 // Normalisation of the gain factor resulting for the fits
4481 // Calcul of the mean of choosen method by fFitChargeNDB
4483 //printf("total number of entries %d\n",fVectorFitCH->GetEntriesFast());
4484 for (Int_t k = 0; k < (Int_t) fVectorFit.GetEntriesFast(); k++) {
4486 Int_t detector = ((AliTRDFitInfo *) fVectorFit.At(k))->GetDetector();
4487 Float_t *coef = ((AliTRDFitInfo *) fVectorFit.At(k))->GetCoef();
4488 //printf("detector %d coef[0] %f\n",detector,coef[0]);
4489 if (GetChamber(detector) == 2) {
4492 if (GetChamber(detector) != 2) {
4495 for (Int_t j = 0; j < total; j++) {
4503 fScaleFitFactor = fScaleFitFactor / sum;
4506 fScaleFitFactor = 1.0;
4509 //methode de boeuf mais bon...
4510 Double_t scalefactor = fScaleFitFactor;
4512 if(fDebugLevel > 1){
4514 if ( !fDebugStreamer ) {
4516 TDirectory *backup = gDirectory;
4517 fDebugStreamer = new TTreeSRedirector("TRDDebugFit.root");
4518 if ( backup ) backup->cd(); //we don't want to be cd'd to the debug streamer
4520 (* fDebugStreamer) << "ScaleFactor"<<
4521 "scalefactor="<<scalefactor<<
4525 //_____________________________________________________________________________
4526 TH1I *AliTRDCalibraFit::ReBin(TH1I *hist) const
4529 // Rebin of the 1D histo for the gain calibration if needed.
4530 // you have to choose fRebin, divider of fNumberBinCharge
4533 TAxis *xhist = hist->GetXaxis();
4534 TH1I *rehist = new TH1I("projrebin","",(Int_t) xhist->GetNbins()/fRebin
4535 ,xhist->GetBinLowEdge(1)
4536 ,xhist->GetBinUpEdge(xhist->GetNbins()));
4538 AliInfo(Form("fRebin: %d",fRebin));
4540 for (Int_t k = 1; k <= (Int_t) xhist->GetNbins()/fRebin; k++) {
4542 for (Int_t ji = i; ji < i+fRebin; ji++) {
4543 sum += hist->GetBinContent(ji);
4546 rehist->SetBinContent(k,sum);
4554 //_____________________________________________________________________________
4555 TH1F *AliTRDCalibraFit::ReBin(TH1F *hist) const
4558 // Rebin of the 1D histo for the gain calibration if needed
4559 // you have to choose fRebin divider of fNumberBinCharge
4562 TAxis *xhist = hist->GetXaxis();
4563 TH1F *rehist = new TH1F("projrebin","",(Int_t) xhist->GetNbins()/fRebin
4564 ,xhist->GetBinLowEdge(1)
4565 ,xhist->GetBinUpEdge(xhist->GetNbins()));
4567 AliInfo(Form("fRebin: %d",fRebin));
4569 for (Int_t k = 1; k <= (Int_t) xhist->GetNbins()/fRebin; k++) {
4571 for (Int_t ji = i; ji < i+fRebin; ji++) {
4572 sum += hist->GetBinContent(ji);
4575 rehist->SetBinContent(k,sum);
4583 //_____________________________________________________________________________
4584 TH1F *AliTRDCalibraFit::CorrectTheError(TGraphErrors *hist)
4587 // In the case of the vectors method the trees contains TGraphErrors for PH and PRF
4588 // to be able to add them after
4589 // We convert it to a TH1F to be able to applied the same fit function method
4590 // After having called this function you can not add the statistics anymore
4595 Int_t nbins = hist->GetN();
4596 Double_t *x = hist->GetX();
4597 Double_t *entries = hist->GetEX();
4598 Double_t *mean = hist->GetY();
4599 Double_t *square = hist->GetEY();
4600 fEntriesCurrent = 0;
4606 Double_t step = x[1] - x[0];
4607 Double_t minvalue = x[0] - step/2;
4608 Double_t maxvalue = x[(nbins-1)] + step/2;
4610 rehist = new TH1F("projcorrecterror","",nbins,minvalue,maxvalue);
4612 for (Int_t k = 0; k < nbins; k++) {
4613 rehist->SetBinContent(k+1,mean[k]);
4614 if (entries[k] > 0.0) {
4615 fEntriesCurrent += (Int_t) entries[k];
4616 Double_t d = TMath::Abs(square[k] - (mean[k]*mean[k]));
4617 rehist->SetBinError(k+1,TMath::Sqrt(d/entries[k]));
4620 rehist->SetBinError(k+1,0.0);
4624 if(fEntriesCurrent > 0) fNumberEnt++;
4630 //____________Some basic geometry function_____________________________________
4633 //_____________________________________________________________________________
4634 Int_t AliTRDCalibraFit::GetPlane(Int_t d) const
4637 // Reconstruct the plane number from the detector number
4640 return ((Int_t) (d % 6));
4644 //_____________________________________________________________________________
4645 Int_t AliTRDCalibraFit::GetChamber(Int_t d) const
4648 // Reconstruct the chamber number from the detector number
4652 return ((Int_t) (d % 30) / fgkNplan);
4656 //_____________________________________________________________________________
4657 Int_t AliTRDCalibraFit::GetSector(Int_t d) const
4660 // Reconstruct the sector number from the detector number
4664 return ((Int_t) (d / fg));
4669 //____________Fill and Init tree Gain, PRF, Vdrift and T0______________________
4671 //_______________________________________________________________________________
4672 void AliTRDCalibraFit::ResetVectorFit()
4674 fVectorFit.SetOwner();
4676 fVectorFit2.SetOwner();
4677 fVectorFit2.Clear();
4681 //____________Private Functions________________________________________________
4684 //_____________________________________________________________________________
4685 Double_t AliTRDCalibraFit::PH(Double_t *x, Double_t *par)
4688 // Function for the fit
4691 //TF1 *fAsymmGauss = new TF1("fAsymmGauss",AsymmGauss,0,4,6);
4693 //PARAMETERS FOR FIT PH
4695 //fAsymmGauss->SetParameter(0,0.113755);
4696 //fAsymmGauss->SetParameter(1,0.350706);
4697 //fAsymmGauss->SetParameter(2,0.0604244);
4698 //fAsymmGauss->SetParameter(3,7.65596);
4699 //fAsymmGauss->SetParameter(4,1.00124);
4700 //fAsymmGauss->SetParameter(5,0.870597); // No tail cancelation
4708 Double_t dx = 0.005;
4709 Double_t xs = par[1];
4711 Double_t paras[2] = { 0.0, 0.0 };
4714 if ((xs >= par[1]) &&
4715 (xs < (par[1]+par[2]))) {
4716 //fAsymmGauss->SetParameter(0,par[0]);
4717 //fAsymmGauss->SetParameter(1,xs);
4718 //ss += fAsymmGauss->Eval(xx);
4721 ss += AsymmGauss(&xx,paras);
4723 if ((xs >= (par[1]+par[2])) &&
4724 (xs < (par[1]+par[2]+par[3]))) {
4725 //fAsymmGauss->SetParameter(0,par[0]*par[4]);
4726 //fAsymmGauss->SetParameter(1,xs);
4727 //ss += fAsymmGauss->Eval(xx);
4728 paras[0] = par[0]*par[4];
4730 ss += AsymmGauss(&xx,paras);
4739 //_____________________________________________________________________________
4740 Double_t AliTRDCalibraFit::AsymmGauss(Double_t *x, Double_t *par)
4743 // Function for the fit
4746 //par[0] = normalization
4754 Double_t par1save = par[1];
4755 //Double_t par2save = par[2];
4756 Double_t par2save = 0.0604244;
4757 //Double_t par3save = par[3];
4758 Double_t par3save = 7.65596;
4759 //Double_t par5save = par[5];
4760 Double_t par5save = 0.870597;
4761 Double_t dx = x[0] - par1save;
4763 Double_t sigma2 = par2save*par2save;
4764 Double_t sqrt2 = TMath::Sqrt(2.0);
4765 Double_t exp1 = par3save * TMath::Exp(-par3save * (dx - 0.5 * par3save * sigma2))
4766 * (1.0 - TMath::Erf((par3save * sigma2 - dx) / (sqrt2 * par2save)));
4767 Double_t exp2 = par5save * TMath::Exp(-par5save * (dx - 0.5 * par5save * sigma2))
4768 * (1.0 - TMath::Erf((par5save * sigma2 - dx) / (sqrt2 * par2save)));
4770 //return par[0]*(exp1+par[4]*exp2);
4771 return par[0] * (exp1 + 1.00124 * exp2);
4775 //_____________________________________________________________________________
4776 Double_t AliTRDCalibraFit::FuncLandauGaus(Double_t *x, Double_t *par)
4779 // Sum Landau + Gaus with identical mean
4782 Double_t valLandau = par[0] * TMath::Landau(x[0],par[1],par[2]);
4783 //Double_t valGaus = par[3] * TMath::Gaus(x[0],par[4],par[5]);
4784 Double_t valGaus = par[3] * TMath::Gaus(x[0],par[1],par[4]);
4785 Double_t val = valLandau + valGaus;
4791 //_____________________________________________________________________________
4792 Double_t AliTRDCalibraFit::LanGauFun(Double_t *x, Double_t *par)
4795 // Function for the fit
4798 // par[0]=Width (scale) parameter of Landau density
4799 // par[1]=Most Probable (MP, location) parameter of Landau density
4800 // par[2]=Total area (integral -inf to inf, normalization constant)
4801 // par[3]=Width (sigma) of convoluted Gaussian function
4803 // In the Landau distribution (represented by the CERNLIB approximation),
4804 // the maximum is located at x=-0.22278298 with the location parameter=0.
4805 // This shift is corrected within this function, so that the actual
4806 // maximum is identical to the MP parameter.
4809 // Numeric constants
4810 Double_t invsq2pi = 0.3989422804014; // (2 pi)^(-1/2)
4811 Double_t mpshift = -0.22278298; // Landau maximum location
4813 // Control constants
4814 Double_t np = 100.0; // Number of convolution steps
4815 Double_t sc = 5.0; // Convolution extends to +-sc Gaussian sigmas
4827 // MP shift correction
4828 mpc = par[1] - mpshift * par[0];
4830 // Range of convolution integral
4831 xlow = x[0] - sc * par[3];
4832 xupp = x[0] + sc * par[3];
4834 step = (xupp - xlow) / np;
4836 // Convolution integral of Landau and Gaussian by sum
4837 for (i = 1.0; i <= np/2; i++) {
4839 xx = xlow + (i-.5) * step;
4840 fland = TMath::Landau(xx,mpc,par[0]) / par[0];
4841 sum += fland * TMath::Gaus(x[0],xx,par[3]);
4843 xx = xupp - (i-.5) * step;
4844 fland = TMath::Landau(xx,mpc,par[0]) / par[0];
4845 sum += fland * TMath::Gaus(x[0],xx,par[3]);
4849 return (par[2] * step * sum * invsq2pi / par[3]);
4852 //_____________________________________________________________________________
4853 TF1 *AliTRDCalibraFit::LanGauFit(TH1 *his, Double_t *fitrange, Double_t *startvalues
4854 , Double_t *parlimitslo, Double_t *parlimitshi
4855 , Double_t *fitparams, Double_t *fiterrors
4856 , Double_t *chiSqr, Int_t *ndf)
4859 // Function for the fit
4863 Char_t funname[100];
4865 TF1 *ffitold = (TF1 *) gROOT->GetListOfFunctions()->FindObject(funname);
4870 TF1 *ffit = new TF1(funname,LanGauFun,fitrange[0],fitrange[1],4);
4871 ffit->SetParameters(startvalues);
4872 ffit->SetParNames("Width","MP","Area","GSigma");
4874 for (i = 0; i < 4; i++) {
4875 ffit->SetParLimits(i,parlimitslo[i],parlimitshi[i]);
4878 his->Fit(funname,"RB0"); // Fit within specified range, use ParLimits, do not plot
4880 ffit->GetParameters(fitparams); // Obtain fit parameters
4881 for (i = 0; i < 4; i++) {
4882 fiterrors[i] = ffit->GetParError(i); // Obtain fit parameter errors
4884 chiSqr[0] = ffit->GetChisquare(); // Obtain chi^2
4885 ndf[0] = ffit->GetNDF(); // Obtain ndf
4887 return (ffit); // Return fit function
4891 //_____________________________________________________________________________
4892 Int_t AliTRDCalibraFit::LanGauPro(Double_t *params, Double_t &maxx, Double_t &fwhm)
4895 // Function for the fit
4908 Int_t maxcalls = 10000;
4910 // Search for maximum
4911 p = params[1] - 0.1 * params[0];
4912 step = 0.05 * params[0];
4916 while ((l != lold) && (i < maxcalls)) {
4920 l = LanGauFun(&x,params);
4922 step = -step / 10.0;
4927 if (i == maxcalls) {
4933 // Search for right x location of fy
4934 p = maxx + params[0];
4940 while ( (l != lold) && (i < maxcalls) ) {
4945 l = TMath::Abs(LanGauFun(&x,params) - fy);
4959 // Search for left x location of fy
4961 p = maxx - 0.5 * params[0];
4967 while ((l != lold) && (i < maxcalls)) {
4971 l = TMath::Abs(LanGauFun(&x,params) - fy);
4973 step = -step / 10.0;
4978 if (i == maxcalls) {
4987 //_____________________________________________________________________________
4988 Double_t AliTRDCalibraFit::GausConstant(Double_t *x, Double_t *par)
4991 // Gaus with identical mean
4994 Double_t Gauss = par[0] * TMath::Gaus(x[0],0.0,par[1])+par[2];