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>
71 #include "AliCDBManager.h"
73 #include "AliTRDCalibraFit.h"
74 #include "AliTRDCalibraMode.h"
75 #include "AliTRDCalibraVector.h"
76 #include "AliTRDcalibDB.h"
77 #include "AliTRDgeometry.h"
78 #include "AliTRDCommonParam.h"
79 #include "./Cal/AliTRDCalROC.h"
80 #include "./Cal/AliTRDCalPad.h"
81 #include "./Cal/AliTRDCalDet.h"
84 ClassImp(AliTRDCalibraFit)
86 AliTRDCalibraFit* AliTRDCalibraFit::fgInstance = 0;
87 Bool_t AliTRDCalibraFit::fgTerminated = kFALSE;
89 //_____________singleton implementation_________________________________________________
90 AliTRDCalibraFit *AliTRDCalibraFit::Instance()
93 // Singleton implementation
96 if (fgTerminated != kFALSE) {
100 if (fgInstance == 0) {
101 fgInstance = new AliTRDCalibraFit();
108 //______________________________________________________________________________________
109 void AliTRDCalibraFit::Terminate()
112 // Singleton implementation
113 // Deletes the instance of this class
116 fgTerminated = kTRUE;
118 if (fgInstance != 0) {
125 //______________________________________________________________________________________
126 AliTRDCalibraFit::AliTRDCalibraFit()
131 ,fFitLagrPolOn(kFALSE)
132 ,fTakeTheMaxPH(kFALSE)
135 ,fBeginFitCharge(0.0)
141 ,fMeanChargeOn(kFALSE)
142 ,fFitChargeBisOn(kFALSE)
143 ,fFitChargeOn(kFALSE)
150 ,fNumberFitSuccess(0)
167 ,fScaleFitFactor(0.0)
173 // Default constructor
176 fCalibraMode = new AliTRDCalibraMode();
179 for (Int_t i = 0; i < 3; i++) {
180 fWriteCoef[i] = kFALSE;
184 for (Int_t k = 0; k < 3; k++) {
188 for (Int_t i = 0; i < 3; i++) {
197 //______________________________________________________________________________________
198 AliTRDCalibraFit::AliTRDCalibraFit(const AliTRDCalibraFit &c)
203 ,fFitLagrPolOn(kFALSE)
204 ,fTakeTheMaxPH(kFALSE)
207 ,fBeginFitCharge(0.0)
213 ,fMeanChargeOn(kFALSE)
214 ,fFitChargeBisOn(kFALSE)
215 ,fFitChargeOn(kFALSE)
222 ,fNumberFitSuccess(0)
239 ,fScaleFitFactor(0.0)
250 //____________________________________________________________________________________
251 AliTRDCalibraFit::~AliTRDCalibraFit()
254 // AliTRDCalibraFit destructor
261 //_____________________________________________________________________________
262 void AliTRDCalibraFit::Destroy()
275 //_____________________________________________________________________________
276 void AliTRDCalibraFit::ClearTree()
301 //_____________________________________________________________________________
302 void AliTRDCalibraFit::Init()
305 // Init some default values
309 fWriteNameCoef = "TRD.coefficient.root";
313 fBeginFitCharge = 3.5;
318 // Internal variables
320 // Variables in the loop
321 for (Int_t k = 0; k < 4; k++) {
322 fChargeCoef[k] = 1.0;
323 fVdriftCoef[k] = 1.5;
326 fChargeCoef[4] = 1.0;
327 for (Int_t i = 0; i < 3; i++) {
331 // Local database to be changed
336 //____________Functions fit Online CH2d________________________________________
337 Bool_t AliTRDCalibraFit::FitCHOnline(TH2I *ch)
340 // Fit the 1D histos, projections of the 2D ch on the Xaxis, for each
341 // calibration group normalized the resulted coefficients (to 1 normally)
342 // and write the results in a tree
346 if((fFitChargeNDB == 0) && (!fFitChargeOn)){
347 AliInfo("You have choosen to write the default fit method but it is not on!");
350 if((fFitChargeNDB == 1) && (!fMeanChargeOn)){
351 AliInfo("You have choosen to write the mean method but it is not on!");
354 if((fFitChargeNDB == 2) && (!fFitChargeBisOn)){
355 AliInfo("You have choosen to write the second fit method but it is not on!");
358 if((fFitChargeNDB == 4) && (!fFitMeanWOn)){
359 AliInfo("You have choosen to write the mean w method but it is not on!");
364 // Number of Xbins (detectors or groups of pads)
365 TAxis *xch = ch->GetXaxis();
366 Int_t nbins = xch->GetNbins();
367 TAxis *yph = ch->GetYaxis();
368 Int_t nybins = yph->GetNbins();
369 if (!InitFit(nbins,0)) {
372 fStatisticMean = 0.0;
374 fNumberFitSuccess = 0;
377 // Init fCountDet and fCount
378 InitfCountDetAndfCount(0);
380 // Beginning of the loop betwwen dect1 and dect2
381 for (Int_t idect = fDect1[0]; idect < fDect2[0]; idect++) {
383 TH1I *projch = (TH1I *) ch->ProjectionY("projch",idect+1,idect+1,(Option_t *)"e");
384 projch->SetDirectory(0);
386 // Determination of fNnZ, fNnRphi, fNfragZ and fNfragRphi
387 UpdatefCountDetAndfCount(idect,0);
389 // Reconstruction of the row and pad group: rowmin, row max ...
390 ReconstructFitRowMinRowMax(idect, 0);
392 // Number of entries for this calibration group
393 Double_t nentries = 0.0;
395 for (Int_t k = 0; k < nybins; k++) {
396 nentries += ch->GetBinContent(ch->GetBin(idect+1,k+1));
397 mean += projch->GetBinCenter(k+1)*projch->GetBinContent(k+1);
405 // Rebin and statistic stuff
408 projch = ReBin((TH1I *) projch);
410 // This detector has not enough statistics or was off
411 if (nentries < fMinEntries) {
412 // Fill with the default infos
413 NotEnoughStatistic(idect,0);
421 // Statistics of the group fitted
422 AliInfo(Form("For the group number %d there are %f stats",idect,nentries));
423 fStatisticMean += nentries;
427 // Method Mean and fit
428 // idect is egal for fDebug = 0 and 2, only to fill the hist
429 fChargeCoef[1] = mean;
431 FitMean((TH1 *) projch,(Int_t) (idect-fDect1[0]),nentries);
434 FitCH((TH1 *) projch,(Int_t) (idect-fDect1[0]));
436 if(fFitChargeBisOn) {
437 FitBisCH((TH1 *) projch,(Int_t) (idect-fDect1[0]));
440 FitMeanW((TH1 *) projch,(Int_t) (idect-fDect1[0]));
443 // Visualise the detector for fDebug 3 or 4
444 // Here is the reconstruction of the pad and row group is used!
449 FillInfosFit(idect,0);
466 // 0 no plot, 1 and 4 error plot, 3 and 4 DB plot
477 if (fNumberFit > 0) {
478 AliInfo(Form("There are %d with at least one entries.",fNumberEnt));
479 AliInfo(Form("%d fits have been proceeded (sucessfully or not...).",fNumberFit));
480 AliInfo(Form("There is a mean statistic of: %d over these fitted histograms and %d successfulled fits"
481 , (Int_t) fStatisticMean/fNumberFit, fNumberFitSuccess));
482 fStatisticMean = fStatisticMean / fNumberFit;
485 AliInfo(Form("There are %d with at least one entries. There is no fit!",fNumberEnt));
489 ConvertVectorFitCHTree();
498 //____________Functions fit Online CH2d________________________________________
499 Bool_t AliTRDCalibraFit::FitCHOnline()
502 // Reconstruct a 1D histo from the vectorCH for each calibration group,
503 // fit the histo, normalized the resulted coefficients (to 1 normally)
504 // and write the results in a tree
508 if((fFitChargeNDB == 0) && (!fFitChargeOn)){
509 AliInfo("You have choosen to write the default fit method but it is not on!");
512 if((fFitChargeNDB == 1) && (!fMeanChargeOn)){
513 AliInfo("You have choosen to write the mean method but it is not on!");
516 if((fFitChargeNDB == 2) && (!fFitChargeBisOn)){
517 AliInfo("You have choosen to write the second fit method but it is not on!");
520 if((fFitChargeNDB == 4) && (!fFitMeanWOn)){
521 AliInfo("You have choosen to write the mean w method but it is not on!");
527 if (!fCalibraVector) {
528 AliError("You have first to set the calibravector before using this function!");
532 // Number of Xbins (detectors or groups of pads)
536 fStatisticMean = 0.0;
538 fNumberFitSuccess = 0;
541 // Init fCountDet and fCount
542 InitfCountDetAndfCount(0);
544 // Beginning of the loop between dect1 and dect2
545 for (Int_t idect = fDect1[0]; idect < fDect2[0]; idect++) {
547 // Search if the group is in the VectorCH
548 Int_t place = fCalibraVector->SearchInVector(idect,0);
555 projch = fCalibraVector->ConvertVectorCTHisto(place,(const char *) name);
556 projch->SetDirectory(0);
559 // Determination of fNnZ, fNnRphi, fNfragZ and fNfragRphi
560 UpdatefCountDetAndfCount(idect,0);
562 // Reconstruction of the row and pad group: rowmin, row max ...
563 ReconstructFitRowMinRowMax(idect,0);
565 // Number of entries and mean
566 Double_t nentries = 0.0;
569 for (Int_t k = 0; k < fCalibraVector->GetNumberBinCharge(); k++) {
570 nentries += projch->GetBinContent(k+1);
571 mean += projch->GetBinCenter(k+1)*projch->GetBinContent(k+1);
579 // Rebin and statistic stuff
583 projch = ReBin((TH1F *) projch);
586 // This detector has not enough statistics or was not found in VectorCH
589 (nentries < fMinEntries))) {
591 // Fill with the default infos
592 NotEnoughStatistic(idect,0);
603 // Statistic of the histos fitted
604 AliInfo(Form("For the group number %d there are %f stats",idect,nentries));
605 fStatisticMean += nentries;
609 // Method Mean and fit
610 // idect is egal for fDebug = 0 and 2, only to fill the hist
611 fChargeCoef[1] = mean;
613 FitMean((TH1 *) projch,(Int_t) (idect-fDect1[0]),nentries);
616 FitCH((TH1 *) projch,(Int_t) (idect-fDect1[0]));
618 if(fFitChargeBisOn) {
619 FitBisCH((TH1 *) projch,(Int_t) (idect-fDect1[0]));
622 FitMeanW((TH1 *) projch,(Int_t) (idect-fDect1[0]));
625 // Visualise the detector for fDebug 3 or 4
626 // Here is the reconstruction of the pad and row group is used!
632 FillInfosFit(idect,0);
649 // 0 no plot, 1 and 4 error plot, 3 and 4 DB plot
660 if (fNumberFit > 0) {
661 AliInfo(Form("There are %d with at least one entries.",fNumberEnt));
662 AliInfo(Form("%d fits have been proceeded (sucessfully or not...).",fNumberFit));
663 AliInfo(Form("There is a mean statistic of: %d over these fitted histograms and %d successfulled fits"
664 ,(Int_t) fStatisticMean/fNumberFit, fNumberFitSuccess));
665 fStatisticMean = fStatisticMean / fNumberFit;
668 AliInfo(Form("There are %d with at least one entries. There is no fit!",fNumberEnt));
672 ConvertVectorFitCHTree();
681 //____________Functions fit Online CH2d________________________________________
682 Bool_t AliTRDCalibraFit::FitCHOnline(TTree *tree)
685 // Look if the calibration group can be found in the tree, if yes take the
686 // histo, fit it, normalized the resulted coefficients (to 1 normally) and
687 // write the results in a tree
691 if((fFitChargeNDB == 0) && (!fFitChargeOn)){
692 AliInfo("You have choosen to write the default fit method but it is not on!");
695 if((fFitChargeNDB == 1) && (!fMeanChargeOn)){
696 AliInfo("You have choosen to write the mean method but it is not on!");
699 if((fFitChargeNDB == 2) && (!fFitChargeBisOn)){
700 AliInfo("You have choosen to write the second fit method but it is not on!");
703 if((fFitChargeNDB == 4) && (!fFitMeanWOn)){
704 AliInfo("You have choosen to write the mean w method but it is not on!");
709 // Number of Xbins (detectors or groups of pads)
713 fStatisticMean = 0.0;
715 fNumberFitSuccess = 0;
719 fCalibraVector = new AliTRDCalibraVector();
722 // Init fCountDet and fCount
723 InitfCountDetAndfCount(0);
725 tree->SetBranchAddress("histo",&projch);
726 TObjArray *vectorplace = fCalibraVector->ConvertTreeVector(tree);
728 // Beginning of the loop between dect1 and dect2
729 for (Int_t idect = fDect1[0]; idect < fDect2[0]; idect++) {
731 //Search if the group is in the VectorCH
732 Int_t place = fCalibraVector->SearchInTreeVector(vectorplace,idect);
737 tree->GetEntry(place);
740 // Determination of fNnZ, fNnRphi, fNfragZ and fNfragRphi
741 UpdatefCountDetAndfCount(idect,0);
743 // Reconstruction of the row and pad group: rowmin, row max ...
744 ReconstructFitRowMinRowMax(idect,0);
746 // Number of entries and mean
747 Double_t nentries = 0.0;
750 for (Int_t k = 0; k < projch->GetXaxis()->GetNbins(); k++) {
751 nentries += projch->GetBinContent(k+1);
752 mean += projch->GetBinCenter(k+1)*projch->GetBinContent(k+1);
761 // Rebin and statistic stuff
765 projch = ReBin((TH1F *) projch);
768 // This detector has not enough statistics or was not found in VectorCH
771 (nentries < fMinEntries))) {
773 // Fill with the default infos
774 NotEnoughStatistic(idect,0);
780 // Statistics of the group fitted
781 AliInfo(Form("For the group number %d there are %f stats",idect,nentries));
783 fStatisticMean += nentries;
785 // Method Mean and fit
786 // idect is egal for fDebug = 0 and 2, only to fill the hist
787 fChargeCoef[1] = mean;
789 FitMean((TH1 *) projch,(Int_t) (idect-fDect1[0]),nentries);
792 FitCH((TH1 *) projch,(Int_t) (idect-fDect1[0]));
794 if(fFitChargeBisOn) {
795 FitBisCH((TH1 *) projch,(Int_t) (idect-fDect1[0]));
798 FitMeanW((TH1 *) projch,(Int_t) (idect-fDect1[0]));
801 // Visualise the detector for fDebug 3 or 4
802 // Here is the reconstruction of the pad and row group is used!
808 FillInfosFit(idect,0);
820 // 0 no plot, 1 and 4 error plot, 3 and 4 DB plot
831 if (fNumberFit > 0) {
832 AliInfo(Form("There are %d with at least one entries.",fNumberEnt));
833 AliInfo(Form("%d fits have been proceeded (sucessfully or not...).",fNumberFit));
834 AliInfo(Form("There is a mean statistic of: %d over these fitted histograms and %d successfulled fits"
835 ,(Int_t) fStatisticMean/fNumberFit,fNumberFitSuccess));
836 fStatisticMean = fStatisticMean / fNumberFit;
839 AliInfo(Form("There are %d with at least one entries. There is no fit!",fNumberEnt));
843 ConvertVectorFitCHTree();
853 //________________functions fit Online PH2d____________________________________
854 Bool_t AliTRDCalibraFit::FitPHOnline(TProfile2D *ph)
857 // Take the 1D profiles (average pulse height), projections of the 2D PH
858 // on the Xaxis, for each calibration group
859 // Fit or use the slope of the average pulse height to reconstruct the
860 // drift velocity write the results in a tree
861 // A first calibration of T0 is also made using the same method (slope method)
865 if((fFitPHNDB == 0) && (!fFitPHOn)){
866 AliInfo("You have choosen to write the fit method but it is not on!");
869 if((fFitPHNDB == 1) && (!fFitPol2On)){
870 AliInfo("You have choosen to write the Pol2 method but it is not on!");
873 if((fFitPHNDB == 3) && (!fFitLagrPolOn)){
874 AliInfo("You have choosen to write the LagrPol2 method but it is not on!");
878 // Number of Xbins (detectors or groups of pads)
879 TAxis *xph = ph->GetXaxis();
880 TAxis *yph = ph->GetYaxis();
881 Int_t nbins = xph->GetNbins();
882 Int_t nybins = yph->GetNbins();
883 if (!InitFit(nbins,1)) {
886 fStatisticMean = 0.0;
888 fNumberFitSuccess = 0;
891 // Init fCountDet and fCount
892 InitfCountDetAndfCount(1);
894 // Beginning of the loop
895 for (Int_t idect = fDect1[1]; idect < fDect2[1]; idect++) {
897 TH1D *projph = (TH1D *) ph->ProjectionY("projph",idect+1,idect+1,(Option_t *) "e");
898 projph->SetDirectory(0);
900 // Number of entries for this calibration group
901 Double_t nentries = 0;
902 for (Int_t k = 0; k < nybins; k++) {
903 nentries += ph->GetBinEntries(ph->GetBin(idect+1,k+1));
909 // Determination of fNnZ, fNnRphi, fNfragZ and fNfragRphi
910 UpdatefCountDetAndfCount(idect,1);
912 // Reconstruction of the row and pad group: rowmin, row max ...
913 ReconstructFitRowMinRowMax(idect,1);
915 // Rebin and statistic stuff
916 // This detector has not enough statistics or was off
917 if (nentries < fMinEntries) {
919 // Fill with the default values
920 NotEnoughStatistic(idect,1);
931 // Statistics of the histos fitted
932 AliInfo(Form("For the group number %d there are %f stats",idect,nentries));
934 fStatisticMean += nentries;
936 // Calcul of "real" coef
937 CalculVdriftCoefMean(fCountDet[1],(Int_t) (idect - fDect1[1]));
938 CalculT0CoefMean(fCountDet[1],(Int_t) (idect - fDect1[1]));
940 // Method Mean and fit
941 // idect is egal for fDebug = 0 and 2, only to fill the hist
943 FitPente((TH1 *) projph,(Int_t) (idect - fDect1[1]));
946 FitLagrangePoly((TH1 *) projph,(Int_t) (idect - fDect1[1]));
949 FitPH((TH1 *) projph,(Int_t) (idect - fDect1[1]));
952 // Visualise the detector for fDebug 3 or 4
953 // Here is the reconstruction of the pad and row group is used!
958 // Fill the tree if end of a detector or only the pointer to the branch!!!
959 FillInfosFit(idect,1);
969 // 0 no plot, 1 and 4 error plot, 3 and 4 DB plot
982 if (fNumberFit > 0) {
983 AliInfo(Form("There are %d with at least one entries.",fNumberEnt));
984 AliInfo(Form("%d fits have been proceeded (sucessfully or not...).",fNumberFit));
985 AliInfo(Form("There is a mean statistic of: %d over these fitted histograms and %d successfulled fits"
986 ,(Int_t) fStatisticMean/fNumberFit,fNumberFitSuccess));
987 fStatisticMean = fStatisticMean / fNumberFit;
990 AliInfo(Form("There are %d with at least one entries. There is no fit!",fNumberEnt));
1002 //____________Functions fit Online PH2d________________________________________
1003 Bool_t AliTRDCalibraFit::FitPHOnline()
1006 // Reconstruct the average pulse height from the vectorPH for each
1007 // calibration group
1008 // Fit or use the slope of the average pulse height to reconstruct the
1009 // drift velocity write the results in a tree
1010 // A first calibration of T0 is also made using the same method (slope method)
1014 if((fFitPHNDB == 0) && (!fFitPHOn)){
1015 AliInfo("You have choosen to write the fit method but it is not on!");
1018 if((fFitPHNDB == 1) && (!fFitPol2On)){
1019 AliInfo("You have choosen to write the Pol2 method but it is not on!");
1022 if((fFitPHNDB == 3) && (!fFitLagrPolOn)){
1023 AliInfo("You have choosen to write the LagrPol2 method but it is not on!");
1028 if (!fCalibraVector) {
1029 AliError("You have first to set the calibravector before using this function!");
1034 // Number of Xbins (detectors or groups of pads)
1035 if (!InitFit(0,1)) {
1038 fStatisticMean = 0.0;
1040 fNumberFitSuccess = 0;
1043 // Init fCountDet and fCount
1044 InitfCountDetAndfCount(1);
1046 // Beginning of the loop
1047 for (Int_t idect = fDect1[1]; idect < fDect2[1]; idect++) {
1049 // Search if the group is in the VectorCH
1050 Int_t place = fCalibraVector->SearchInVector(idect,1);
1059 projph = CorrectTheError((TGraphErrors *) (fCalibraVector->ConvertVectorPHisto(place,(const char *) name)));
1060 projph->SetDirectory(0);
1063 // Determination of fNnZ, fNnRphi, fNfragZ and fNfragRphi
1064 UpdatefCountDetAndfCount(idect,1);
1066 // Reconstruction of the row and pad group: rowmin, row max ...
1067 ReconstructFitRowMinRowMax(idect,1);
1069 // Rebin and statistic stuff
1070 // This detector has not enough statistics or was off
1071 if ((place == -1) ||
1073 (fEntriesCurrent < fMinEntries))) {
1075 // Fill with the default values
1076 NotEnoughStatistic(idect,1);
1087 // Statistic of the histos fitted
1088 AliInfo(Form("For the group number %d there are %d stats",idect,fEntriesCurrent));
1090 fStatisticMean += fEntriesCurrent;
1092 // Calcul of "real" coef
1093 CalculVdriftCoefMean(fCountDet[1],(Int_t) (idect - fDect1[1]));
1094 CalculT0CoefMean(fCountDet[1],(Int_t) (idect - fDect1[1]));
1096 // Method Mean and fit
1097 // idect is egal for fDebug = 0 and 2, only to fill the hist
1099 FitPente((TH1 *) projph,(Int_t) (idect - fDect1[1]));
1102 FitLagrangePoly((TH1 *) projph,(Int_t) (idect - fDect1[1]));
1105 FitPH((TH1 *) projph,(Int_t) (idect - fDect1[1]));
1108 // Visualise the detector for fDebug 3 or 4
1109 // Here is the reconstruction of the pad and row group is used!
1115 // Fill the tree if end of a detector or only the pointer to the branch!!!
1116 FillInfosFit(idect,1);
1127 // 0 no plot, 1 and 4 error plot, 3 and 4 DB plot
1128 if ((fDebug == 1) ||
1133 if ((fDebug == 4) ||
1140 if (fNumberFit > 0) {
1141 AliInfo(Form("There are %d with at least one entries.",fNumberEnt));
1142 AliInfo(Form("%d fits have been proceeded (sucessfully or not...).",fNumberFit));
1143 AliInfo(Form("There is a mean statistic of: %d over these fitted histograms and %d successfulled fits"
1144 ,(Int_t) fStatisticMean/fNumberFit,fNumberFitSuccess));
1145 fStatisticMean = fStatisticMean / fNumberFit;
1148 AliInfo(Form("There are %d with at least one entries. There is no fit!",fNumberEnt));
1151 // Write the things!
1152 if (fWriteCoef[1]) {
1160 //____________Functions fit Online PH2d________________________________________
1161 Bool_t AliTRDCalibraFit::FitPHOnline(TTree *tree)
1164 // Look if the calibration group can be found in the tree, if yes take the
1165 // histo, fit it, and write the results in a tree
1166 // A first calibration of T0 is also made using the same method (slope method)
1170 if ((fFitPHNDB == 0) && (!fFitPHOn)){
1171 AliInfo("You have choosen to write the fit method but it is not on!");
1174 if ((fFitPHNDB == 1) && (!fFitPol2On)){
1175 AliInfo("You have choosen to write the Pol2 method but it is not on!");
1178 if ((fFitPHNDB == 3) && (!fFitLagrPolOn)){
1179 AliInfo("You have choosen to write the LagrPol2 method but it is not on!");
1183 // Number of Xbins (detectors or groups of pads)
1184 if (!InitFit(0,1)) {
1187 fStatisticMean = 0.0;
1189 fNumberFitSuccess = 0;
1193 fCalibraVector = new AliTRDCalibraVector();
1195 // Init fCountDet and fCount
1196 InitfCountDetAndfCount(1);
1197 TGraphErrors *projphtree = 0x0;
1198 tree->SetBranchAddress("histo",&projphtree);
1199 TObjArray *vectorplace = fCalibraVector->ConvertTreeVector(tree);
1201 // Beginning of the loop
1202 for (Int_t idect = fDect1[1]; idect < fDect2[1]; idect++) {
1204 // Search if the group is in the VectorCH
1205 Int_t place = fCalibraVector->SearchInTreeVector(vectorplace,idect);
1213 tree->GetEntry(place);
1214 projph = CorrectTheError(projphtree);
1217 // Determination of fNnZ, fNnRphi, fNfragZ and fNfragRphi
1218 UpdatefCountDetAndfCount(idect,1);
1220 // Reconstruction of the row and pad group: rowmin, row max ...
1221 ReconstructFitRowMinRowMax(idect,1);
1223 // Rebin and statistic stuff
1224 // This detector has not enough statistics or was off
1227 (fEntriesCurrent < fMinEntries))) {
1229 // Fill with the default values
1230 NotEnoughStatistic(idect,1);
1241 // Statistics of the histos fitted
1242 AliInfo(Form("For the group number %d there are %d stats",idect,fEntriesCurrent));
1244 fStatisticMean += fEntriesCurrent;
1246 // Calcul of "real" coef
1247 CalculVdriftCoefMean(fCountDet[1],(Int_t) (idect - fDect1[1]));
1248 CalculT0CoefMean(fCountDet[1],(Int_t) (idect - fDect1[1]));
1250 // Method Mean and fit
1251 // idect is egal for fDebug = 0 and 2, only to fill the hist
1253 FitPente((TH1 *) projph,(Int_t) (idect - fDect1[1]));
1256 FitLagrangePoly((TH1 *) projph,(Int_t) (idect - fDect1[1]));
1259 FitPH((TH1 *) projph,(Int_t) (idect - fDect1[1]));
1262 // Visualise the detector for fDebug 3 or 4
1263 // Here is the reconstruction of the pad and row group is used!
1269 // Fill the tree if end of a detector or only the pointer to the branch!!!
1270 FillInfosFit(idect,1);
1280 // 0 no plot, 1 and 4 error plot, 3 and 4 DB plot
1281 if ((fDebug == 1) ||
1286 if ((fDebug == 4) ||
1293 if (fNumberFit > 0) {
1294 AliInfo(Form("There are %d with at least one entries.",fNumberEnt));
1295 AliInfo(Form("%d fits have been proceeded (sucessfully or not...).",fNumberFit));
1296 AliInfo(Form("There is a mean statistic of: %d over these fitted histograms and %d successfulled fits"
1297 ,(Int_t) fStatisticMean/fNumberFit,fNumberFitSuccess));
1298 fStatisticMean = fStatisticMean / fNumberFit;
1301 AliInfo(Form("There are %d with at least one entries. There is no fit!",fNumberEnt));
1304 // Write the things!
1305 if (fWriteCoef[1]) {
1313 //____________Functions fit Online PRF2d_______________________________________
1314 Bool_t AliTRDCalibraFit::FitPRFOnline(TProfile2D *prf)
1317 // Take the 1D profiles (pad response function), projections of the 2D PRF
1318 // on the Xaxis, for each calibration group
1319 // Fit with a gaussian to reconstruct the sigma of the pad response function
1320 // write the results in a tree
1324 if ((fFitPRFNDB == 2) && (!fRMSPRFOn)){
1325 AliInfo("You have choosen to write the RMS method but it is not on!");
1328 if ((fFitPRFNDB == 0) && (!fFitPRFOn)){
1329 AliInfo("You have choosen to write the fit method but it is not on!");
1333 // Number of Xbins (detectors or groups of pads)
1334 TAxis *xprf = prf->GetXaxis();
1335 TAxis *yprf = prf->GetYaxis();
1336 Int_t nybins = yprf->GetNbins();
1337 Int_t nbins = xprf->GetNbins();
1338 if (!InitFit(nbins,2)) {
1341 fStatisticMean = 0.0;
1343 fNumberFitSuccess = 0;
1346 // Init fCountDet and fCount
1347 InitfCountDetAndfCount(2);
1349 // Beginning of the loop
1350 for (Int_t idect = fDect1[2]; idect < fDect2[2]; idect++) {
1352 TH1D *projprf = (TH1D *) prf->ProjectionY("projprf",idect+1,idect+1,(Option_t *) "e");
1353 projprf->SetDirectory(0);
1355 // Number of entries for this calibration group
1356 Double_t nentries = 0;
1357 for (Int_t k = 0; k < nybins; k++) {
1358 nentries += prf->GetBinEntries(prf->GetBin(idect+1,k+1));
1360 if(nentries > 0) fNumberEnt++;
1362 // Determination of fNnZ, fNnRphi, fNfragZ and fNfragRphi
1363 UpdatefCountDetAndfCount(idect,2);
1365 // Reconstruction of the row and pad group: rowmin, row max ...
1366 ReconstructFitRowMinRowMax(idect,2);
1368 // Rebin and statistic stuff
1369 // This detector has not enough statistics or was off
1370 if (nentries < fMinEntries) {
1372 // Fill with the default values
1373 NotEnoughStatistic(idect,2);
1384 // Statistics of the histos fitted
1385 AliInfo(Form("For the group number %d there are %f stats",idect,nentries));
1387 fStatisticMean += nentries;
1389 // Calcul of "real" coef
1390 if ((fDebug == 1) ||
1392 CalculPRFCoefMean(fCountDet[2],(Int_t) (idect - fDect1[2]));
1395 // Method Mean and fit
1396 // idect is egal for fDebug = 0 and 2, only to fill the hist
1398 FitPRF((TH1 *) projprf,(Int_t) (idect - fDect1[2]));
1401 RmsPRF((TH1 *) projprf,(Int_t) (idect - fDect1[2]));
1405 // Visualise the detector for fDebug 3 or 4
1406 // Here is the reconstruction of the pad and row group is used!
1411 // Fill the tree if end of a detector or only the pointer to the branch!!!
1412 FillInfosFit(idect,2);
1422 // No plot, 1 and 4 error plot, 3 and 4 DB plot
1423 if ((fDebug == 1) ||
1427 if ((fDebug == 4) ||
1433 if (fNumberFit > 0) {
1434 AliInfo(Form("There are %d with at least one entries.",fNumberEnt));
1435 AliInfo(Form("%d fits have been proceeded (sucessfully or not...).",fNumberFit));
1436 AliInfo(Form("There is a mean statistic of: %d over these fitted histograms and %d successfulled fits"
1437 ,(Int_t) fStatisticMean/fNumberFit,fNumberFitSuccess));
1438 fStatisticMean = fStatisticMean / fNumberFit;
1441 AliInfo(Form("There are %d with at least one entries. There is no fit!",fNumberEnt));
1444 // Write the things!
1445 if (fWriteCoef[2]) {
1453 //____________Functions fit Online PRF2d_______________________________________
1454 Bool_t AliTRDCalibraFit::FitPRFOnline()
1457 // Reconstruct the 1D histo (pad response function) from the vectorPRD for
1458 // each calibration group
1459 // Fit with a gaussian to reconstruct the sigma of the pad response function
1460 // write the results in a tree
1464 if ((fFitPRFNDB == 2) && (!fRMSPRFOn)){
1465 AliInfo("You have choosen to write the RMS method but it is not on!");
1468 if ((fFitPRFNDB == 0) && (!fFitPRFOn)){
1469 AliInfo("You have choosen to write the fit method but it is not on!");
1474 if (!fCalibraVector) {
1475 AliError("You have first to set the calibravector before using this function!");
1479 // Number of Xbins (detectors or groups of pads)
1480 if (!InitFit(0,2)) {
1483 fStatisticMean = 0.0;
1485 fNumberFitSuccess = 0;
1488 // Init fCountDet and fCount
1489 InitfCountDetAndfCount(2);
1491 // Beginning of the loop
1492 for (Int_t idect = fDect1[2]; idect < fDect2[2]; idect++) {
1494 // Search if the group is in the VectorCH
1495 Int_t place = fCalibraVector->SearchInVector(idect,2);
1498 TH1F *projprf = 0x0;
1499 TString name("PRF");
1504 projprf = CorrectTheError((TGraphErrors *) (fCalibraVector->ConvertVectorPHisto(place,(const char *)name)));
1505 projprf->SetDirectory(0);
1508 // Determination of fNnZ, fNnRphi, fNfragZ and fNfragRphi
1509 UpdatefCountDetAndfCount(idect,2);
1511 // Reconstruction of the row and pad group: rowmin, row max ...
1512 ReconstructFitRowMinRowMax(idect,2);
1514 // Rebin and statistic stuff
1515 // This detector has not enough statistics or was off
1516 if ((place == -1) ||
1518 (fEntriesCurrent < fMinEntries))) {
1520 // Fill with the default values
1521 NotEnoughStatistic(idect,2);
1532 // Statistic of the histos fitted
1533 AliInfo(Form("For the group number %d there are %d stats", idect,fEntriesCurrent));
1535 fStatisticMean += fEntriesCurrent;
1537 // Calcul of "real" coef
1538 if ((fDebug == 1) ||
1540 CalculPRFCoefMean(fCountDet[2],(Int_t) (idect-fDect1[2]));
1543 // Method Mean and fit
1544 // idect is egal for fDebug = 0 and 2, only to fill the hist
1546 FitPRF((TH1 *) projprf,(Int_t) (idect-fDect1[2]));
1549 RmsPRF((TH1 *) projprf,(Int_t) (idect-fDect1[2]));
1552 // Visualise the detector for fDebug 3 or 4
1553 // Here is the reconstruction of the pad and row group is used!
1557 // Fill the tree if end of a detector or only the pointer to the branch!!!
1558 FillInfosFit(idect,2);
1568 // No plot, 1 and 4 error plot, 3 and 4 DB plot
1569 if ((fDebug == 1) ||
1573 if ((fDebug == 4) ||
1579 if (fNumberFit > 0) {
1580 AliInfo(Form("There are %d with at least one entries.",fNumberEnt));
1581 AliInfo(Form("%d fits have been proceeded (sucessfully or not...).",fNumberFit));
1582 AliInfo(Form("There is a mean statistic of: %d over these fitted histograms and %d successfulled fits"
1583 ,(Int_t) fStatisticMean/fNumberFit,fNumberFitSuccess));
1586 AliInfo(Form("There are %d with at least one entries. There is no fit!",fNumberEnt));
1589 // Write the things!
1590 if (fWriteCoef[2]) {
1598 //____________Functions fit Online PRF2d_______________________________________
1599 Bool_t AliTRDCalibraFit::FitPRFOnline(TTree *tree)
1602 // Look if the calibration group can be found in the tree, if yes take
1603 // the histo, fit it, and write the results in a tree
1607 if ((fFitPRFNDB == 2) && (!fRMSPRFOn)){
1608 AliInfo("You have choosen to write the RMS method but it is not on!");
1611 if ((fFitPRFNDB == 0) && (!fFitPRFOn)){
1612 AliInfo("You have choosen to write the fit method but it is not on!");
1616 // Number of Xbins (detectors or groups of pads)
1617 if (!InitFit(0,2)) {
1620 fStatisticMean = 0.0;
1622 fNumberFitSuccess = 0;
1626 fCalibraVector = new AliTRDCalibraVector();
1628 // Init fCountDet and fCount
1629 InitfCountDetAndfCount(2);
1630 TGraphErrors *projprftree = 0x0;
1631 tree->SetBranchAddress("histo",&projprftree);
1632 TObjArray *vectorplace = fCalibraVector->ConvertTreeVector(tree);
1634 // Beginning of the loop
1635 for (Int_t idect = fDect1[2]; idect < fDect2[2]; idect++) {
1637 // Search if the group is in the VectorCH
1638 Int_t place = fCalibraVector->SearchInTreeVector(vectorplace,idect);
1641 TH1F *projprf = 0x0;
1646 tree->GetEntry(place);
1647 projprf = CorrectTheError(projprftree);
1650 // Determination of fNnZ, fNnRphi, fNfragZ and fNfragRphi
1651 UpdatefCountDetAndfCount(idect,2);
1653 // Reconstruction of the row and pad group: rowmin, row max ...
1654 ReconstructFitRowMinRowMax(idect,2);
1656 // Rebin and statistic stuff
1657 // This detector has not enough statistics or was off
1658 if ((place == -1) ||
1660 (fEntriesCurrent < fMinEntries))) {
1662 // Fill with the default values
1663 NotEnoughStatistic(idect,2);
1674 // Statistics of the histos fitted
1675 AliInfo(Form("For the group number %d there are %d stats",idect,fEntriesCurrent));
1677 fStatisticMean += fEntriesCurrent;
1679 // Calcul of "real" coef
1680 if ((fDebug == 1) ||
1682 CalculPRFCoefMean(fCountDet[2],(Int_t) (idect - fDect1[2]));
1685 // Method Mean and fit
1686 // idect is egal for fDebug = 0 and 2, only to fill the hist
1688 FitPRF((TH1 *) projprf,(Int_t) (idect - fDect1[2]));
1691 RmsPRF((TH1 *) projprf,(Int_t) (idect - fDect1[2]));
1694 // Visualise the detector for fDebug 3 or 4
1695 // Here is the reconstruction of the pad and row group is used!
1699 // Fill the tree if end of a detector or only the pointer to the branch!!!
1700 FillInfosFit(idect,2);
1710 // No plot, 1 and 4 error plot, 3 and 4 DB plot
1711 if ((fDebug == 1) ||
1715 if ((fDebug == 4) ||
1721 if (fNumberFit > 0) {
1722 AliInfo(Form("There are %d with at least one entries.",fNumberEnt));
1723 AliInfo(Form("%d fits have been proceeded (sucessfully or not...).",fNumberFit));
1724 AliInfo(Form("There is a mean statistic of: %d over these fitted histograms and %d successfulled fits"
1725 ,(Int_t) fStatisticMean/fNumberFit,fNumberFitSuccess));
1726 fStatisticMean = fStatisticMean / fNumberFit;
1729 AliInfo(Form("There are %d with at least one entries. There is no fit!",fNumberEnt));
1732 // Write the things!
1733 if (fWriteCoef[2]) {
1741 //____________Functions for seeing if the pad is really okey___________________
1743 //_____________________________________________________________________________
1744 Bool_t AliTRDCalibraFit::SetModeCalibrationFromTObject(TObject *object, Int_t i)
1747 // Set fNz[i] and fNrphi[i] of the AliTRDCalibraFit::Instance()
1748 // corresponding to the given TObject
1751 const char *nametitle = object->GetTitle();
1754 const Char_t *patternz0 = "Nz0";
1755 const Char_t *patternz1 = "Nz1";
1756 const Char_t *patternz2 = "Nz2";
1757 const Char_t *patternz3 = "Nz3";
1758 const Char_t *patternz4 = "Nz4";
1759 const Char_t *patternrphi0 = "Nrphi0";
1760 const Char_t *patternrphi1 = "Nrphi1";
1761 const Char_t *patternrphi2 = "Nrphi2";
1762 const Char_t *patternrphi3 = "Nrphi3";
1763 const Char_t *patternrphi4 = "Nrphi4";
1764 const Char_t *patternrphi5 = "Nrphi5";
1765 const Char_t *patternrphi6 = "Nrphi6";
1768 UShort_t testrphi = 0;
1771 if (strstr(nametitle,patternz0)) {
1773 fCalibraMode->SetNz(i, 0);
1775 if (strstr(nametitle,patternz1)) {
1777 fCalibraMode->SetNz(i ,1);
1779 if (strstr(nametitle,patternz2)) {
1781 fCalibraMode->SetNz(i ,2);
1783 if (strstr(nametitle,patternz3)) {
1785 fCalibraMode->SetNz(i ,3);
1787 if (strstr(nametitle,patternz4)) {
1789 fCalibraMode->SetNz(i ,4);
1793 if (strstr(nametitle,patternrphi0)) {
1795 fCalibraMode->SetNrphi(i ,0);
1797 if (strstr(nametitle,patternrphi1)) {
1799 fCalibraMode->SetNrphi(i, 1);
1801 if (strstr(nametitle,patternrphi2)) {
1803 fCalibraMode->SetNrphi(i, 2);
1805 if (strstr(nametitle,patternrphi3)) {
1807 fCalibraMode->SetNrphi(i, 3);
1809 if (strstr(nametitle,patternrphi4)) {
1811 fCalibraMode->SetNrphi(i, 4);
1813 if (strstr(nametitle,patternrphi5)) {
1815 fCalibraMode->SetNrphi(i, 5);
1817 if (strstr(nametitle,patternrphi6)) {
1819 fCalibraMode->SetNrphi(i, 6);
1822 // Look if all is okey
1828 fCalibraMode->SetNrphi(i ,0);
1829 fCalibraMode->SetNz(i ,0);
1835 //_____________________________________________________________________________
1836 AliTRDCalDet *AliTRDCalibraFit::CreateDetObjectTree(TTree *tree, Int_t i)
1839 // It creates the AliTRDCalDet object from the tree of the coefficient
1840 // for the calibration i (i != 2)
1841 // It takes the mean value of the coefficients per detector
1842 // This object has to be written in the database
1845 // Create the DetObject
1846 AliTRDCalDet *object = 0x0;
1848 object = new AliTRDCalDet("ChamberGainFactor","GainFactor (detector value)");
1851 object = new AliTRDCalDet("ChamberVdrift","TRD drift velocities (detector value)");
1854 object = new AliTRDCalDet("ChamberT0","T0 (detector value)");
1858 Int_t detector = -1;
1859 Float_t values[2304];
1860 tree->SetBranchAddress("detector",&detector);
1862 tree->SetBranchAddress("gainPad",values);
1865 tree->SetBranchAddress("vdrift" ,values);
1868 tree->SetBranchAddress("t0" ,values);
1871 // For calculating the mean
1874 Int_t numberofentries = tree->GetEntries();
1876 if (numberofentries != 540) {
1877 AliInfo("The tree is not complete");
1880 for (Int_t det = 0; det < numberofentries; ++det) {
1881 tree->GetEntry(det);
1882 if (GetChamber(detector) == 2) {
1890 for (Int_t k = 0; k < nto; k++) {
1891 mean += TMath::Abs(values[k]) / nto;
1895 for (Int_t k = 0; k < nto; k++) {
1896 if(k == 0) mean = values[k];
1897 if(mean > values[k]) mean = values[k];
1900 object->SetValue(detector,mean);
1907 //_____________________________________________________________________________
1908 TObject *AliTRDCalibraFit::CreatePadObjectTree(TTree *tree, Int_t i
1909 , AliTRDCalDet *detobject)
1912 // It Creates the AliTRDCalPad object from the tree of the
1913 // coefficient for the calibration i (i != 2)
1914 // You need first to create the object for the detectors,
1915 // where the mean value is put.
1916 // This object has to be written in the database
1919 // Create the DetObject
1920 AliTRDCalPad *object = 0x0;
1922 object = new AliTRDCalPad("GainFactor","GainFactor (local variations)");
1925 object = new AliTRDCalPad("LocalVdrift","TRD drift velocities (local variations)");
1928 object = new AliTRDCalPad("LocalT0","T0 (local variations)");
1932 Int_t detector = -1;
1933 Float_t values[2304];
1934 tree->SetBranchAddress("detector",&detector);
1936 tree->SetBranchAddress("gainPad",values);
1939 tree->SetBranchAddress("vdrift" ,values);
1942 tree->SetBranchAddress("t0" ,values);
1947 Int_t numberofentries = tree->GetEntries();
1949 if (numberofentries != 540) {
1950 AliInfo("The tree is not complete");
1953 for (Int_t det = 0; det < numberofentries; ++det) {
1954 tree->GetEntry(det);
1955 AliTRDCalROC *calROC = object->GetCalROC(detector);
1956 mean = detobject->GetValue(detector);
1957 if ((mean == 0) && (i != 3)) {
1960 Int_t rowMax = calROC->GetNrows();
1961 Int_t colMax = calROC->GetNcols();
1962 for (Int_t row = 0; row < rowMax; ++row) {
1963 for (Int_t col = 0; col < colMax; ++col) {
1964 if(i != 3) calROC->SetValue(col,row,TMath::Abs(values[(Int_t) (col*rowMax+row)])/mean);
1965 else calROC->SetValue(col,row,values[(Int_t) (col*rowMax+row)]-mean);
1975 //_____________________________________________________________________________
1976 TObject *AliTRDCalibraFit::CreatePadObjectTree(TTree *tree)
1979 // It Creates the AliTRDCalPad object from the tree of the
1980 // coefficient for the calibration PRF (i = 2)
1981 // This object has to be written in the database
1984 // Create the DetObject
1985 AliTRDCalPad *object = new AliTRDCalPad("PRFWidth","PRFWidth");
1988 Int_t detector = -1;
1989 Float_t values[2304];
1990 tree->SetBranchAddress("detector",&detector);
1991 tree->SetBranchAddress("width" ,values);
1994 Int_t numberofentries = tree->GetEntries();
1996 if (numberofentries != 540) {
1997 AliInfo("The tree is not complete");
2000 for (Int_t det = 0; det < numberofentries; ++det) {
2001 tree->GetEntry(det);
2002 AliTRDCalROC *calROC = object->GetCalROC(detector);
2003 Int_t rowMax = calROC->GetNrows();
2004 Int_t colMax = calROC->GetNcols();
2005 for (Int_t row = 0; row < rowMax; ++row) {
2006 for (Int_t col = 0; col < colMax; ++col) {
2007 calROC->SetValue(col,row,TMath::Abs(values[(Int_t) (col*rowMax+row)]));
2016 //_____________________________________________________________________________
2017 void AliTRDCalibraFit::SetPeriodeFitPH(Int_t periodeFitPH)
2020 // Set FitPH if 1 then each detector will be fitted
2023 if (periodeFitPH > 0) {
2024 fFitPHPeriode = periodeFitPH;
2027 AliInfo("periodeFitPH must be higher than 0!");
2032 //_____________________________________________________________________________
2033 void AliTRDCalibraFit::SetFitPRFNDB(Int_t fitPRFNDB)
2036 // TO choose the method that you write into the database
2039 if ((fitPRFNDB >= 3) || (fitPRFNDB == 1)) {
2040 AliInfo("fitPRFNDB is not a correct number!");
2043 fFitPRFNDB = fitPRFNDB;
2048 //_____________________________________________________________________________
2049 void AliTRDCalibraFit::SetFitChargeNDB(Int_t fitChargeNDB)
2052 // To choose the method that you write into the database
2054 if ((fitChargeNDB >= 5) || (fitChargeNDB == 3)) {
2055 AliInfo("fitChargeNDB is not a correct number!");
2058 fFitChargeNDB = fitChargeNDB;
2063 //_____________________________________________________________________________
2064 void AliTRDCalibraFit::SetFitPHNDB(Int_t fitPHNDB)
2067 // To choose the method that you write into the database
2070 if ((fitPHNDB >= 4) || (fitPHNDB == 2)) {
2071 AliInfo("fitPHNDB is not a correct number!");
2074 fFitPHNDB = fitPHNDB;
2079 //_____________________________________________________________________________
2080 void AliTRDCalibraFit::SetBeginFitCharge(Float_t beginFitCharge)
2083 // The fit of the deposited charge distribution begins at
2084 // histo->Mean()/beginFitCharge
2085 // You can here set beginFitCharge
2088 if (beginFitCharge > 0) {
2089 fBeginFitCharge = beginFitCharge;
2092 AliInfo("beginFitCharge must be strict positif!");
2097 //_____________________________________________________________________________
2098 void AliTRDCalibraFit::SetT0Shift(Float_t t0Shift)
2101 // The t0 calculated with the maximum positif slope is shift from t0Shift
2102 // You can here set t0Shift
2109 AliInfo("t0Shift must be strict positif!");
2114 //_____________________________________________________________________________
2115 void AliTRDCalibraFit::SetRangeFitPRF(Float_t rangeFitPRF)
2118 // The fit of the PRF is from -rangeFitPRF to rangeFitPRF
2119 // You can here set rangeFitPRF
2122 if ((rangeFitPRF > 0) &&
2123 (rangeFitPRF <= 1.5)) {
2124 fRangeFitPRF = rangeFitPRF;
2127 AliInfo("rangeFitPRF must be between 0 and 1.0");
2132 //_____________________________________________________________________________
2133 void AliTRDCalibraFit::SetRebin(Short_t rebin)
2136 // Rebin with rebin time less bins the Ch histo
2137 // You can set here rebin that should divide the number of bins of CH histo
2142 AliInfo("You have to be sure that fRebin divides fNumberBinCharge used!");
2145 AliInfo("You have to choose a positiv value!");
2150 //____________Pad Calibration Public___________________________________________
2152 //____________Protected Functions______________________________________________
2153 //____________Create the 2D histo to be filled online__________________________
2155 //____________Fit______________________________________________________________
2156 //____________Create histos if fDebug == 1 or fDebug >= 3______________________
2158 //_____________________________________________________________________________
2159 void AliTRDCalibraFit::InitArrayFitPH()
2162 // Initialise fCoefVdrift[3] and fCoefVdriftE[2] to the right dimension
2165 Int_t nbins = fDect2[1]-fDect1[1];
2167 fCoefVdrift[2] = new Double_t[nbins];
2169 // Init the pointer to nbins
2171 fCoefVdrift[0] = new Double_t[nbins];
2172 fCoefVdriftE[0] = new Double_t[nbins];
2173 for(Int_t k = 0; k < nbins; k++){
2174 fCoefVdriftE[0][k] = 0.0;
2180 fCoefVdrift[1] = new Double_t[nbins];
2181 fCoefVdriftE[1] = new Double_t[nbins];
2182 for(Int_t k = 0; k < nbins; k++){
2183 fCoefVdriftE[1][k] = 0.0;
2187 fCoefVdrift[3] = new Double_t[nbins];
2188 fCoefVdriftE[2] = new Double_t[nbins];
2189 for(Int_t k = 0; k < nbins; k++){
2190 fCoefVdriftE[2][k] = 0.0;
2196 //_____________________________________________________________________________
2197 void AliTRDCalibraFit::InitArrayFitT0()
2200 // Initialise fCoefT0[3] and fCoefT0E[2] to the right dimension
2203 Int_t nbins = fDect2[1]-fDect1[1];
2205 fCoefT0[2] = new Double_t[nbins];
2207 // Init the pointer to nbins
2209 fCoefT0[0] = new Double_t[nbins];
2210 fCoefT0E[0] = new Double_t[nbins];
2211 for(Int_t k = 0; k < nbins; k++){
2212 fCoefT0E[0][k] = 0.0;
2216 fCoefT0[1] = new Double_t[nbins];
2217 fCoefT0E[1] = new Double_t[nbins];
2218 for(Int_t k = 0; k < nbins; k++){
2219 fCoefT0E[1][k] = 0.0;
2223 fCoefT0[3] = new Double_t[nbins];
2224 fCoefT0E[2] = new Double_t[nbins];
2225 for(Int_t k = 0; k < nbins; k++){
2226 fCoefT0E[2][k] = 0.0;
2232 //_____________________________________________________________________________
2233 void AliTRDCalibraFit::InitArrayFitCH()
2236 // Initialise fCoefCharge[4] and fCoefChargeE[3] to the right dimension
2239 Int_t nbins = fDect2[0]-fDect1[0];
2241 //Init the pointer to nbins
2243 fCoefCharge[1] = new Double_t[nbins];
2244 fCoefChargeE[1] = new Double_t[nbins];
2245 for(Int_t k = 0; k < nbins; k++){
2246 fCoefChargeE[1][k] = 0.0;
2250 fCoefCharge[4] = new Double_t[nbins];
2251 fCoefChargeE[3] = new Double_t[nbins];
2252 for(Int_t k = 0; k < nbins; k++){
2253 fCoefChargeE[3][k] = 0.0;
2257 fCoefCharge[0] = new Double_t[nbins];
2258 fCoefChargeE[0] = new Double_t[nbins];
2259 for(Int_t k = 0; k < nbins; k++){
2260 fCoefChargeE[0][k] = 0.0;
2264 if(fFitChargeBisOn){
2265 fCoefCharge[2] = new Double_t[nbins];
2266 fCoefChargeE[2] = new Double_t[nbins];
2267 for(Int_t k = 0; k < nbins; k++){
2268 fCoefChargeE[2][k] = 0.0;
2272 fCoefCharge[3] = new Double_t[nbins];
2276 //_____________________________________________________________________________
2277 void AliTRDCalibraFit::InitArrayFitPRF()
2280 // Initialise fCoefPRF[2] and fCoefPRFE to the right dimension
2283 Int_t nbins = fDect2[2]-fDect1[2];
2284 fCoefPRF[1] = new Double_t[nbins];
2286 //Init the pointer to nbins
2288 fCoefPRF[0] = new Double_t[nbins];
2289 fCoefPRFE[0] = new Double_t[nbins];
2290 for(Int_t k = 0; k < nbins; k++){
2291 fCoefPRFE[0][k] = 0.0;
2295 fCoefPRF[2] = new Double_t[nbins];
2296 fCoefPRFE[1] = new Double_t[nbins];
2297 for(Int_t k = 0; k < nbins; k++){
2298 fCoefPRFE[1][k] = 0.0;
2303 //_____________________________________________________________________________
2304 void AliTRDCalibraFit::CreateFitHistoPRFDB(Int_t rowMax, Int_t colMax)
2307 // Create the histos for fDebug = 3 and fDebug = 4 (Fit functions)
2310 fCoefPRFDB[0] = new TH2F("coefPRF0","",rowMax,0,rowMax,colMax,0,colMax);
2311 fCoefPRFDB[0]->SetStats(0);
2312 fCoefPRFDB[0]->SetXTitle("row Number");
2313 fCoefPRFDB[0]->SetYTitle("col Number");
2314 fCoefPRFDB[0]->SetZTitle("PRF width [pad width units]");
2315 fCoefPRFDB[0]->SetFillColor(6);
2316 fCoefPRFDB[0]->SetLineColor(6);
2319 fCoefPRFDB[1] = new TH2F("coefPRF1","",rowMax,0,rowMax,colMax,0,colMax);
2320 fCoefPRFDB[1]->SetStats(0);
2321 fCoefPRFDB[1]->SetXTitle("row Number");
2322 fCoefPRFDB[1]->SetYTitle("col Number");
2323 fCoefPRFDB[1]->SetZTitle("PRF width [pad width units]");
2324 fCoefPRFDB[1]->SetFillColor(1);
2325 fCoefPRFDB[1]->SetLineColor(1);
2329 //_____________________________________________________________________________
2330 void AliTRDCalibraFit::CreateFitHistoCHDB(Int_t rowMax, Int_t colMax)
2333 // Create the histos for fDebug = 3 and fDebug = 4 (Fit functions)
2337 fCoefChargeDB[0] = new TH2F("coefchargedb0","",rowMax,0,rowMax,colMax,0,colMax);
2338 fCoefChargeDB[0]->SetStats(0);
2339 fCoefChargeDB[0]->SetXTitle("row Number");
2340 fCoefChargeDB[0]->SetYTitle("col Number");
2341 fCoefChargeDB[0]->SetZTitle("f_{g} Fit method");
2342 fCoefChargeDB[0]->SetFillColor(6);
2343 fCoefChargeDB[0]->SetLineColor(6);
2345 if(fFitChargeBisOn){
2346 fCoefChargeDB[2] = new TH2F("coefchargedb2","",rowMax,0,rowMax,colMax,0,colMax);
2347 fCoefChargeDB[2]->SetStats(0);
2348 fCoefChargeDB[2]->SetXTitle("row Number");
2349 fCoefChargeDB[2]->SetYTitle("col Number");
2350 fCoefChargeDB[2]->SetZTitle("f_{g} Fitbis method");
2351 fCoefChargeDB[2]->SetFillColor(8);
2352 fCoefChargeDB[2]->SetLineColor(8);
2355 fCoefChargeDB[3] = new TH2F("coefchargedb3","",rowMax,0,rowMax,colMax,0,colMax);
2356 fCoefChargeDB[3]->SetStats(0);
2357 fCoefChargeDB[3]->SetXTitle("row Number");
2358 fCoefChargeDB[3]->SetYTitle("col Number");
2359 fCoefChargeDB[3]->SetFillColor(1);
2360 fCoefChargeDB[3]->SetLineColor(1);
2364 fCoefChargeDB[1] = new TH2F("coefchargedb1","",rowMax,0,rowMax,colMax,0,colMax);
2365 fCoefChargeDB[1]->SetStats(0);
2366 fCoefChargeDB[1]->SetXTitle("row Number");
2367 fCoefChargeDB[1]->SetYTitle("col Number");
2368 fCoefChargeDB[1]->SetZTitle("f_{g} Mean method");
2369 fCoefChargeDB[1]->SetFillColor(2);
2370 fCoefChargeDB[1]->SetLineColor(2);
2375 //_____________________________________________________________________________
2376 void AliTRDCalibraFit::CreateFitHistoPHDB(Int_t rowMax, Int_t colMax)
2379 // Create the histos for fDebug = 3 and fDebug = 4 (Fit functions)
2383 fCoefVdriftDB[0] = new TH2F("coefvdriftdb0","",rowMax,0,rowMax,colMax,0,colMax);
2384 fCoefVdriftDB[0]->SetStats(0);
2385 fCoefVdriftDB[0]->SetXTitle("row Number");
2386 fCoefVdriftDB[0]->SetYTitle("col Number");
2387 fCoefVdriftDB[0]->SetZTitle("v_{drift} Fit method");
2388 fCoefVdriftDB[0]->SetFillColor(6);
2389 fCoefVdriftDB[0]->SetLineColor(6);
2393 fCoefVdriftDB[1] = new TH2F("coefvdriftdb1","",rowMax,0,rowMax,colMax,0,colMax);
2394 fCoefVdriftDB[1]->SetStats(0);
2395 fCoefVdriftDB[1]->SetXTitle("row Number");
2396 fCoefVdriftDB[1]->SetYTitle("col Number");
2397 fCoefVdriftDB[1]->SetZTitle("v_{drift} slope method");
2398 fCoefVdriftDB[1]->SetFillColor(2);
2399 fCoefVdriftDB[1]->SetLineColor(2);
2402 fCoefVdriftDB[2] = new TH2F("coefvdriftdb1","",rowMax,0,rowMax,colMax,0,colMax);
2403 fCoefVdriftDB[2]->SetStats(0);
2404 fCoefVdriftDB[2]->SetXTitle("row Number");
2405 fCoefVdriftDB[2]->SetYTitle("col Number");
2406 fCoefVdriftDB[2]->SetZTitle("v_{drift} slope method");
2407 fCoefVdriftDB[2]->SetFillColor(1);
2408 fCoefVdriftDB[2]->SetLineColor(1);
2413 //_____________________________________________________________________________
2414 void AliTRDCalibraFit::CreateFitHistoT0DB(Int_t rowMax, Int_t colMax)
2417 // Create the histos for fDebug = 3 and fDebug = 4 (Fit functions)
2421 fCoefT0DB[0] = new TH2F("coefT0db0","",rowMax,0,rowMax,colMax,0,colMax);
2422 fCoefT0DB[0]->SetStats(0);
2423 fCoefT0DB[0]->SetXTitle("row Number");
2424 fCoefT0DB[0]->SetYTitle("col Number");
2425 fCoefT0DB[0]->SetZTitle("t0 Fit method");
2426 fCoefT0DB[0]->SetFillColor(6);
2427 fCoefT0DB[0]->SetLineColor(6);
2430 fCoefT0DB[1] = new TH2F("coefT0db1","",rowMax,0,rowMax,colMax,0,colMax);
2431 fCoefT0DB[1]->SetStats(0);
2432 fCoefT0DB[1]->SetXTitle("row Number");
2433 fCoefT0DB[1]->SetYTitle("col Number");
2434 fCoefT0DB[1]->SetZTitle("t0 slope method");
2435 fCoefT0DB[1]->SetFillColor(2);
2436 fCoefT0DB[1]->SetLineColor(2);
2439 fCoefT0DB[2] = new TH2F("coefT0db1","",rowMax,0,rowMax,colMax,0,colMax);
2440 fCoefT0DB[2]->SetStats(0);
2441 fCoefT0DB[2]->SetXTitle("row Number");
2442 fCoefT0DB[2]->SetYTitle("col Number");
2443 fCoefT0DB[2]->SetZTitle("t0 slope method");
2444 fCoefT0DB[2]->SetFillColor(1);
2445 fCoefT0DB[2]->SetLineColor(1);
2450 //_____________________________________________________________________________
2451 Bool_t AliTRDCalibraFit::FillVectorFitCH(Int_t countdet)
2454 // For the Fit functions fill the vector FitCH special for the gain calibration
2457 AliTRDFitCHInfo *fitCHInfo = new AliTRDFitCHInfo();
2460 if (GetChamber(countdet) == 2) {
2467 //printf("For the detector %d , ntotal %d and fCoefCH[0] %f\n",countdet,ntotal,fCoefCH[0]);
2468 Float_t *coef = new Float_t[ntotal];
2469 for (Int_t i = 0; i < ntotal; i++) {
2470 coef[i] = fCoefCH[i];
2473 Int_t detector = countdet;
2475 fitCHInfo->SetCoef(coef);
2476 fitCHInfo->SetDetector(detector);
2477 fVectorFitCH->Add((TObject *) fitCHInfo);
2483 //____________Functions for initialising the AliTRDCalibraFit in the code_________
2484 Bool_t AliTRDCalibraFit::InitFit(Int_t nbins, Int_t i)
2487 // Init the calibration mode (Nz, Nrphi), the histograms for
2488 // debugging the fit methods if fDebug > 0,
2491 gStyle->SetPalette(1);
2492 gStyle->SetOptStat(1111);
2493 gStyle->SetPadBorderMode(0);
2494 gStyle->SetCanvasColor(10);
2495 gStyle->SetPadLeftMargin(0.13);
2496 gStyle->SetPadRightMargin(0.01);
2498 // Get the parameter object
2499 AliTRDCommonParam *parCom = AliTRDCommonParam::Instance();
2501 AliInfo("Could not get CommonParam");
2505 // Mode groups of pads: the total number of bins!
2506 Int_t numberofbinsexpected = 0;
2507 fCalibraMode->ModePadCalibration(2,i);
2508 fCalibraMode->ModePadFragmentation(0,2,0,i);
2509 fCalibraMode->SetDetChamb2(i);
2511 AliInfo(Form("For the chamber 2: %d",fCalibraMode->GetDetChamb2(i)));
2513 numberofbinsexpected += 6 * 18 * fCalibraMode->GetDetChamb2(i);
2514 fCalibraMode->ModePadCalibration(0,i);
2515 fCalibraMode->ModePadFragmentation(0,0,0,i);
2516 fCalibraMode->SetDetChamb0(i);
2518 AliInfo(Form("For the other chamber 0: %d",fCalibraMode->GetDetChamb0(i)));
2520 numberofbinsexpected += 6 * 4 * 18 * fCalibraMode->GetDetChamb0(i);
2522 // Quick verification that we have the good pad calibration mode if 2D histos!
2524 if (numberofbinsexpected != nbins) {
2525 AliInfo("It doesn't correspond to the mode of pad group calibration!");
2530 // Security for fDebug 3 and 4
2531 if ((fDebug >= 3) &&
2535 AliInfo("This detector doesn't exit!");
2539 // Determine fDet1 and fDet2
2543 fDect1[i] = fFitVoir;
2544 fDect2[i] = fDect1[i] +1;
2548 fDect2[i] = numberofbinsexpected;
2551 fCalibraMode->CalculXBins(AliTRDgeometry::GetDetector(fDet[0],fDet[1],fDet[2]),i);
2552 fDect1[i] = fCalibraMode->GetXbins(i);
2553 fCalibraMode->CalculXBins((AliTRDgeometry::GetDetector(fDet[0],fDet[1],fDet[2])+1),i);
2554 fDect2[i] = fCalibraMode->GetXbins(i);
2557 // Create the histos for debugging
2562 // Init the VectorFitCH
2563 fVectorFitCH = new TObjArray();
2564 fCoefCH = new Float_t[2304];
2565 for (Int_t k = 0; k < 2304; k++) {
2568 fScaleFitFactor = 0.0;
2570 // Number of Xbins(detectors or groups of pads) if Vector2d
2571 // Quick verification that we are not out of range!
2572 if (fCalibraVector) {
2574 (fCalibraVector->GetVectorCH()->GetEntriesFast() > 0) &&
2575 ((Int_t) fCalibraVector->GetPlaCH()->GetEntriesFast() > 0)) {
2576 if ((Int_t) fCalibraVector->GetVectorCH()->GetEntriesFast() > numberofbinsexpected) {
2577 AliInfo("ch doesn't correspond to the mode of pad group calibration!");
2580 if ((Int_t) fCalibraVector->GetVectorCH()->GetEntriesFast() !=
2581 (Int_t) fCalibraVector->GetPlaCH()->GetEntriesFast()) {
2582 AliInfo("VectorCH doesn't correspond to PlaCH!");
2589 // Debugging: Create the histos
2592 // fDebug == 0 nothing
2599 // fDebug == 2 and fFitVoir no histo
2601 if (fFitVoir < numberofbinsexpected) {
2602 AliInfo(Form("We will see the fit of the object %d",fFitVoir));
2605 AliInfo("fFitVoir is out of range of the histo!");
2610 // fDebug == 3 or 4 and fDet
2612 if ((fCalibraMode->GetNz(0) == 0) && (fCalibraMode->GetNrphi(0) == 0)) {
2613 AliInfo("Do you really want to see one detector without pad groups?");
2617 AliInfo(Form("You will see the detector: iPlane %d, iChamb %d, iSect %d"
2618 ,fDet[0],fDet[1],fDet[2]));
2619 // A little geometry:
2620 Int_t rowMax = parCom->GetRowMax(fDet[0],fDet[1],fDet[2]);
2621 Int_t colMax = parCom->GetColMax(fDet[0]);
2622 // Create the histos to visualise
2623 CreateFitHistoCHDB(rowMax,colMax);
2635 // Number of Xbins (detectors or groups of pads) if vector2d
2636 // Quick verification that we are not out of range!
2637 if (fCalibraVector) {
2639 (fCalibraVector->GetVectorPH()->GetEntriesFast() > 0) &&
2640 ((Int_t) fCalibraVector->GetPlaPH()->GetEntriesFast() > 0)) {
2641 if ((Int_t) fCalibraVector->GetVectorPH()->GetEntriesFast() > numberofbinsexpected) {
2642 AliInfo("ph doesn't correspond to the mode of pad group calibration!");
2645 if ((Int_t) fCalibraVector->GetVectorPH()->GetEntriesFast() !=
2646 (Int_t) fCalibraVector->GetPlaPH()->GetEntriesFast()) {
2647 AliInfo("VectorPH doesn't correspond to PlaPH!");
2658 // Debugging: Create the histos
2661 // fDebug == 0 nothing
2665 // Create the histos replique de ph
2670 // fDebug == 2 and fFitVoir no histo
2672 if (fFitVoir < numberofbinsexpected) {
2673 AliInfo(Form("We will see the fit of the object %d",fFitVoir));
2676 AliInfo("fFitVoir is out of range of the histo!");
2681 // fDebug == 3 or 4 and fDet
2683 if ((fCalibraMode->GetNz(1) == 0) &&
2684 (fCalibraMode->GetNrphi(1) == 0)) {
2685 AliInfo("Do you really want to see one detector without pad groups?");
2689 AliInfo(Form("You will see the detector: iPlane %d, iChamb %d, iSect %d"
2690 ,fDet[0],fDet[1],fDet[2]));
2691 // A little geometry:
2692 Int_t rowMax = parCom->GetRowMax(fDet[0],fDet[1],fDet[2]);
2693 Int_t colMax = parCom->GetColMax(fDet[0]);
2694 // Create the histos to visualise
2695 CreateFitHistoPHDB(rowMax,colMax);
2696 CreateFitHistoT0DB(rowMax,colMax);
2709 // Number of Xbins(detectors or groups of pads) if vector2d
2710 if (fCalibraVector){
2712 (fCalibraVector->GetVectorPRF()->GetEntriesFast() > 0) &&
2713 (fCalibraVector->GetPlaPRF()->GetEntriesFast() > 0)) {
2714 // Quick verification that we are not out of range!
2715 if ((Int_t) fCalibraVector->GetVectorPRF()->GetEntriesFast() > numberofbinsexpected) {
2716 AliInfo("ch doesn't correspond to the mode of pad group calibration!");
2719 if ((Int_t) fCalibraVector->GetVectorPRF()->GetEntriesFast() !=
2720 (Int_t) fCalibraVector->GetPlaPRF()->GetEntriesFast()) {
2721 AliInfo("VectorPRF doesn't correspond to PlaCH!");
2731 // Debugging: Create the histos
2734 // fDebug == 0 nothing
2738 // Create the histos replique de ch
2742 // fDebug == 2 and fFitVoir no histo
2744 if (fFitVoir < numberofbinsexpected) {
2745 AliInfo(Form("We will see the fit of the object %d",fFitVoir));
2748 AliInfo("fFitVoir is out of range of the histo!");
2753 // fDebug == 3 or 4 and fDet
2755 if ((fCalibraMode->GetNz(2) == 0) &&
2756 (fCalibraMode->GetNrphi(2) == 0)) {
2757 AliInfo("Do you really want to see one detector without pad groups?");
2761 AliInfo(Form("You will see the detector: iPlane %d, iChamb %d, iSect %d"
2762 ,fDet[0],fDet[1],fDet[2]));
2763 // A little geometry:
2764 Int_t rowMax = parCom->GetRowMax(fDet[0],fDet[1],fDet[2]);
2765 Int_t colMax = parCom->GetColMax(fDet[0]);
2766 // Create the histos to visualise
2767 CreateFitHistoPRFDB(rowMax,colMax);
2780 //____________Functions for initialising the AliTRDCalibraFit in the code_________
2781 void AliTRDCalibraFit::InitfCountDetAndfCount(Int_t i)
2784 // Init the current detector where we are fCountDet and the
2785 // next fCount for the functions Fit...
2788 // Loop on the Xbins of ch!!
2789 fCountDet[i] = -1; // Current detector
2790 fCount[i] = 0; // To find the next detector
2795 // Set countdet to the detector
2796 fCountDet[i] = AliTRDgeometry::GetDetector(fDet[0],fDet[1],fDet[2]);
2798 // Determination of fNnZ, fNnRphi, fNfragZ and fNfragRphi
2799 fCalibraMode->ModePadCalibration(fDet[1],i);
2800 fCalibraMode->ModePadFragmentation(fDet[0],fDet[1],fDet[2],i);
2802 // Set counter to write at the end of the detector
2803 fCount[i] = fDect1[i] + fCalibraMode->GetNfragZ(i)*fCalibraMode->GetNfragRphi(i);
2809 //____________Functions for initialising the AliTRDCalibraFit in the code_________
2810 void AliTRDCalibraFit::UpdatefCountDetAndfCount(Int_t idect, Int_t i)
2813 // See if we are in a new detector and update the
2814 // variables fNfragZ and fNfragRphi if yes
2817 // Determination of fNnZ, fNnRphi, fNfragZ and fNfragRphi
2818 // If fDebug == 1 or 0
2819 if ((fDebug == 0) ||
2822 if (fCount[i] == idect) {
2824 // On en est au detector
2827 // Determination of fNnZ, fNnRphi, fNfragZ and fNfragRphi
2828 fCalibraMode->ModePadCalibration((Int_t) GetChamber(fCountDet[i]),i);
2829 fCalibraMode->ModePadFragmentation((Int_t) GetPlane(fCountDet[i])
2830 ,(Int_t) GetChamber(fCountDet[i])
2831 ,(Int_t) GetSector(fCountDet[i]),i);
2833 // Set for the next detector
2834 fCount[i] += fCalibraMode->GetNfragZ(i)*fCalibraMode->GetNfragRphi(i);
2842 //____________Functions for initialising the AliTRDCalibraFit in the code_________
2843 void AliTRDCalibraFit::ReconstructFitRowMinRowMax(Int_t idect, Int_t i)
2846 // Reconstruct the min pad row, max pad row, min pad col and
2847 // max pad col of the calibration group for the Fit functions
2851 fCalibraMode->ReconstructionRowPadGroup((Int_t) (idect-(fCount[i]-(fCalibraMode->GetNfragZ(i)
2852 *fCalibraMode->GetNfragRphi(i)))),i);
2855 fCalibraMode->ReconstructionRowPadGroup((Int_t) (idect-fDect1[i]),i);
2860 //____________Functions for initialising the AliTRDCalibraFit in the code_________
2861 Bool_t AliTRDCalibraFit::NotEnoughStatistic(Int_t idect, Int_t i)
2864 // For the case where there are not enough entries in the histograms
2865 // of the calibration group, the value present in the choosen database
2866 // will be put. A negativ sign enables to know that a fit was not possible.
2869 // Get the parameter object
2870 AliTRDCommonParam *parCom = AliTRDCommonParam::Instance();
2872 AliInfo("Could not get CommonParam Manager");
2877 AliTRDcalibDB *cal = AliTRDcalibDB::Instance();
2879 AliInfo("Could not get calibDB");
2884 AliInfo(Form("The element %d in this detector %d has not enough statistic to be fitted"
2885 ,idect-(fCount[i]-(fCalibraMode->GetNfragZ(i)*fCalibraMode->GetNfragRphi(i))),fCountDet[i]));
2888 AliInfo("The element has not enough statistic to be fitted");
2891 if ((i == 0) && (fDebug != 2)) {
2893 // Calcul the coef from the database choosen
2894 CalculChargeCoefMean(fCountDet[0],(Int_t) (idect-fDect1[0]),kFALSE);
2896 // Fill the coefCH[2304] with negative value to say: not fitted
2897 AliInfo(Form("The row min %d, the row max %d, the colmin %d and the col max %d"
2898 ,fCalibraMode->GetRowMin(0)
2899 ,fCalibraMode->GetRowMax(0)
2900 ,fCalibraMode->GetColMin(0)
2901 ,fCalibraMode->GetColMax(0)));
2902 for (Int_t k = fCalibraMode->GetRowMin(0); k < fCalibraMode->GetRowMax(0); k++) {
2903 for (Int_t j = fCalibraMode->GetColMin(0); j < fCalibraMode->GetColMax(0); j++) {
2904 if (GetChamber(fCountDet[0]) == 2) {
2905 fCoefCH[(Int_t)(j*12+k)] = -TMath::Abs(fChargeCoef[3]);
2907 if (GetChamber(fCountDet[0]) != 2) {
2908 fCoefCH[(Int_t)(j*16+k)] = -TMath::Abs(fChargeCoef[3]);
2913 // Put the default value negative
2914 if ((fDebug == 1) ||
2917 if (fFitChargeBisOn) {
2918 fCoefCharge[2][idect-fDect1[0]]=-TMath::Abs(fChargeCoef[3]);
2920 if (fMeanChargeOn) {
2921 fCoefCharge[1][idect-fDect1[0]]=-TMath::Abs(fChargeCoef[3]);
2924 fCoefCharge[0][idect-fDect1[0]]=-TMath::Abs(fChargeCoef[3]);
2929 // End of one detector
2930 if ((idect == (fCount[0]-1))) {
2931 FillVectorFitCH((Int_t) fCountDet[0]);
2933 for (Int_t k = 0; k < 2304; k++) {
2940 if ((i == 1) && (fDebug != 2)) {
2942 CalculVdriftCoefMean(fCountDet[1],(Int_t) (idect-fDect1[1]));
2943 CalculT0CoefMean(fCountDet[1],(Int_t) (idect-fDect1[1]));
2945 // Put the default value (time0 can be negativ, so we stay with + )
2946 if ((fDebug == 1) ||
2950 fCoefVdrift[0][(idect-fDect1[1])] = -fVdriftCoef[2];
2951 fCoefT0[0][(idect-fDect1[1])] = fT0Coef[2];
2955 fCoefVdrift[1][(idect-fDect1[1])] = -fVdriftCoef[2];
2956 fCoefT0[1][(idect-fDect1[1])] = fT0Coef[2];
2959 fCoefVdrift[3][(idect-fDect1[1])] = -fVdriftCoef[2];
2960 fCoefT0[3][(idect-fDect1[1])] = fT0Coef[2];
2965 // Put the default value
2968 fVdriftCoef[0] = fVdriftCoef[2];
2969 fT0Coef[0] = fT0Coef[2];
2972 fVdriftCoef[1] = fVdriftCoef[2];
2973 fT0Coef[1] = fT0Coef[2];
2976 fVdriftCoef[3] = fVdriftCoef[2];
2977 fT0Coef[3] = fT0Coef[2];
2983 // Fill the tree if end of a detector.
2984 // The pointer to the branch stays with the default value negative!!!
2986 // Pointer to the branch
2987 for (Int_t k = fCalibraMode->GetRowMin(1); k < fCalibraMode->GetRowMax(1); k++) {
2988 for (Int_t j = fCalibraMode->GetColMin(1); j < fCalibraMode->GetColMax(1); j++) {
2989 if (GetChamber(fCountDet[1]) == 2) {
2990 fVdriftPad[(Int_t)(j*12+k)] = -TMath::Abs(fVdriftCoef[2]);
2992 if (GetChamber(fCountDet[1]) != 2) {
2993 fVdriftPad[(Int_t)(j*16+k)] = -TMath::Abs(fVdriftCoef[2]);
2998 // End of one detector
2999 if ((idect == (fCount[1]-1)) && (fDebug != 2)) {
3000 FillTreeVdrift((Int_t) fCountDet[1]);
3004 // Fill the tree if end of a detector.
3005 // The pointer to the branch stays with the default value positive!!!
3006 // Pointer to the branch
3007 for (Int_t k = fCalibraMode->GetRowMin(1); k < fCalibraMode->GetRowMax(1); k++) {
3008 for (Int_t j = fCalibraMode->GetColMin(1); j < fCalibraMode->GetColMax(1); j++) {
3009 if (GetChamber(fCountDet[1]) == 2) {
3010 fT0Pad[(Int_t)(j*12+k)] = fT0Coef[2];
3012 if (GetChamber(fCountDet[1]) != 2) {
3013 fT0Pad[(Int_t)(j*16+k)] = fT0Coef[2];
3018 // End of one detector
3019 if ((idect == (fCount[1]-1)) && (fDebug != 2)) {
3020 FillTreeT0((Int_t) fCountDet[1]);
3025 if ((i == 2) && (fDebug != 2)) {
3027 CalculPRFCoefMean(fCountDet[2],(Int_t) (idect-fDect1[2]));
3029 if ((fDebug == 1) ||
3032 fCoefPRF[0][(idect-fDect1[2])] = -fPRFCoef[1];
3035 fCoefPRF[2][(idect-fDect1[2])] = -fPRFCoef[1];
3041 fPRFCoef[0] = fPRFCoef[1];
3044 fPRFCoef[2] = fPRFCoef[1];
3049 // Fill the tree if end of a detector.
3050 // The pointer to the branch stays with the default value 1.5!!!
3051 // Pointer to the branch
3052 for (Int_t k = fCalibraMode->GetRowMin(2); k < fCalibraMode->GetRowMax(2); k++) {
3053 for (Int_t j = fCalibraMode->GetColMin(2); j < fCalibraMode->GetColMax(2); j++) {
3054 if((parCom->GetColMax(GetPlane(fCountDet[2])) != (j+1)) && (j != 0)){
3055 if (GetChamber(fCountDet[2]) == 2) {
3056 fPRFPad[(Int_t)(j*12+k)] = -fPRFCoef[1];
3058 if (GetChamber(fCountDet[2]) != 2) {
3059 fPRFPad[(Int_t)(j*16+k)] = -fPRFCoef[1];
3064 if (GetChamber(fCountDet[2]) == 2) {
3065 fPRFPad[(Int_t)(j*12+k)] = -((Float_t) cal->GetPRFWidth(fCountDet[2],j,k));
3067 if (GetChamber(fCountDet[2]) != 2) {
3068 fPRFPad[(Int_t)(j*16+k)] = -((Float_t) cal->GetPRFWidth(fCountDet[2],j,k));
3072 if (GetChamber(fCountDet[2]) == 2) {
3073 fPRFPad[(Int_t)(j*12+k)] = -((Float_t) GetPRFDefault(GetPlane(fCountDet[2])));
3075 if (GetChamber(fCountDet[2]) != 2) {
3076 fPRFPad[(Int_t)(j*16+k)] = -((Float_t) GetPRFDefault(GetPlane(fCountDet[2])));
3083 // End of one detector
3084 if ((idect == (fCount[2]-1)) && (fDebug != 2)) {
3085 FillTreePRF((Int_t) fCountDet[2]);
3094 //____________Functions for initialising the AliTRDCalibraFit in the code_________
3095 Bool_t AliTRDCalibraFit::FillInfosFit(Int_t idect, Int_t i)
3098 // Fill the coefficients found with the fits or other
3099 // methods from the Fit functions
3102 // Get the parameter object
3103 AliTRDCommonParam *parCom = AliTRDCommonParam::Instance();
3105 AliInfo("Could not get CommonParam Manager");
3110 AliTRDcalibDB *cal = AliTRDcalibDB::Instance();
3112 AliInfo("Could not get calibDB");
3116 if ((i == 0) && (fDebug != 2)) {
3117 // Fill the coefCH[2304] with fChargeCoef[0]
3118 // that would be negativ only if the fit failed totally
3119 //printf("for fCountDet %d we have %f\n",fCountDet[0],fChargeCoef[fFitChargeNDB]);
3120 //printf("RowMin %d RowMax %d ColMin %d ColMax %d\n",fCalibraMode->GetRowMin(0),fCalibraMode->GetRowMax(0),fCalibraMode->GetColMin(0),fCalibraMode->GetColMax(0));
3121 for (Int_t k = fCalibraMode->GetRowMin(0); k < fCalibraMode->GetRowMax(0); k++) {
3122 for (Int_t j = fCalibraMode->GetColMin(0); j < fCalibraMode->GetColMax(0); j++) {
3123 if (GetChamber(fCountDet[0]) == 2) {
3124 fCoefCH[(Int_t)(j*12+k)] = fChargeCoef[fFitChargeNDB];
3126 if (GetChamber(fCountDet[0]) != 2) {
3127 fCoefCH[(Int_t)(j*16+k)] = fChargeCoef[fFitChargeNDB];
3131 // End of one detector
3132 if ((idect == (fCount[0]-1))) {
3133 FillVectorFitCH((Int_t) fCountDet[0]);
3135 for (Int_t k = 0; k < 2304; k++) {
3141 if ((i == 1) && (fDebug != 2)) {
3144 // Pointer to the branch: fVdriftCoef[1] will ne negativ only if the fit failed totally
3145 for (Int_t k = fCalibraMode->GetRowMin(1); k < fCalibraMode->GetRowMax(1); k++) {
3146 for (Int_t j = fCalibraMode->GetColMin(1); j < fCalibraMode->GetColMax(1); j++) {
3147 if (GetChamber(fCountDet[1]) == 2) {
3148 fVdriftPad[(Int_t)(j*12+k)]=fVdriftCoef[fFitPHNDB];
3150 if (GetChamber(fCountDet[1]) != 2) {
3151 fVdriftPad[(Int_t)(j*16+k)]=fVdriftCoef[fFitPHNDB];
3155 // End of one detector
3156 if ((idect == (fCount[1]-1)) && (fDebug != 2)) {
3157 FillTreeVdrift((Int_t) fCountDet[1]);
3161 // Pointer to the branch: fT0Coef[1] will ne negativ only if the fit failed totally
3162 for (Int_t k = fCalibraMode->GetRowMin(1); k < fCalibraMode->GetRowMax(1); k++) {
3163 for (Int_t j = fCalibraMode->GetColMin(1); j < fCalibraMode->GetColMax(1); j++) {
3164 if (GetChamber(fCountDet[1]) == 2) {
3165 fT0Pad[(Int_t)(j*12+k)]=fT0Coef[fFitPHNDB];
3167 if (GetChamber(fCountDet[1]) != 2) {
3168 fT0Pad[(Int_t)(j*16+k)]=fT0Coef[fFitPHNDB];
3172 // End of one detector
3173 if ((idect == (fCount[1]-1)) && (fDebug != 2)) {
3174 FillTreeT0((Int_t) fCountDet[1]);
3179 if ((i == 2) && (fDebug != 2)) {
3180 // Pointer to the branch
3181 for (Int_t k = fCalibraMode->GetRowMin(2); k < fCalibraMode->GetRowMax(2); k++) {
3182 for (Int_t j = fCalibraMode->GetColMin(2); j < fCalibraMode->GetColMax(2); j++) {
3183 if ((parCom->GetColMax(GetPlane(fCountDet[2])) != (j+1)) && (j != 0)) {
3184 if (GetChamber(fCountDet[2]) == 2) {
3185 fPRFPad[(Int_t)(j*12+k)] = fPRFCoef[fFitPRFNDB];
3187 if (GetChamber(fCountDet[2]) != 2) {
3188 fPRFPad[(Int_t)(j*16+k)] = fPRFCoef[fFitPRFNDB];
3193 if (GetChamber(fCountDet[2]) == 2) {
3194 fPRFPad[(Int_t)(j*12+k)] = (Float_t) cal->GetPRFWidth(fCountDet[2],j,k);
3196 if (GetChamber(fCountDet[2]) != 2) {
3197 fPRFPad[(Int_t)(j*16+k)] = (Float_t) cal->GetPRFWidth(fCountDet[2],j,k);
3201 if (GetChamber(fCountDet[2]) == 2) {
3202 fPRFPad[(Int_t)(j*12+k)] = (Float_t) GetPRFDefault(GetPlane(fCountDet[2]));
3204 if (GetChamber(fCountDet[2]) != 2) {
3205 fPRFPad[(Int_t)(j*16+k)] = (Float_t) GetPRFDefault(GetPlane(fCountDet[2]));
3211 // End of one detector
3212 if ((idect == (fCount[2]-1)) && (fDebug != 2)) {
3213 FillTreePRF((Int_t) fCountDet[2]);
3221 //____________Functions for initialising the AliTRDCalibraFit in the code_________
3222 Bool_t AliTRDCalibraFit::WriteFitInfos(Int_t i)
3225 // In the case the user wants to write a file with a tree of the found
3226 // coefficients for the calibration before putting them in the database
3229 TFile *fout = TFile::Open(fWriteNameCoef,"UPDATE");
3230 // Check if the file could be opened
3231 if (!fout || !fout->IsOpen()) {
3232 AliInfo("No File found!");
3236 if ((i == 0) && (fDebug != 2)) {
3238 if ((fDebug == 4) ||
3243 fout->WriteTObject(fGain,fGain->GetName(),(Option_t *) "writedelete");
3246 if ((i == 1) && (fDebug != 2)) {
3248 if ((fDebug == 4) ||
3253 fout->WriteTObject(fVdrift,fVdrift->GetName(),(Option_t *) "writedelete");
3255 if ((fDebug == 4) ||
3260 fout->WriteTObject(fT0,fT0->GetName(),(Option_t *) "writedelete");
3263 if ((i == 2) && (fDebug != 2)) {
3265 if ((fDebug == 4) ||
3270 fout->WriteTObject(fPRF,fPRF->GetName(),(Option_t *) "writedelete");
3280 //____________Fill Coef DB in case of visualisation of one detector____________
3283 //_____________________________________________________________________________
3284 void AliTRDCalibraFit::FillCoefVdriftDB()
3287 // Fill the histos for fDebug = 3 and fDebug = 4 to visualise the detector
3290 for (Int_t row = fCalibraMode->GetRowMin(1); row < fCalibraMode->GetRowMax(1); row++) {
3291 for (Int_t col = fCalibraMode->GetColMin(1); col < fCalibraMode->GetColMax(1); col++) {
3293 fCoefVdriftDB[1]->SetBinContent(row+1,col+1,TMath::Abs(fVdriftCoef[1]));
3296 fCoefVdriftDB[0]->SetBinContent(row+1,col+1,TMath::Abs(fVdriftCoef[0]));
3298 if (fFitLagrPolOn ) {
3299 fCoefVdriftDB[2]->SetBinContent(row+1,col+1,TMath::Abs(fVdriftCoef[3]));
3306 //_____________________________________________________________________________
3307 void AliTRDCalibraFit::FillCoefT0DB()
3310 // Fill the histos for fDebug = 3 and fDebug = 4 to visualise the detector
3313 for (Int_t row = fCalibraMode->GetRowMin(1); row < fCalibraMode->GetRowMax(1); row++) {
3314 for (Int_t col = fCalibraMode->GetColMin(1); col < fCalibraMode->GetColMax(1); col++) {
3316 fCoefT0DB[1]->SetBinContent(row+1,col+1,TMath::Abs(fT0Coef[1]));
3319 fCoefT0DB[0]->SetBinContent(row+1,col+1,TMath::Abs(fT0Coef[0]));
3321 if (fFitLagrPolOn) {
3322 fCoefT0DB[2]->SetBinContent(row+1,col+1,TMath::Abs(fT0Coef[3]));
3329 //_____________________________________________________________________________
3330 void AliTRDCalibraFit::FillCoefChargeDB()
3333 // Fill the histos for fDebug = 3 and fDebug = 4 to visualise the detector
3336 for (Int_t row = fCalibraMode->GetRowMin(0); row < fCalibraMode->GetRowMax(0); row++) {
3337 for (Int_t col = fCalibraMode->GetColMin(0); col < fCalibraMode->GetColMax(0); col++) {
3338 if (fMeanChargeOn) {
3339 fCoefChargeDB[1]->SetBinContent(row+1,col+1,TMath::Abs(fChargeCoef[1]));
3341 if (fFitChargeBisOn) {
3342 fCoefChargeDB[2]->SetBinContent(row+1,col+1,TMath::Abs(fChargeCoef[2]));
3345 fCoefChargeDB[0]->SetBinContent(row+1,col+1,TMath::Abs(fChargeCoef[0]));
3348 fCoefChargeDB[3]->SetBinContent(row+1,col+1,TMath::Abs(fChargeCoef[4]));
3355 //_____________________________________________________________________________
3356 void AliTRDCalibraFit::FillCoefPRFDB()
3359 // Fill the histos for fDebug = 3 and fDebug = 4 to visualise the detector
3362 for (Int_t row = fCalibraMode->GetRowMin(2); row < fCalibraMode->GetRowMax(2); row++) {
3363 for (Int_t col = fCalibraMode->GetColMin(2); col < fCalibraMode->GetColMax(2); col++) {
3364 fCoefPRFDB[0]->SetBinContent(row+1,col+1,fPRFCoef[0]);
3369 for (Int_t row = fCalibraMode->GetRowMin(2); row < fCalibraMode->GetRowMax(2); row++) {
3370 for (Int_t col = fCalibraMode->GetColMin(2); col < fCalibraMode->GetColMax(2); col++) {
3371 fCoefPRFDB[1]->SetBinContent(row+1,col+1,fPRFCoef[2]);
3379 //____________Plot histos CoefPRF....__________________________________________
3382 //_____________________________________________________________________________
3383 void AliTRDCalibraFit::PlotWriteCH()
3386 // Scale the coefficients to one, create the graph errors and write them if wanted
3389 //TObjArray of the grapherrors and so on
3390 TObjArray *listofgraphs = new TObjArray();
3392 Int_t nbins = fDect2[0]-fDect1[0];
3395 // We will check fScaleFitFactor for the fFitChargeNDB, otherwise we calculate and normalise to 1
3396 // It can be that fScaleFitFactor is different from scale if we have taken a no default database as reference
3401 counter[0] = 0; //how many groups are fitted for 0
3402 counter[1] = 0; //how many groups are with mean for 1
3403 counter[2] = 0; //how many groups are fitted for 2
3404 counter[3] = 0; //how many groups are fitted for 4
3406 Double_t scale = 1.0;
3409 // Is -1 if no fit or mean, is 1 if fit or mean
3410 Double_t *xValuesFitted = new Double_t[nbins];
3411 Double_t *xValuesFittedMean = new Double_t[nbins];
3412 Double_t *xValuesFittedBis = new Double_t[nbins];
3413 Double_t *xValuesFittedMeanW = new Double_t[nbins];
3414 for(Int_t k = 0; k < nbins; k ++){
3415 xValuesFitted[k] = -1;
3416 xValuesFittedMean[k] = -1;
3417 xValuesFittedMeanW[k] = -1;
3418 xValuesFittedBis[k] = -1;
3423 for(Int_t l = 0; l < nbins; l++){
3424 if(fCoefCharge[0][l] > 0){
3425 sum += fCoefCharge[0][l];
3426 xValuesFitted[counter[0]]= l;
3431 if(sum > 0.0) scale = counter[0]/sum;
3432 if(fFitChargeNDB == 0){
3433 if(scale != fScaleFitFactor){
3434 AliInfo(Form("The normalisation is different from a nomalisation to one."));
3435 AliInfo(Form("For one we have %f and here %f",scale,fScaleFitFactor));
3436 if(!fAccCDB) AliInfo(Form("It is not normal because we didn't choose a reference database!"))
3438 scale = fScaleFitFactor;
3440 for(Int_t l = 0; l < nbins; l++){
3441 if(fCoefCharge[0][l] > 0){
3442 fCoefCharge[0][l]=fCoefCharge[0][l]*scale;
3443 fCoefChargeE[0][l]=fCoefChargeE[0][l]*scale;
3450 for(Int_t l = 0; l < nbins; l++){
3451 if(fCoefCharge[4][l] > 0){
3452 sum += fCoefCharge[4][l];
3453 xValuesFittedMeanW[counter[3]]= l;
3458 if(sum > 0.0) scale = counter[3]/sum;
3459 if(fFitChargeNDB == 4){
3460 if(scale != fScaleFitFactor){
3461 AliInfo(Form("The normalisation is different from a nomalisation to one."));
3462 AliInfo(Form("For one we have %f and here %f",scale,fScaleFitFactor));
3463 if(!fAccCDB) AliInfo(Form("It is not normal because we didn't choose a reference database!"))
3465 scale = fScaleFitFactor;
3467 for(Int_t l = 0; l < nbins; l++){
3468 if(fCoefCharge[4][l] > 0){
3469 fCoefCharge[4][l]=fCoefCharge[4][l]*scale;
3470 fCoefChargeE[3][l]=fCoefChargeE[3][l]*scale;
3477 for(Int_t l = 0; l < nbins; l++){
3478 if(fCoefCharge[1][l] > 0){
3479 sum += fCoefCharge[1][l];
3480 xValuesFittedMean[counter[1]]= l;
3485 if(sum > 0.0) scale = counter[1]/sum;
3486 if(fFitChargeNDB == 1){
3487 if(scale != fScaleFitFactor){
3488 AliInfo(Form("The normalisation is different from a nomalisation to one."));
3489 AliInfo(Form("For one we have %f and here %f",scale,fScaleFitFactor));
3490 if(!fAccCDB) AliInfo(Form("It is not normal because we didn't choose a reference database!"))
3492 scale = fScaleFitFactor;
3494 for(Int_t l = 0; l < nbins; l++){
3495 if(fCoefCharge[1][l] > 0){
3496 fCoefCharge[1][l]=fCoefCharge[1][l]*scale;
3497 fCoefChargeE[1][l]=fCoefChargeE[1][l]*scale;
3502 if(fFitChargeBisOn){
3504 for(Int_t l = 0; l < nbins; l++){
3505 if(fCoefCharge[2][l] > 0){
3506 sum += fCoefCharge[2][l];
3507 xValuesFittedBis[counter[2]]= l;
3512 if(sum > 0.0) scale = counter[2]/sum;
3513 if(fFitChargeNDB == 0){
3514 if(scale != fScaleFitFactor){
3515 AliInfo(Form("The normalisation is different from a nomalisation to one."));
3516 AliInfo(Form("For one we have %f and here %f",scale,fScaleFitFactor));
3517 if(!fAccCDB) AliInfo(Form("It is not normal because we didn't choose a reference database!"))
3519 scale = fScaleFitFactor;
3521 for(Int_t l = 0; l < nbins; l++){
3522 if(fCoefCharge[2][l] > 0){
3523 fCoefCharge[2][l]=fCoefCharge[2][l]*scale;
3524 fCoefChargeE[2][l]=fCoefChargeE[2][l]*scale;
3529 // Create the X and Xerror
3530 Double_t *xValues = new Double_t[nbins];
3531 Double_t *xValuesE = new Double_t[nbins];
3532 for(Int_t k = 0; k < nbins; k ++){
3537 // Create the graph erros and plot them
3538 TCanvas *cch1 = new TCanvas("cch1","",50,50,600,800);
3540 TLegend *legch1 = new TLegend(0.4,0.6,0.89,0.89);
3542 TGraph *graphCharge3 = new TGraph(nbins,xValues,fCoefCharge[3]);
3543 graphCharge3->SetName("coefcharge3");
3544 graphCharge3->SetTitle("");
3545 graphCharge3->GetXaxis()->SetTitle("Det/Pad groups");
3546 graphCharge3->GetYaxis()->SetTitle("gain factor");
3547 graphCharge3->SetLineColor(4);
3548 graphCharge3->SetMarkerStyle(25);
3549 graphCharge3->SetMarkerColor(4);
3550 listofgraphs->Add((TObject *)graphCharge3);
3551 legch1->AddEntry(graphCharge3,"f_{g} simulated","p");
3552 graphCharge3->Draw("AP");
3555 TGraphErrors *graphCharge0 = new TGraphErrors(nbins,xValues,fCoefCharge[0],xValuesE,fCoefChargeE[0]);
3556 graphCharge0->SetName("coefcharge0");
3557 graphCharge0->SetTitle("");
3558 graphCharge0->GetXaxis()->SetTitle("Det/Pad groups");
3559 graphCharge0->GetYaxis()->SetTitle("gain factor");
3560 graphCharge0->SetMarkerColor(6);
3561 graphCharge0->SetLineColor(6);
3562 graphCharge0->SetMarkerStyle(26);
3563 listofgraphs->Add((TObject *)graphCharge0);
3564 legch1->AddEntry(graphCharge0,"f_{g} fit","p");
3565 graphCharge0->Draw("P");
3568 TGraphErrors *graphCharge4 = new TGraphErrors(nbins,xValues,fCoefCharge[4],xValuesE,fCoefChargeE[3]);
3569 graphCharge4->SetName("coefcharge4");
3570 graphCharge4->SetTitle("");
3571 graphCharge4->GetXaxis()->SetTitle("Det/Pad groups");
3572 graphCharge4->GetYaxis()->SetTitle("gain factor");
3573 graphCharge4->SetMarkerColor(1);
3574 graphCharge4->SetLineColor(1);
3575 graphCharge4->SetMarkerStyle(30);
3576 listofgraphs->Add((TObject *)graphCharge4);
3577 legch1->AddEntry(graphCharge4,"f_{g} Mean W","p");
3578 graphCharge4->Draw("P");
3580 if (fMeanChargeOn) {
3581 TGraphErrors *graphCharge1 = new TGraphErrors(nbins,xValues,fCoefCharge[1],xValuesE,fCoefChargeE[1]);
3582 graphCharge1->SetName("coefcharge1");
3583 graphCharge1->GetXaxis()->SetTitle("Det/Pad groups");
3584 graphCharge1->GetYaxis()->SetTitle("gain factor");
3585 graphCharge1->SetTitle("");
3586 graphCharge1->SetMarkerColor(2);
3587 graphCharge1->SetLineColor(2);
3588 graphCharge1->SetMarkerStyle(24);
3589 legch1->AddEntry(graphCharge1,"f_{g} mean","p");
3590 graphCharge1->Draw("P");
3591 listofgraphs->Add((TObject *)graphCharge1);
3593 if (fFitChargeBisOn ) {
3594 TGraphErrors *graphCharge2 = new TGraphErrors(nbins,xValues,fCoefCharge[2],xValuesE,fCoefChargeE[2]);
3595 graphCharge2->SetName("coefcharge2");
3596 graphCharge2->SetTitle("");
3597 graphCharge2->GetXaxis()->SetTitle("Det/Pad groups");
3598 graphCharge2->GetYaxis()->SetTitle("gain factor");
3599 graphCharge2->SetMarkerColor(8);
3600 graphCharge2->SetLineColor(8);
3601 graphCharge2->SetMarkerStyle(25);
3602 legch1->AddEntry(graphCharge2,"f_{g} fitbis","p");
3603 graphCharge2->Draw("P");
3604 listofgraphs->Add((TObject *)graphCharge2);
3606 legch1->Draw("same");
3608 //Create the arrays and the graphs for the delta
3610 TCanvas *cch2 = new TCanvas("cch2","",50,50,600,800);
3612 TLegend *legch3 = new TLegend(0.4,0.6,0.89,0.89);
3613 TLegend *legch2 = new TLegend(0.4,0.6,0.89,0.89);
3617 Double_t *yValuesDelta = new Double_t[counter[0]];
3618 for(Int_t k = 0; k < counter[0]; k++){
3619 if (fCoefCharge[3][(Int_t)(xValuesFitted[k])] > 0.0) {
3620 yValuesDelta[k] = (fCoefCharge[0][(Int_t)xValuesFitted[k]]-fCoefCharge[3][(Int_t)xValuesFitted[k]])
3621 / fCoefCharge[3][(Int_t)xValuesFitted[k]];
3624 yValuesDelta[k] = 0.0;
3627 TGraph *graphDeltaCharge0 = new TGraph(counter[0],&xValuesFitted[0],yValuesDelta);
3628 graphDeltaCharge0->SetName("deltacharge0");
3629 graphDeltaCharge0->GetXaxis()->SetTitle("Det/Pad groups");
3630 graphDeltaCharge0->GetYaxis()->SetTitle("#Deltag/g_{sim}");
3631 graphDeltaCharge0->SetMarkerColor(6);
3632 graphDeltaCharge0->SetTitle("");
3633 graphDeltaCharge0->SetLineColor(6);
3634 graphDeltaCharge0->SetMarkerStyle(26);
3635 listofgraphs->Add((TObject *)graphDeltaCharge0);
3636 legch3->AddEntry(graphDeltaCharge0,"fit","p");
3637 graphDeltaCharge0->Draw("AP");
3640 TH1I *histoErrorCharge0 = new TH1I("errorcharge0","",100 ,-0.10,0.10);
3641 histoErrorCharge0->SetXTitle("#Deltag/g_{sim}");
3642 histoErrorCharge0->SetYTitle("counts");
3643 histoErrorCharge0->SetLineColor(6);
3644 histoErrorCharge0->SetLineStyle(1);
3645 histoErrorCharge0->SetStats(0);
3646 Double_t maxvalue = 0.0;
3647 for(Int_t k = 0; k < counter[0]; k++){
3648 histoErrorCharge0->Fill(yValuesDelta[k]);
3649 if(k == 0) maxvalue = TMath::Abs(yValuesDelta[k]);
3650 if(maxvalue < (TMath::Abs(yValuesDelta[k]))) maxvalue = TMath::Abs(yValuesDelta[k]);
3652 AliInfo(Form("The maximum deviation found dor the fit method is %f",maxvalue));
3653 legch2->AddEntry(histoErrorCharge0,"f_{g} fit","l");
3654 histoErrorCharge0->Draw();
3655 listofgraphs->Add((TObject *)histoErrorCharge0);
3661 Double_t *yValuesDelta = new Double_t[counter[3]];
3662 for(Int_t k = 0; k < counter[3]; k++){
3663 if (fCoefCharge[3][(Int_t)(xValuesFittedMeanW[k])] > 0.0) {
3664 yValuesDelta[k] = (fCoefCharge[4][(Int_t)xValuesFittedMeanW[k]]-fCoefCharge[3][(Int_t)xValuesFittedMeanW[k]])
3665 / fCoefCharge[3][(Int_t)xValuesFittedMeanW[k]];
3668 yValuesDelta[k] = 0.0;
3671 TGraph *graphDeltaCharge4 = new TGraph(counter[3],&xValuesFittedMeanW[0],yValuesDelta);
3672 graphDeltaCharge4->SetName("deltacharge4");
3673 graphDeltaCharge4->GetXaxis()->SetTitle("Det/Pad groups");
3674 graphDeltaCharge4->GetYaxis()->SetTitle("#Deltag/g_{sim}");
3675 graphDeltaCharge4->SetMarkerColor(1);
3676 graphDeltaCharge4->SetTitle("");
3677 graphDeltaCharge4->SetLineColor(1);
3678 graphDeltaCharge4->SetMarkerStyle(30);
3679 listofgraphs->Add((TObject *)graphDeltaCharge4);
3680 legch3->AddEntry(graphDeltaCharge4,"Mean W","p");
3682 graphDeltaCharge4->Draw("AP");
3685 graphDeltaCharge4->Draw("P");
3689 TH1I *histoErrorCharge4 = new TH1I("errorcharge4","",100 ,-0.10,0.10);
3690 histoErrorCharge4->SetXTitle("#Deltag/g_{sim}");
3691 histoErrorCharge4->SetYTitle("counts");
3692 histoErrorCharge4->SetLineColor(1);
3693 histoErrorCharge4->SetLineStyle(1);
3694 histoErrorCharge4->SetStats(0);
3695 Double_t maxvalue = 0.0;
3696 for(Int_t k = 0; k < counter[3]; k++){
3697 histoErrorCharge4->Fill(yValuesDelta[k]);
3698 if(k == 0) maxvalue = yValuesDelta[k];
3699 if(maxvalue < (TMath::Abs(yValuesDelta[k]))) maxvalue = TMath::Abs(yValuesDelta[k]);
3701 AliInfo(Form("The maximum deviation found for the meanW method is %f",maxvalue));
3702 legch2->AddEntry(histoErrorCharge4,"f_{g} Mean W","l");
3704 histoErrorCharge4->Draw();
3707 histoErrorCharge4->Draw("same");
3709 listofgraphs->Add((TObject *)histoErrorCharge4);
3713 if (fMeanChargeOn) {
3715 Double_t *yValuesDeltaMean = new Double_t[counter[1]];
3716 for (Int_t k = 0; k < counter[1]; k++){
3717 if (fCoefCharge[3][(Int_t)xValuesFittedMean[k]] > 0.0) {
3718 yValuesDeltaMean[k] = (fCoefCharge[1][(Int_t)xValuesFittedMean[k]]-fCoefCharge[3][(Int_t)xValuesFittedMean[k]])
3719 / fCoefCharge[3][(Int_t)xValuesFittedMean[k]];
3722 yValuesDeltaMean[k] = 0.0;
3725 TGraph *graphDeltaCharge1 = new TGraph(counter[1],&xValuesFittedMean[0],yValuesDeltaMean);
3726 graphDeltaCharge1->SetName("deltacharge1");
3727 graphDeltaCharge1->GetXaxis()->SetTitle("Det/Pad groups");
3728 graphDeltaCharge1->GetYaxis()->SetTitle("#Deltag/g_{sim}");
3729 graphDeltaCharge1->SetMarkerColor(2);
3730 graphDeltaCharge1->SetMarkerStyle(24);
3731 graphDeltaCharge1->SetLineColor(2);
3732 graphDeltaCharge1->SetTitle("");
3733 legch3->AddEntry(graphDeltaCharge1,"mean","p");
3735 graphDeltaCharge1->Draw("AP");
3738 graphDeltaCharge1->Draw("P");
3740 listofgraphs->Add((TObject *)graphDeltaCharge1);
3743 TH1I *histoErrorCharge1 = new TH1I("errorcharge1","",100 ,-0.10,0.10);
3744 histoErrorCharge1->SetXTitle("#Deltag/g_{sim}");
3745 histoErrorCharge1->SetYTitle("counts");
3746 histoErrorCharge1->SetLineColor(2);
3747 histoErrorCharge1->SetLineStyle(2);
3748 histoErrorCharge1->SetStats(0);
3749 Double_t maxvalue = 0.0;
3750 for(Int_t k = 0; k < counter[1]; k++){
3751 histoErrorCharge1->Fill(yValuesDeltaMean[k]);
3752 if(k == 0) maxvalue = TMath::Abs(yValuesDeltaMean[k]);
3753 if(maxvalue < (TMath::Abs(yValuesDeltaMean[k]))) maxvalue = TMath::Abs(yValuesDeltaMean[k]);
3755 AliInfo(Form("The maximum deviation found for the mean method is %f",maxvalue));
3756 legch2->AddEntry(histoErrorCharge1,"f_{g} mean","l");
3758 histoErrorCharge1->Draw();
3761 histoErrorCharge1->Draw("same");
3763 listofgraphs->Add((TObject *)histoErrorCharge1);
3767 if (fFitChargeBisOn) {
3769 Double_t *yValuesDeltaBis = new Double_t[counter[2]];
3770 for(Int_t k = 0; k < counter[2]; k++){
3771 if (fCoefCharge[3][(Int_t)xValuesFittedBis[k]] > 0.0) {
3772 yValuesDeltaBis[k] = (fCoefCharge[2][(Int_t)xValuesFittedBis[k]]-fCoefCharge[3][(Int_t)xValuesFittedBis[k]])
3773 / fCoefCharge[3][(Int_t)xValuesFittedBis[k]];
3776 yValuesDeltaBis[k] = 0.0;
3779 TGraph *graphDeltaCharge2 = new TGraph(counter[2],&xValuesFittedBis[0],yValuesDeltaBis);
3780 graphDeltaCharge2->SetName("deltacharge2");
3781 graphDeltaCharge2->GetXaxis()->SetTitle("Det/Pad groups");
3782 graphDeltaCharge2->GetYaxis()->SetTitle("#Deltag/g_{sim}");
3783 graphDeltaCharge2->SetMarkerColor(8);
3784 graphDeltaCharge2->SetLineColor(8);
3785 graphDeltaCharge2->SetMarkerStyle(25);
3786 legch3->AddEntry(graphDeltaCharge2,"fit","p");
3787 graphDeltaCharge2->SetTitle("");
3789 graphDeltaCharge2->Draw("AP");
3792 graphDeltaCharge2->Draw("P");
3794 listofgraphs->Add((TObject *)graphDeltaCharge2);
3797 TH1I *histoErrorCharge2 = new TH1I("errorcharge2","",100 ,-0.10, 0.10);
3798 histoErrorCharge2->SetXTitle("#Deltag/g_{sim}");
3799 histoErrorCharge2->SetYTitle("counts");
3800 histoErrorCharge2->SetLineColor(8);
3801 histoErrorCharge2->SetLineStyle(5);
3802 histoErrorCharge2->SetLineWidth(3);
3803 histoErrorCharge2->SetStats(0);
3804 Double_t maxvalue = 0.0;
3805 for(Int_t k = 0; k < counter[2]; k++){
3806 histoErrorCharge2->Fill(yValuesDeltaBis[k]);
3807 if(k == 0) maxvalue = TMath::Abs(yValuesDeltaBis[k]);
3808 if(maxvalue < (TMath::Abs(yValuesDeltaBis[k]))) maxvalue = TMath::Abs(yValuesDeltaBis[k]);
3810 AliInfo(Form("The maximum deviation found for the fit bis method is %f",maxvalue));
3811 legch2->AddEntry(histoErrorCharge2,"f_{g} fitbis","l");
3813 histoErrorCharge2->Draw();
3816 histoErrorCharge2->Draw("same");
3818 listofgraphs->Add((TObject *)histoErrorCharge2);
3819 //it doesn't matter anymore but...
3824 legch3->Draw("same");
3826 legch2->Draw("same");
3830 TFile *fout = TFile::Open(fWriteNameCoef,"UPDATE");
3831 // Check if the file could be opened
3832 if (!fout || !fout->IsOpen()) {
3833 AliInfo("No File found!");
3837 for(Int_t k = 0; k < listofgraphs->GetEntriesFast(); k++){
3838 fout->WriteTObject((TObject *) listofgraphs->At(k),((TObject *)listofgraphs->At(k))->GetName()
3839 ,(Option_t *) "OverWrite");
3847 //_____________________________________________________________________________
3848 void AliTRDCalibraFit::PlotWritePH()
3851 // create the graph errors and write them if wanted
3854 //TObjArray of the grapherrors and so on
3855 TObjArray *listofgraphs = new TObjArray();
3857 Int_t nbins = fDect2[1]-fDect1[1];
3859 //See the number of fitted for delta
3867 Double_t *xValuesFitted = new Double_t[nbins];
3868 Double_t *xValuesFittedPH = new Double_t[nbins];
3869 Double_t *xValuesFittedLP = new Double_t[nbins];
3870 for(Int_t k = 0; k < nbins; k ++){
3871 xValuesFitted[k] = -1;
3872 xValuesFittedPH[k] = -1;
3873 xValuesFittedLP[k] = -1;
3877 for(Int_t l = 0; l < nbins; l++){
3878 if(fCoefVdrift[1][l] > 0){
3879 xValuesFitted[counter[1]]=l;
3885 for(Int_t l = 0; l < nbins; l++){
3886 if(fCoefVdrift[3][l] > 0){
3887 xValuesFittedLP[counter[2]]=l;
3893 for(Int_t l = 0; l < nbins; l++){
3894 if(fCoefVdrift[0][l] > 0){
3895 xValuesFittedPH[counter[0]]= l;
3901 //Create the X and Xerror
3902 Double_t *xValues = new Double_t[nbins];
3903 Double_t *xValuesE = new Double_t[nbins];
3904 for(Int_t k = 0; k < nbins; k ++){
3909 //Create the graph erros and plot them
3910 TCanvas *cph1 = new TCanvas("cph1","",50,50,600,800);
3913 TGraph *graphVdrift2 = new TGraph(nbins,xValues,fCoefVdrift[2]);
3914 graphVdrift2->SetName("coefvdrift2");
3915 graphVdrift2->SetTitle("");
3916 graphVdrift2->GetXaxis()->SetTitle("Det/Pad groups");
3917 graphVdrift2->GetYaxis()->SetTitle("Vdrift [cm/#mus]");
3918 graphVdrift2->SetLineColor(4);
3919 listofgraphs->Add((TObject *)graphVdrift2);
3920 TLegend *legph1 = new TLegend(0.4,0.6,0.89,0.89);
3921 legph1->AddEntry(graphVdrift2,"Vdrift simulated","l");
3922 graphVdrift2->Draw("AL");
3925 TGraphErrors *graphVdrift1 = new TGraphErrors(nbins,xValues,fCoefVdrift[1],xValuesE,fCoefVdriftE[1]);
3926 graphVdrift1->SetName("coefvdrift1");
3927 graphVdrift1->SetTitle("");
3928 graphVdrift1->GetXaxis()->SetTitle("Det/Pad groups");
3929 graphVdrift1->GetYaxis()->SetTitle("Vdrift [cm/#mus]");
3930 graphVdrift1->SetMarkerColor(6);
3931 graphVdrift1->SetLineColor(6);
3932 graphVdrift1->SetMarkerStyle(26);
3933 listofgraphs->Add((TObject *)graphVdrift1);
3934 legph1->AddEntry(graphVdrift1,"Vdrift fit","p");
3935 graphVdrift1->Draw("P");
3938 TGraphErrors *graphVdrift0 = new TGraphErrors(nbins,xValues,fCoefVdrift[0],xValuesE,fCoefVdriftE[0]);
3939 graphVdrift0->SetName("coefVdrift0");
3940 graphVdrift0->GetXaxis()->SetTitle("Det/Pad groups");
3941 graphVdrift0->GetYaxis()->SetTitle("Vdrift [cm/#mus]");
3942 graphVdrift0->SetTitle("");
3943 graphVdrift0->SetMarkerColor(2);
3944 graphVdrift0->SetLineColor(2);
3945 graphVdrift0->SetMarkerStyle(24);
3946 legph1->AddEntry(graphVdrift0,"v_{fit PH}","p");
3947 graphVdrift0->Draw("P");
3948 listofgraphs->Add((TObject *)graphVdrift0);
3950 if (fFitLagrPolOn) {
3951 TGraphErrors *graphVdrift3 = new TGraphErrors(nbins,xValues,fCoefVdrift[3],xValuesE,fCoefVdriftE[2]);
3952 graphVdrift3->SetName("coefVdrift3");
3953 graphVdrift3->GetXaxis()->SetTitle("Det/Pad groups");
3954 graphVdrift3->GetYaxis()->SetTitle("Vdrift [cm/#mus]");
3955 graphVdrift3->SetTitle("");
3956 graphVdrift3->SetMarkerColor(1);
3957 graphVdrift3->SetLineColor(1);
3958 graphVdrift3->SetMarkerStyle(28);
3959 legph1->AddEntry(graphVdrift3,"v_{LagrPol}","p");
3960 graphVdrift3->Draw("P");
3961 listofgraphs->Add((TObject *)graphVdrift3);
3963 legph1->Draw("same");
3965 //Create the arrays and the graphs for the delta
3966 TCanvas *cph2 = new TCanvas("cph2","",50,50,600,800);
3968 TLegend *legph3 = new TLegend(0.4,0.6,0.89,0.89);
3969 TLegend *legph2 = new TLegend(0.4,0.6,0.89,0.89);
3974 Double_t *yValuesDelta = new Double_t[counter[1]];
3975 for (Int_t k = 0; k < counter[1]; k++){
3976 if (fCoefVdrift[2][(Int_t)(xValuesFitted[k])] > 0.0) {
3977 yValuesDelta[k] = (fCoefVdrift[1][(Int_t)xValuesFitted[k]]-fCoefVdrift[2][(Int_t)xValuesFitted[k]])
3978 / fCoefVdrift[2][(Int_t)xValuesFitted[k]];
3981 yValuesDelta[k] = 0.0;
3984 TGraph *graphDeltaVdrift1 = new TGraph(counter[1],&xValuesFitted[0],yValuesDelta);
3985 graphDeltaVdrift1->SetName("deltavdrift1");
3986 graphDeltaVdrift1->GetXaxis()->SetTitle("Det/Pad groups");
3987 graphDeltaVdrift1->GetYaxis()->SetTitle("#Deltav/v_{sim}");
3988 graphDeltaVdrift1->SetMarkerColor(6);
3989 graphDeltaVdrift1->SetTitle("");
3990 graphDeltaVdrift1->SetLineColor(6);
3991 graphDeltaVdrift1->SetMarkerStyle(26);
3992 listofgraphs->Add((TObject *)graphDeltaVdrift1);
3993 legph3->AddEntry(graphDeltaVdrift1,"v_{slope method}","p");
3994 graphDeltaVdrift1->Draw("AP");
3997 TH1I *histoErrorVdrift1 = new TH1I("errorvdrift1","",100 ,-0.2,0.2);
3998 histoErrorVdrift1->SetXTitle("#Deltav/v_{sim}");
3999 histoErrorVdrift1->SetYTitle("counts");
4000 histoErrorVdrift1->SetLineColor(6);
4001 histoErrorVdrift1->SetLineStyle(1);
4002 histoErrorVdrift1->SetStats(0);
4003 Double_t maxvalue = 0.0;
4004 for(Int_t k = 0; k < counter[1]; k++){
4005 histoErrorVdrift1->Fill(yValuesDelta[k]);
4006 if(k == 0) maxvalue = yValuesDelta[k];
4007 if(maxvalue < (TMath::Abs(yValuesDelta[k]))) maxvalue = TMath::Abs(yValuesDelta[k]);
4009 AliInfo(Form("The maximum deviation found for the Pol2 method is %f",maxvalue));
4010 legph2->AddEntry(histoErrorVdrift1,"v_{slope method}","l");
4011 histoErrorVdrift1->Draw();
4012 listofgraphs->Add((TObject *)histoErrorVdrift1);
4018 Double_t *yValuesDeltaPH = new Double_t[counter[0]];
4019 for(Int_t k = 0; k < counter[0]; k++){
4020 if(fCoefVdrift[2][(Int_t)xValuesFittedPH[k]] > 0.0) {
4021 yValuesDeltaPH[k] = (fCoefVdrift[0][(Int_t)xValuesFittedPH[k]]-fCoefVdrift[2][(Int_t)xValuesFittedPH[k]])/fCoefVdrift[2][(Int_t)xValuesFittedPH[k]];
4023 else yValuesDeltaPH[k] = 0.0;
4025 TGraph *graphDeltaVdrift0 = new TGraph(counter[0],&xValuesFittedPH[0],yValuesDeltaPH);
4026 graphDeltaVdrift0->SetName("deltavdrift0");
4027 graphDeltaVdrift0->GetXaxis()->SetTitle("Det/Pad groups");
4028 graphDeltaVdrift0->GetYaxis()->SetTitle("#Deltav/v_{sim}");
4029 graphDeltaVdrift0->SetMarkerColor(2);
4030 graphDeltaVdrift0->SetMarkerStyle(24);
4031 graphDeltaVdrift0->SetLineColor(2);
4032 graphDeltaVdrift0->SetTitle("");
4033 legph3->AddEntry(graphDeltaVdrift0,"v_{fit PH}","p");
4035 graphDeltaVdrift0->Draw("P");
4038 graphDeltaVdrift0->Draw("AP");
4040 listofgraphs->Add((TObject *)graphDeltaVdrift0);
4042 TH1I *histoErrorVdrift0 = new TH1I("errorvdrift0","",100 ,-0.2,0.2);
4043 histoErrorVdrift0->SetXTitle("#Deltav/v_{sim}");
4044 histoErrorVdrift0->SetYTitle("counts");
4045 histoErrorVdrift0->SetLineColor(2);
4046 histoErrorVdrift0->SetLineStyle(2);
4047 histoErrorVdrift0->SetStats(0);
4048 Double_t maxvalue = 0.0;
4049 for(Int_t k = 0; k < counter[0]; k++){
4050 histoErrorVdrift0->Fill(yValuesDeltaPH[k]);
4051 if(k == 0) maxvalue = yValuesDeltaPH[k];
4052 if(maxvalue < (TMath::Abs(yValuesDeltaPH[k]))) maxvalue = TMath::Abs(yValuesDeltaPH[k]);
4054 AliInfo(Form("The maximum deviation found for the fit method is %f",maxvalue));
4055 legph2->AddEntry(histoErrorVdrift0,"v_{fit PH}","l");
4057 histoErrorVdrift0->Draw("same");
4060 histoErrorVdrift0->Draw();
4062 listofgraphs->Add((TObject *)histoErrorVdrift0);
4066 if (fFitLagrPolOn) {
4068 Double_t *yValuesDeltaPH = new Double_t[counter[2]];
4069 for (Int_t k = 0; k < counter[2]; k++){
4070 if (fCoefVdrift[2][(Int_t)xValuesFittedLP[k]] > 0.0) {
4071 yValuesDeltaPH[k] = (fCoefVdrift[3][(Int_t)xValuesFittedLP[k]]-fCoefVdrift[2][(Int_t)xValuesFittedLP[k]])
4072 / fCoefVdrift[2][(Int_t)xValuesFittedLP[k]];
4075 yValuesDeltaPH[k] = 0.0;
4078 TGraph *graphDeltaVdrift3 = new TGraph(counter[2],&xValuesFittedLP[0],yValuesDeltaPH);
4079 graphDeltaVdrift3->SetName("deltavdrift3");
4080 graphDeltaVdrift3->GetXaxis()->SetTitle("Det/Pad groups");
4081 graphDeltaVdrift3->GetYaxis()->SetTitle("#Deltav/v_{sim}");
4082 graphDeltaVdrift3->SetMarkerColor(1);
4083 graphDeltaVdrift3->SetMarkerStyle(28);
4084 graphDeltaVdrift3->SetLineColor(1);
4085 graphDeltaVdrift3->SetTitle("");
4086 legph3->AddEntry(graphDeltaVdrift3,"v_{LagrPol}","p");
4088 graphDeltaVdrift3->Draw("P");
4091 graphDeltaVdrift3->Draw("AP");
4093 listofgraphs->Add((TObject *)graphDeltaVdrift3);
4095 TH1I *histoErrorVdrift3 = new TH1I("errorvdrift3","",100 ,-0.2,0.2);
4096 histoErrorVdrift3->SetXTitle("#Deltav/v_{sim}");
4097 histoErrorVdrift3->SetYTitle("counts");
4098 histoErrorVdrift3->SetLineColor(1);
4099 histoErrorVdrift3->SetLineStyle(1);
4100 histoErrorVdrift3->SetStats(0);
4101 Double_t maxvalue = 0.0;
4102 for(Int_t k = 0; k < counter[2]; k++){
4103 histoErrorVdrift3->Fill(yValuesDeltaPH[k]);
4104 if(k == 0) maxvalue = yValuesDeltaPH[k];
4105 if(maxvalue < (TMath::Abs(yValuesDeltaPH[k]))) maxvalue = TMath::Abs(yValuesDeltaPH[k]);
4107 AliInfo(Form("The maximum deviation found for the LagrPol method is %f",maxvalue));
4108 legph2->AddEntry(histoErrorVdrift3,"v_{LagrPol}","l");
4110 histoErrorVdrift3->Draw("same");
4113 histoErrorVdrift3->Draw();
4115 listofgraphs->Add((TObject *)histoErrorVdrift3);
4119 legph3->Draw("same");
4121 legph2->Draw("same");
4125 TFile *fout = TFile::Open(fWriteNameCoef,"UPDATE");
4126 // Check if the file could be opened
4127 if (!fout || !fout->IsOpen()) {
4128 AliInfo("No File found!");
4132 for(Int_t k = 0; k < listofgraphs->GetEntriesFast(); k++){
4133 fout->WriteTObject((TObject *) listofgraphs->At(k),((TObject *)listofgraphs->At(k))->GetName()
4134 ,(Option_t *) "OverWrite");
4142 //_____________________________________________________________________________
4143 void AliTRDCalibraFit::PlotWriteT0()
4146 // create the graph errors and write them if wanted
4149 //TObjArray of the grapherrors and so on
4150 TObjArray *listofgraphs = new TObjArray();
4152 Int_t nbins = fDect2[1]-fDect1[1];
4154 //See the number of fitted for delta: here T0 can be negative, we don't use the sign but the error
4155 //and the grapherrors of the coefficients contained the no fitted with error 0.0
4163 Double_t *xValuesFitted = new Double_t[nbins];
4164 Double_t *xValuesFittedPH = new Double_t[nbins];
4165 Double_t *xValuesFittedLP = new Double_t[nbins];
4166 for(Int_t k = 0; k < nbins; k ++){
4167 xValuesFitted[k] = -1;
4168 xValuesFittedPH[k] = -1;
4169 xValuesFittedLP[k] = -1;
4173 for(Int_t l = 0; l < nbins; l++){
4174 if(fCoefT0E[1][l] != 0.0){
4175 xValuesFitted[counter[1]]=l;
4182 for(Int_t l = 0; l < nbins; l++){
4183 if(fCoefT0E[0][l] != 0.0){
4184 xValuesFittedPH[counter[0]]= l;
4191 for(Int_t l = 0; l < nbins; l++){
4192 if(fCoefT0E[2][l] == 1.0){
4193 xValuesFittedLP[counter[2]]= l;
4199 //Create the X and Xerror
4200 Double_t *xValues = new Double_t[nbins];
4201 Double_t *xValuesE = new Double_t[nbins];
4202 for(Int_t k = 0; k < nbins; k ++){
4207 //Create the graph erros and plot them
4208 TCanvas *ct01 = new TCanvas("ct01","",50,50,600,800);
4210 TLegend *legt01 = new TLegend(0.4,0.6,0.89,0.89);
4212 TGraph *graphT02 = new TGraph(nbins,xValues,fCoefT0[2]);
4213 graphT02->SetName("coeft02");
4214 graphT02->SetTitle("");
4215 graphT02->GetXaxis()->SetTitle("Det/Pad groups");
4216 graphT02->GetYaxis()->SetTitle("T0 [time bins]");
4217 graphT02->SetLineColor(4);
4218 listofgraphs->Add((TObject *)graphT02);
4219 legt01->AddEntry(graphT02,"T0 simulated","l");
4220 graphT02->Draw("AL");
4223 TGraphErrors *graphT01 = new TGraphErrors(nbins,xValues,fCoefT0[1],xValuesE,fCoefT0E[1]);
4224 graphT01->SetName("coeft01");
4225 graphT01->SetTitle("");
4226 graphT01->GetXaxis()->SetTitle("Det/Pad groups");
4227 graphT01->GetYaxis()->SetTitle("T0 [time bins]");
4228 graphT01->SetMarkerColor(6);
4229 graphT01->SetLineColor(6);
4230 graphT01->SetMarkerStyle(26);
4231 listofgraphs->Add((TObject *)graphT01);
4232 legt01->AddEntry(graphT01,"T0 slope method","p");
4233 graphT01->Draw("P");
4236 TGraphErrors *graphT00 = new TGraphErrors(nbins,xValues,fCoefT0[0],xValuesE,fCoefT0E[0]);
4237 graphT00->SetName("coeft00");
4238 graphT00->GetXaxis()->SetTitle("Det/Pad groups");
4239 graphT00->GetYaxis()->SetTitle("T0 [time bins]");
4240 graphT00->SetTitle("");
4241 graphT00->SetMarkerColor(2);
4242 graphT00->SetLineColor(2);
4243 graphT00->SetMarkerStyle(24);
4244 legt01->AddEntry(graphT00,"T0 fit","p");
4245 graphT00->Draw("P");
4246 listofgraphs->Add((TObject *)graphT00);
4248 if (fFitLagrPolOn) {
4249 TGraphErrors *graphT03 = new TGraphErrors(nbins,xValues,fCoefT0[3],xValuesE,xValuesE);
4250 graphT03->SetName("coeft03");
4251 graphT03->GetXaxis()->SetTitle("Det/Pad groups");
4252 graphT03->GetYaxis()->SetTitle("T0 [time bins]");
4253 graphT03->SetTitle("");
4254 graphT03->SetMarkerColor(1);
4255 graphT03->SetLineColor(1);
4256 graphT03->SetMarkerStyle(28);
4257 legt01->AddEntry(graphT03,"T0 LagrPol","p");
4258 graphT03->Draw("P");
4259 listofgraphs->Add((TObject *)graphT03);
4261 legt01->Draw("same");
4263 //Create the arrays and the graphs for the delta
4264 TCanvas *ct02 = new TCanvas("ct02","",50,50,600,800);
4266 TLegend *legt03 = new TLegend(0.4,0.6,0.89,0.89);
4267 TLegend *legt02 = new TLegend(0.4,0.6,0.89,0.89);
4272 Double_t *yValuesDelta = new Double_t[counter[1]];
4273 for(Int_t k = 0; k < counter[1]; k++){
4274 yValuesDelta[k] = (fCoefT0[1][(Int_t)xValuesFitted[k]]-fCoefT0[2][(Int_t)xValuesFitted[k]]);
4276 TGraph *graphDeltaT01 = new TGraph(counter[1],&xValuesFitted[0],yValuesDelta);
4277 graphDeltaT01->SetName("deltat01");
4278 graphDeltaT01->GetXaxis()->SetTitle("Det/Pad groups");
4279 graphDeltaT01->GetYaxis()->SetTitle("#Deltat0 [time bins]");
4280 graphDeltaT01->SetMarkerColor(6);
4281 graphDeltaT01->SetTitle("");
4282 graphDeltaT01->SetLineColor(6);
4283 graphDeltaT01->SetMarkerStyle(26);
4284 listofgraphs->Add((TObject *)graphDeltaT01);
4285 legt03->AddEntry(graphDeltaT01,"T0_{slope method}","p");
4286 graphDeltaT01->Draw("AP");
4289 TH1I *histoErrorT01 = new TH1I("errort01","",100 ,-0.2,0.2);
4290 histoErrorT01->SetXTitle("#Deltat0 [time bins]");
4291 histoErrorT01->SetYTitle("counts");
4292 histoErrorT01->SetLineColor(6);
4293 histoErrorT01->SetLineStyle(1);
4294 histoErrorT01->SetStats(0);
4295 Double_t maxvalue = 0.0;
4296 for(Int_t k = 0; k < counter[1]; k++){
4297 histoErrorT01->Fill(yValuesDelta[k]);
4298 if(k == 0) maxvalue = yValuesDelta[k];
4299 if(maxvalue < (TMath::Abs(yValuesDelta[k]))) maxvalue = (TMath::Abs(yValuesDelta[k]));
4301 AliInfo(Form("The maximum deviation found for the Pol2 method is %f",maxvalue));
4302 legt02->AddEntry(histoErrorT01,"T0_{slope method}","l");
4303 histoErrorT01->Draw();
4304 listofgraphs->Add((TObject *)histoErrorT01);
4309 Double_t *yValuesDeltaPH = new Double_t[counter[0]];
4310 for(Int_t k = 0; k < counter[0]; k++){
4311 yValuesDeltaPH[k] = (fCoefT0[0][(Int_t)xValuesFittedPH[k]]-fCoefT0[2][(Int_t)xValuesFittedPH[k]]);
4313 TGraph *graphDeltaT00 = new TGraph(counter[0],&xValuesFittedPH[0],yValuesDeltaPH);
4314 graphDeltaT00->SetName("deltat00");
4315 graphDeltaT00->GetXaxis()->SetTitle("Det/Pad groups");
4316 graphDeltaT00->GetYaxis()->SetTitle("#Deltat0 [time bins]");
4317 graphDeltaT00->SetMarkerColor(2);
4318 graphDeltaT00->SetMarkerStyle(24);
4319 graphDeltaT00->SetLineColor(2);
4320 graphDeltaT00->SetTitle("");
4321 legt03->AddEntry(graphDeltaT00,"T0_{fit PH}","p");
4323 graphDeltaT00->Draw("P");
4326 graphDeltaT00->Draw("AP");
4328 listofgraphs->Add((TObject *)graphDeltaT00);
4330 TH1I *histoErrorT00 = new TH1I("errort00","",100 ,-0.2,0.2);
4331 histoErrorT00->SetXTitle("#Deltat0 [time bins]");
4332 histoErrorT00->SetYTitle("counts");
4333 histoErrorT00->SetLineColor(2);
4334 histoErrorT00->SetLineStyle(2);
4335 histoErrorT00->SetStats(0);
4336 Double_t maxvalue = 0.0;
4337 for(Int_t k = 0; k < counter[0]; k++){
4338 histoErrorT00->Fill(yValuesDeltaPH[k]);
4339 if(k == 0) maxvalue = yValuesDeltaPH[k];
4340 if(maxvalue < (TMath::Abs(yValuesDeltaPH[k]))) maxvalue = (TMath::Abs(yValuesDeltaPH[k]));
4342 AliInfo(Form("The maximum deviation found for the fit method is %f",maxvalue));
4343 legt02->AddEntry(histoErrorT00,"T0_{fit PH}","l");
4345 histoErrorT00->Draw("same");
4348 histoErrorT00->Draw();
4350 listofgraphs->Add((TObject *)histoErrorT00);
4354 if (fFitLagrPolOn) {
4356 Double_t *yValuesDeltaPH = new Double_t[counter[2]];
4357 for(Int_t k = 0; k < counter[2]; k++){
4358 yValuesDeltaPH[k] = (fCoefT0[3][(Int_t)xValuesFittedLP[k]]-fCoefT0[2][(Int_t)xValuesFittedLP[k]]);
4360 TGraph *graphDeltaT03 = new TGraph(counter[2],&xValuesFittedLP[0],yValuesDeltaPH);
4361 graphDeltaT03->SetName("deltat03");
4362 graphDeltaT03->GetXaxis()->SetTitle("Det/Pad groups");
4363 graphDeltaT03->GetYaxis()->SetTitle("#Deltat0 [time bins]");
4364 graphDeltaT03->SetMarkerColor(1);
4365 graphDeltaT03->SetMarkerStyle(28);
4366 graphDeltaT03->SetLineColor(1);
4367 graphDeltaT03->SetTitle("");
4368 legt03->AddEntry(graphDeltaT03,"T0_{LagrPol}","p");
4370 graphDeltaT03->Draw("P");
4373 graphDeltaT03->Draw("AP");
4375 listofgraphs->Add((TObject *)graphDeltaT03);
4377 TH1I *histoErrorT03 = new TH1I("errort03","",100 ,-0.2,0.2);
4378 histoErrorT03->SetXTitle("#Deltat0 [time bins]");
4379 histoErrorT03->SetYTitle("counts");
4380 histoErrorT03->SetLineColor(1);
4381 histoErrorT03->SetLineStyle(1);
4382 histoErrorT03->SetStats(0);
4383 Double_t maxvalue = 0.0;
4384 for(Int_t k = 0; k < counter[2]; k++){
4385 histoErrorT03->Fill(yValuesDeltaPH[k]);
4386 if(k == 0) maxvalue = yValuesDeltaPH[k];
4387 if(maxvalue < (TMath::Abs(yValuesDeltaPH[k]))) maxvalue = (TMath::Abs(yValuesDeltaPH[k]));
4389 AliInfo(Form("The maximum deviation found for the LagrPol method is %f",maxvalue));
4390 legt02->AddEntry(histoErrorT03,"T0_{LagrPol}","l");
4392 histoErrorT03->Draw("same");
4395 histoErrorT03->Draw();
4397 listofgraphs->Add((TObject *)histoErrorT03);
4402 legt03->Draw("same");
4404 legt02->Draw("same");
4408 TFile *fout = TFile::Open(fWriteNameCoef,"UPDATE");
4409 // Check if the file could be opened
4410 if (!fout || !fout->IsOpen()) {
4411 AliInfo("No File found!");
4415 for(Int_t k = 0; k < listofgraphs->GetEntriesFast(); k++){
4416 fout->WriteTObject((TObject *) listofgraphs->At(k),((TObject *)listofgraphs->At(k))->GetName()
4417 ,(Option_t *) "OverWrite");
4425 //_____________________________________________________________________________
4426 void AliTRDCalibraFit::PlotWritePRF()
4429 // create the graph errors and write them if wanted
4432 //TObjArray of the grapherrors and so on
4433 TObjArray *listofgraphs = new TObjArray();
4435 Int_t nbins = fDect2[2]-fDect1[2];
4437 //See the number of fitted for delta
4444 Double_t *xValuesFitted = new Double_t[nbins];
4445 for(Int_t k = 0; k < nbins; k ++){
4446 xValuesFitted[k] = -1;
4448 Double_t *xValuesRMS = new Double_t[nbins];
4449 for(Int_t k = 0; k < nbins; k ++){
4454 for(Int_t l = 0; l < nbins; l++){
4455 if(fCoefPRF[0][l] > 0){
4456 xValuesFitted[counter[0]]=l;
4462 for(Int_t l = 0; l < nbins; l++){
4463 if(fCoefPRF[2][l] > 0){
4464 xValuesRMS[counter[1]]=l;
4471 //Create the X and Xerror
4472 Double_t *xValues = new Double_t[nbins];
4473 Double_t *xValuesE = new Double_t[nbins];
4474 for(Int_t k = 0; k < nbins; k ++){
4479 //Create the graph erros and plot them
4480 TCanvas *cprf1 = new TCanvas("cprf1","",50,50,600,800);
4482 TLegend *legprf1 = new TLegend(0.4,0.6,0.89,0.89);
4484 TGraph *graphPRF1 = new TGraph(nbins,xValues,fCoefPRF[1]);
4485 graphPRF1->SetName("coefprf1");
4486 graphPRF1->SetTitle("");
4487 graphPRF1->GetXaxis()->SetTitle("Det/Pad groups");
4488 graphPRF1->GetYaxis()->SetTitle("PRF width [p.u]");
4489 graphPRF1->SetLineColor(4);
4490 graphPRF1->SetMarkerColor(4);
4491 graphPRF1->SetMarkerStyle(25);
4492 graphPRF1->SetMarkerSize(0.7);
4493 listofgraphs->Add((TObject *)graphPRF1);
4494 legprf1->AddEntry(graphPRF1,"PRF width simulated","p");
4495 graphPRF1->Draw("AP");
4498 TGraphErrors *graphPRF0 = new TGraphErrors(nbins,xValues,fCoefPRF[0],xValuesE,fCoefPRFE[0]);
4499 graphPRF0->SetName("coefprf0");
4500 graphPRF0->SetTitle("");
4501 graphPRF0->GetXaxis()->SetTitle("Det/Pad groups");
4502 graphPRF0->GetYaxis()->SetTitle("PRF Width [p.u]");
4503 graphPRF0->SetMarkerColor(6);
4504 graphPRF0->SetLineColor(6);
4505 graphPRF0->SetMarkerStyle(26);
4506 listofgraphs->Add((TObject *)graphPRF0);
4507 legprf1->AddEntry(graphPRF0,"PRF fit","p");
4508 graphPRF0->Draw("P");
4511 TGraphErrors *graphPRF2 = new TGraphErrors(nbins,xValues,fCoefPRF[2],xValuesE,fCoefPRFE[1]);
4512 graphPRF2->SetName("coefprf2");
4513 graphPRF2->SetTitle("");
4514 graphPRF2->GetXaxis()->SetTitle("Det/Pad groups");
4515 graphPRF2->GetYaxis()->SetTitle("PRF Width [p.u]");
4516 graphPRF2->SetMarkerColor(1);
4517 graphPRF2->SetLineColor(1);
4518 graphPRF2->SetMarkerStyle(28);
4519 listofgraphs->Add((TObject *)graphPRF2);
4520 legprf1->AddEntry(graphPRF2,"PRF Rms","p");
4521 graphPRF2->Draw("P");
4523 legprf1->Draw("same");
4526 //Create the arrays and the graphs for the delta
4527 TCanvas *cprf2 = new TCanvas("cprf2","",50,50,600,800);
4530 TLegend *legprf3 = new TLegend(0.4,0.6,0.89,0.89);
4531 TLegend *legprf2 = new TLegend(0.4,0.6,0.89,0.89);
4535 Double_t *yValuesDelta = new Double_t[counter[0]];
4536 for(Int_t k = 0; k < counter[0]; k++){
4537 if(fCoefPRF[1][(Int_t)xValuesFitted[k]] > 0.0){
4538 yValuesDelta[k] = (fCoefPRF[0][(Int_t)xValuesFitted[k]]-fCoefPRF[1][(Int_t)xValuesFitted[k]])
4539 / (fCoefPRF[1][(Int_t)xValuesFitted[k]]);
4542 TGraph *graphDeltaPRF0 = new TGraph(counter[0],&xValuesFitted[0],yValuesDelta);
4543 graphDeltaPRF0->SetName("deltaprf0");
4544 graphDeltaPRF0->GetXaxis()->SetTitle("Det/Pad groups");
4545 graphDeltaPRF0->GetYaxis()->SetTitle("#Delta#sigma/#sigma_{sim}");
4546 graphDeltaPRF0->SetMarkerColor(6);
4547 graphDeltaPRF0->SetTitle("");
4548 graphDeltaPRF0->SetLineColor(6);
4549 graphDeltaPRF0->SetMarkerStyle(26);
4550 listofgraphs->Add((TObject *)graphDeltaPRF0);
4551 legprf3->AddEntry(graphDeltaPRF0,"#sigma_{fit}","p");
4552 graphDeltaPRF0->Draw("AP");
4555 TH1I *histoErrorPRF0 = new TH1I("errorprf10","",100 ,-0.1,0.2);
4556 histoErrorPRF0->SetXTitle("#Delta#sigma/#sigma_{sim}");
4557 histoErrorPRF0->SetYTitle("counts");
4558 histoErrorPRF0->SetLineColor(6);
4559 histoErrorPRF0->SetLineStyle(1);
4560 histoErrorPRF0->SetStats(0);
4561 Double_t maxvalue = 0.0;
4562 for(Int_t k = 0; k < counter[0]; k++){
4563 histoErrorPRF0->Fill(yValuesDelta[k]);
4564 if(k == 0) maxvalue = yValuesDelta[k];
4565 if(maxvalue < (TMath::Abs(yValuesDelta[k]))) maxvalue = (TMath::Abs(yValuesDelta[k]));
4567 AliInfo(Form("The maximum deviation for the fit method is %f",maxvalue));
4568 legprf2->AddEntry(histoErrorPRF0,"#sigma_{fit}","l");
4569 histoErrorPRF0->Draw();
4570 listofgraphs->Add((TObject *)histoErrorPRF0);
4576 Double_t *yValuesDelta = new Double_t[counter[1]];
4577 for(Int_t k = 0; k < counter[1]; k++){
4578 if(fCoefPRF[1][(Int_t)xValuesRMS[k]] > 0.0){
4579 yValuesDelta[k] = (fCoefPRF[2][(Int_t)xValuesRMS[k]]-fCoefPRF[1][(Int_t)xValuesRMS[k]])
4580 / (fCoefPRF[1][(Int_t)xValuesRMS[k]]);
4583 TGraph *graphDeltaPRF2 = new TGraph(counter[1],&xValuesRMS[0],yValuesDelta);
4584 graphDeltaPRF2->SetName("deltaprf2");
4585 graphDeltaPRF2->GetXaxis()->SetTitle("Det/Pad groups");
4586 graphDeltaPRF2->GetYaxis()->SetTitle("#Delta#sigma/#sigma_{sim}");
4587 graphDeltaPRF2->SetMarkerColor(1);
4588 graphDeltaPRF2->SetTitle("");
4589 graphDeltaPRF2->SetLineColor(1);
4590 graphDeltaPRF2->SetMarkerStyle(28);
4591 listofgraphs->Add((TObject *)graphDeltaPRF2);
4592 legprf3->AddEntry(graphDeltaPRF2,"#sigma_{rms}","p");
4594 graphDeltaPRF2->Draw("P");
4597 graphDeltaPRF2->Draw("AP");
4601 TH1I *histoErrorPRF2 = new TH1I("errorprf12","",100 ,-0.1,0.2);
4602 histoErrorPRF2->SetXTitle("#Delta#sigma/#sigma_{sim}");
4603 histoErrorPRF2->SetYTitle("counts");
4604 histoErrorPRF2->SetLineColor(1);
4605 histoErrorPRF2->SetLineStyle(1);
4606 histoErrorPRF2->SetStats(0);
4607 Double_t maxvalue = 0.0;
4608 for(Int_t k = 0; k < counter[1]; k++){
4609 histoErrorPRF2->Fill(yValuesDelta[k]);
4610 if(k == 0) maxvalue = yValuesDelta[k];
4611 if(maxvalue < TMath::Abs(yValuesDelta[k])) maxvalue = TMath::Abs(yValuesDelta[k]);
4613 AliInfo(Form("The maximum deviation for the rms is %f",maxvalue));
4614 legprf2->AddEntry(histoErrorPRF2,"#sigma_{rms}","l");
4616 histoErrorPRF2->Draw("same");
4619 histoErrorPRF2->Draw();
4621 listofgraphs->Add((TObject *)histoErrorPRF2);
4626 legprf3->Draw("same");
4628 legprf2->Draw("same");
4633 TFile *fout = TFile::Open(fWriteNameCoef,"UPDATE");
4634 // Check if the file could be opened
4635 if (!fout || !fout->IsOpen()) {
4636 AliInfo("No File found!");
4640 for(Int_t k = 0; k < listofgraphs->GetEntriesFast(); k++){
4641 fout->WriteTObject((TObject *) listofgraphs->At(k),((TObject *)listofgraphs->At(k))->GetName()
4642 ,(Option_t *) "OverWrite");
4651 //____________Plot histos DB___________________________________________________
4654 //_____________________________________________________________________________
4655 void AliTRDCalibraFit::PlotCHDB()
4658 // Plot the histos for fDebug = 3 and fDebug = 4 to visualise the detector
4661 TCanvas *cchdb = new TCanvas("cchdb","",50,50,600,800);
4663 if(fFitChargeOn) nb++;
4664 if(fFitChargeBisOn) nb++;
4665 if(fMeanChargeOn) nb++;
4666 if(fFitMeanWOn) nb++;
4668 cchdb->Divide(nb,1);
4672 fCoefChargeDB[1]->Draw("LEGO");
4677 fCoefChargeDB[0]->Draw("LEGO");
4682 fCoefChargeDB[3]->Draw("LEGO");
4685 if(fFitChargeBisOn){
4687 fCoefChargeDB[2]->Draw("LEGO");
4688 //it doesn't matter anymore but....
4694 //_____________________________________________________________________________
4695 void AliTRDCalibraFit::PlotPHDB()
4698 // Plot the histos for fDebug = 3 and fDebug = 4 to visualise the detector
4701 TCanvas *cphdb = new TCanvas("cphdb","",50,50,600,800);
4703 if(fFitPol2On) nb++;
4705 if(fFitLagrPolOn) nb++;
4707 cphdb->Divide(nb,1);
4711 fCoefVdriftDB[0]->Draw("LEGO");
4716 fCoefVdriftDB[1]->Draw("LEGO");
4721 fCoefVdriftDB[2]->Draw("LEGO");
4727 //_____________________________________________________________________________
4728 void AliTRDCalibraFit::PlotT0DB()
4731 // Plot the histos for fDebug = 3 and fDebug = 4 to visualise the detector
4733 TCanvas *ct0db = new TCanvas("ct0db","",50,50,600,800);
4735 if(fFitPol2On) nb++;
4737 if(fFitLagrPolOn) nb++;
4739 ct0db->Divide(nb,1);
4743 fCoefT0DB[0]->Draw("LEGO");
4748 fCoefT0DB[1]->Draw("LEGO");
4753 fCoefT0DB[2]->Draw("LEGO");
4759 //_____________________________________________________________________________
4760 void AliTRDCalibraFit::PlotPRFDB()
4763 // Plot the histos for fDebug = 3 and fDebug = 4 to visualise the detector
4766 TCanvas *cprfdb = new TCanvas("cprfdb","",50,50,600,800);
4771 cprfdb->Divide(nb,1);
4775 fCoefPRFDB[0]->Draw("LEGO");
4780 fCoefPRFDB[1]->Draw("LEGO");
4787 //____________Write DB Histos__________________________________________________
4790 //_____________________________________________________________________________
4791 void AliTRDCalibraFit::WriteCHDB(TFile *fout)
4794 // If wanted, write the debug histos for fDebug = 3 and fDebug = 4
4797 fout->WriteTObject(fCoefChargeDB[0],fCoefChargeDB[0]->GetName(),(Option_t *) "OverWrite");
4800 fout->WriteTObject(fCoefChargeDB[3],fCoefChargeDB[0]->GetName(),(Option_t *) "OverWrite");
4802 if (fMeanChargeOn) {
4803 fout->WriteTObject(fCoefChargeDB[1],fCoefChargeDB[1]->GetName(),(Option_t *) "OverWrite");
4805 if (fFitChargeBisOn ) {
4806 fout->WriteTObject(fCoefChargeDB[2],fCoefChargeDB[2]->GetName(),(Option_t *) "OverWrite");
4811 //_____________________________________________________________________________
4812 void AliTRDCalibraFit::WritePHDB(TFile *fout)
4815 // If wanted, write the debug histos for fDebug = 3 and fDebug = 4
4819 fout->WriteTObject(fCoefVdriftDB[0],fCoefVdriftDB[0]->GetName(),(Option_t *) "OverWrite");
4822 fout->WriteTObject(fCoefVdriftDB[1],fCoefVdriftDB[1]->GetName(),(Option_t *) "OverWrite");
4825 fout->WriteTObject(fCoefVdriftDB[2],fCoefVdriftDB[2]->GetName(),(Option_t *) "OverWrite");
4830 //_____________________________________________________________________________
4831 void AliTRDCalibraFit::WriteT0DB(TFile *fout)
4834 // If wanted, write the debug histos for fDebug = 3 and fDebug = 4
4838 fout->WriteTObject(fCoefT0DB[0],fCoefT0DB[0]->GetName(),(Option_t *) "OverWrite");
4841 fout->WriteTObject(fCoefT0DB[1],fCoefT0DB[1]->GetName(),(Option_t *) "OverWrite");
4844 fout->WriteTObject(fCoefT0DB[2],fCoefT0DB[2]->GetName(),(Option_t *) "OverWrite");
4849 //_____________________________________________________________________________
4850 void AliTRDCalibraFit::WritePRFDB(TFile *fout)
4853 // If wanted, write the debug histos for fDebug = 3 and fDebug = 4
4856 fout->WriteTObject(fCoefPRFDB[0],fCoefPRFDB[0]->GetName(),(Option_t *) "OverWrite");
4859 fout->WriteTObject(fCoefPRFDB[1],fCoefPRFDB[1]->GetName(),(Option_t *) "OverWrite");
4865 //____________Calcul Coef Mean_________________________________________________
4868 //_____________________________________________________________________________
4869 Bool_t AliTRDCalibraFit::CalculT0CoefMean(Int_t dect, Int_t idect)
4872 // For the detector Dect calcul the mean time 0
4873 // for the calibration group idect from the choosen database
4876 AliTRDcalibDB *cal = AliTRDcalibDB::Instance();
4878 AliInfo("Could not get calibDB Manager");
4884 if ((fDebug != 2) && fAccCDB) {
4886 for (Int_t row = fCalibraMode->GetRowMin(1); row < fCalibraMode->GetRowMax(1); row++) {
4887 for (Int_t col = fCalibraMode->GetColMin(1); col < fCalibraMode->GetColMax(1); col++) {
4889 if ((fCalibraMode->GetNz(1) > 0) ||
4890 (fCalibraMode->GetNrphi(1) > 0)) {
4891 fT0Coef[2] += (Float_t) cal->GetT0(dect,col,row);
4895 fT0Coef[2] += (Float_t) cal->GetT0Average(dect);
4900 fT0Coef[2] = fT0Coef[2] / ((fCalibraMode->GetColMax(1)-fCalibraMode->GetColMin(1))
4901 * (fCalibraMode->GetRowMax(1)-fCalibraMode->GetRowMin(1)));
4902 if ((fDebug == 1) ||
4904 fCoefT0[2][idect] = fT0Coef[2];
4913 //_____________________________________________________________________________
4914 Bool_t AliTRDCalibraFit::CalculChargeCoefMean(Int_t dect, Int_t idect, Bool_t vrai)
4917 // For the detector Dect calcul the mean gain factor
4918 // for the calibration group idect from the choosen database
4921 AliTRDcalibDB *cal = AliTRDcalibDB::Instance();
4923 AliInfo("Could not get calibDB Manager");
4926 AliTRDCommonParam *parCom = AliTRDCommonParam::Instance();
4928 AliInfo("Could not get CommonParam Manager");
4932 fChargeCoef[3] = 0.0;
4936 for (Int_t row = fCalibraMode->GetRowMin(0); row < fCalibraMode->GetRowMax(0); row++) {
4937 for (Int_t col = fCalibraMode->GetColMin(0); col < fCalibraMode->GetColMax(0); col++) {
4939 if ((fCalibraMode->GetNz(0) > 0) ||
4940 (fCalibraMode->GetNrphi(0) > 0)) {
4942 fChargeCoef[3] += (Float_t) cal->GetGainFactor(dect,col,row);
4944 if (vrai && fAccCDB) {
4945 fScaleFitFactor += (Float_t) cal->GetGainFactor(dect,col,row);
4948 fChargeCoef[3] += 1.0;
4950 if (vrai && (!fAccCDB)) {
4951 fScaleFitFactor += 1.0;
4957 fChargeCoef[3] += (Float_t) cal->GetGainFactorAverage(dect);
4959 if (vrai && fAccCDB) {
4960 fScaleFitFactor += ((Float_t) cal->GetGainFactorAverage(dect));
4963 fChargeCoef[3] += 1.0;
4965 if (vrai && (!fAccCDB)) {
4966 fScaleFitFactor += 1.0;
4972 fChargeCoef[3] = fChargeCoef[3] / ((fCalibraMode->GetColMax(0)-fCalibraMode->GetColMin(0))
4973 * (fCalibraMode->GetRowMax(0)-fCalibraMode->GetRowMin(0)));
4974 if ((fDebug == 1) ||
4976 fCoefCharge[3][idect]=fChargeCoef[3];
4985 //_____________________________________________________________________________
4986 Bool_t AliTRDCalibraFit::CalculPRFCoefMean(Int_t dect, Int_t idect)
4989 // For the detector Dect calcul the mean sigma of pad response
4990 // function for the calibration group idect from the choosen database
4993 AliTRDcalibDB *cal = AliTRDcalibDB::Instance();
4995 AliInfo("Could not get calibDB Manager");
4999 AliTRDCommonParam *parCom = AliTRDCommonParam::Instance();
5001 AliInfo("Could not get CommonParam Manager");
5010 for (Int_t row = fCalibraMode->GetRowMin(2); row < fCalibraMode->GetRowMax(2); row++) {
5011 for (Int_t col = fCalibraMode->GetColMin(2); col < fCalibraMode->GetColMax(2); col++) {
5012 if ((parCom->GetColMax(GetPlane(dect)) != (col+1)) && (col != 0)) {
5015 fPRFCoef[1] += (Float_t) cal->GetPRFWidth(dect,col,row);
5018 fPRFCoef[1] += GetPRFDefault(GetPlane(dect));
5025 fPRFCoef[1] = fPRFCoef[1]/cot;
5026 if ((fDebug == 1) ||
5028 fCoefPRF[1][idect] = fPRFCoef[1];
5032 if ((fDebug == 1) ||
5035 fCoefPRF[1][idect] = cal->GetPRFWidth(dect,fCalibraMode->GetColMin(2),fCalibraMode->GetRowMin(2));
5038 fCoefPRF[1][idect] = GetPRFDefault(GetPlane(dect));
5049 //_____________________________________________________________________________
5050 Bool_t AliTRDCalibraFit::CalculVdriftCoefMean(Int_t dect, Int_t idect)
5053 // For the detector dect calcul the mean drift velocity for the
5054 // calibration group idect from the choosen database
5057 AliTRDcalibDB *cal = AliTRDcalibDB::Instance();
5059 AliInfo("Could not get calibDB Manager");
5063 fVdriftCoef[2] = 0.0;
5066 for (Int_t row = fCalibraMode->GetRowMin(1); row < fCalibraMode->GetRowMax(1); row++) {
5067 for (Int_t col = fCalibraMode->GetColMin(1); col < fCalibraMode->GetColMax(1); col++) {
5069 if ((fCalibraMode->GetNz(1) > 0) ||
5070 (fCalibraMode->GetNrphi(1) > 0)) {
5072 fVdriftCoef[2] += (Float_t) cal->GetVdrift(dect,col,row);
5075 fVdriftCoef[2] += 1.5;
5081 fVdriftCoef[2] += (Float_t) cal->GetVdriftAverage(dect);
5084 fVdriftCoef[2] += 1.5;
5089 fVdriftCoef[2] = fVdriftCoef[2] / ((fCalibraMode->GetColMax(1)-fCalibraMode->GetColMin(1))
5090 * (fCalibraMode->GetRowMax(1)-fCalibraMode->GetRowMin(1)));
5091 if ((fDebug == 1) ||
5093 fCoefVdrift[2][idect] = fVdriftCoef[2];
5101 //_____________________________________________________________________________
5102 Float_t AliTRDCalibraFit::GetPRFDefault(Int_t plane) const
5105 // Default width of the PRF if there is no database as reference
5132 //____________Fit Methods______________________________________________________
5134 //_____________________________________________________________________________
5135 void AliTRDCalibraFit::FitPente(TH1* projPH, Int_t idect)
5138 // Slope methode for the drift velocity
5142 const Float_t kDrWidth = AliTRDgeometry::DrThick();
5149 Double_t vdriftCoefE = 0.0;
5150 Double_t t0CoefE = 0.0;
5151 fVdriftCoef[1] = 0.0;
5153 TLine *line = new TLine();
5156 TAxis *xpph = projPH->GetXaxis();
5157 Int_t nbins = xpph->GetNbins();
5158 Double_t lowedge = xpph->GetBinLowEdge(1);
5159 Double_t upedge = xpph->GetBinUpEdge(xpph->GetNbins());
5160 Double_t widbins = (upedge - lowedge) / nbins;
5161 Double_t limit = upedge + 0.5 * widbins;
5164 // Beginning of the signal
5165 TH1D *pentea = new TH1D("pentea","pentea",projPH->GetNbinsX(),0,(Float_t) limit);
5166 for (Int_t k = 1; k < projPH->GetNbinsX(); k++) {
5167 pentea->SetBinContent(k,(Double_t) (projPH->GetBinContent(k+1) - projPH->GetBinContent(k)));
5170 binmax = (Int_t) pentea->GetMaximumBin();
5171 if(fDebug == 2) AliInfo(Form("maximum positive bin for the positive slope %d",binmax));
5174 AliInfo("Put the binmax from 1 to 2 to enable the fit");
5176 if (binmax >= nbins) {
5179 AliInfo("Put the binmax from nbins-1 to nbins-2 to enable the fit");
5181 pentea->Fit("pol2","0MR","",TMath::Max(pentea->GetBinCenter(binmax-1),0.0),pentea->GetBinCenter(binmax+1));
5182 Float_t l3P1am = pentea->GetFunction("pol2")->GetParameter(1);
5183 Float_t l3P2am = pentea->GetFunction("pol2")->GetParameter(2);
5184 Float_t l3P1amE = pentea->GetFunction("pol2")->GetParError(1);
5185 Float_t l3P2amE = pentea->GetFunction("pol2")->GetParError(2);
5187 fPhd[0] = -(l3P1am / (2 * l3P2am));
5190 if((l3P1am != 0.0) && (l3P2am != 0.0)){
5191 t0CoefE = (l3P1amE/l3P1am + l3P2amE/l3P2am)*fPhd[0];
5194 if(fDebug == 2) AliInfo(Form("maximum extrapolated positive bin for the positive slope %f",fPhd[0]));
5196 // Amplification region
5199 for (Int_t kbin = 1; kbin < projPH->GetNbinsX(); kbin ++) {
5200 if (((projPH->GetBinContent(kbin+1) - projPH->GetBinContent(kbin)) <= 0.0) &&
5202 (kbin > (fPhd[0]/widbins))) {
5207 if(fDebug == 2) AliInfo(Form("For the amplification region binmax %d",binmax));
5210 AliInfo("Put the binmax from 1 to 2 to enable the fit");
5212 if (binmax >= nbins) {
5215 AliInfo("Put the binmax from nbins-1 to nbins-2 to enable the fit");
5217 projPH->Fit("pol2","0MR","",TMath::Max(projPH->GetBinCenter(binmax-1),0.0),projPH->GetBinCenter(binmax+1));
5218 Float_t l3P1amf = projPH->GetFunction("pol2")->GetParameter(1);
5219 Float_t l3P2amf = projPH->GetFunction("pol2")->GetParameter(2);
5220 Float_t l3P1amfE = projPH->GetFunction("pol2")->GetParError(1);
5221 Float_t l3P2amfE = projPH->GetFunction("pol2")->GetParError(2);
5224 fPhd[1] = -(l3P1amf / (2 * l3P2amf));
5226 if((l3P1amf != 0.0) && (l3P2amf != 0.0)){
5227 vdriftCoefE = (l3P1amfE/l3P1amf + l3P2amfE/l3P2amf)*fPhd[1];
5230 t0CoefE = vdriftCoefE;
5232 if(fDebug == 2) AliInfo(Form("For the amplification region extrapolated binmax %f",fPhd[1]));
5235 TH1D *pente = new TH1D("pente","pente",projPH->GetNbinsX(),0,(Float_t) limit);
5236 for (Int_t k = TMath::Min(binmax+4,projPH->GetNbinsX()); k < projPH->GetNbinsX(); k++) {
5237 pente->SetBinContent(k,(Double_t) (projPH->GetBinContent(k+1) - projPH->GetBinContent(k)));
5240 if(pente->GetEntries() > 0) binmin = (Int_t) pente->GetMinimumBin();
5243 AliInfo("Put the binmax from 1 to 2 to enable the fit");
5245 if (binmin >= nbins) {
5248 AliInfo("Put the binmax from nbins-1 to nbins-2 to enable the fit");
5250 if(fDebug == 2) AliInfo(Form("For the drift region binmin %d",binmin));
5254 ,TMath::Max(pente->GetBinCenter(binmin-1), 0.0)
5255 ,TMath::Min(pente->GetBinCenter(binmin+1),(Double_t) limit));
5256 Float_t l3P1dr = pente->GetFunction("pol2")->GetParameter(1);
5257 Float_t l3P2dr = pente->GetFunction("pol2")->GetParameter(2);
5258 Float_t l3P1drE = pente->GetFunction("pol2")->GetParError(1);
5259 Float_t l3P2drE = pente->GetFunction("pol2")->GetParError(2);
5261 fPhd[2] = -(l3P1dr / (2 * l3P2dr));
5263 if((l3P1dr != 0.0) && (l3P2dr != 0.0)){
5264 vdriftCoefE += (l3P1drE/l3P1dr + l3P2drE/l3P2dr)*fPhd[2];
5266 if(fDebug == 2) AliInfo(Form("For the drift region extrapolated binmax %f",fPhd[2]));
5268 Float_t fPhdt0 = 0.0;
5269 if(fTakeTheMaxPH) fPhdt0 = fPhd[1];
5270 else fPhdt0 = fPhd[0];
5272 if ((fPhd[2] > fPhd[0]) &&
5273 (fPhd[2] > fPhd[1]) &&
5274 (fPhd[1] > fPhd[0]) &&
5276 fVdriftCoef[1] = (kDrWidth) / (fPhd[2]-fPhd[1]);
5277 if(fFitPHNDB == 1) fNumberFitSuccess++;
5278 if (fPhdt0 >= 0.0) {
5279 fT0Coef[1] = (fPhdt0 - fT0Shift) / widbins;
5280 if (fT0Coef[1] < -1.0) {
5281 fT0Coef[1] = fT0Coef[2];
5285 fT0Coef[1] = fT0Coef[2];
5289 fVdriftCoef[1] = -TMath::Abs(fVdriftCoef[2]);
5290 fT0Coef[1] = fT0Coef[2];
5293 if ((fDebug == 1) ||
5295 fCoefVdrift[1][idect] = fVdriftCoef[1];
5296 fCoefVdriftE[1] [idect] = vdriftCoefE;
5297 fCoefT0[1][idect] = fT0Coef[1];
5298 fCoefT0E[1][idect] = t0CoefE;
5302 TCanvas *cpentei = new TCanvas("cpentei","cpentei",50,50,600,800);
5305 line->SetLineColor(2);
5306 line->DrawLine(fPhd[0],0,fPhd[0],projPH->GetMaximum());
5307 line->DrawLine(fPhd[1],0,fPhd[1],projPH->GetMaximum());
5308 line->DrawLine(fPhd[2],0,fPhd[2],projPH->GetMaximum());
5309 AliInfo(Form("fPhd[0] (beginning of the signal): %f" ,(Float_t) fPhd[0]));
5310 AliInfo(Form("fPhd[1] (end of the amplification region): %f" ,(Float_t) fPhd[1]));
5311 AliInfo(Form("fPhd[2] (end of the drift region): %f" ,(Float_t) fPhd[2]));
5312 AliInfo(Form("fVriftCoef[1] (with only the drift region(default)): %f",(Float_t) fVdriftCoef[1]));
5313 TCanvas *cpentei2 = new TCanvas("cpentei2","cpentei2",50,50,600,800);
5316 TCanvas *cpentei3 = new TCanvas("cpentei3","cpentei3",50,50,600,800);
5330 //_____________________________________________________________________________
5331 void AliTRDCalibraFit::FitLagrangePoly(TH1* projPH, Int_t idect)
5334 // Slope methode but with polynomes de Lagrange
5338 const Float_t kDrWidth = AliTRDgeometry::DrThick();
5341 Double_t *x = new Double_t[5];
5342 Double_t *y = new Double_t[5];
5357 Double_t vdriftCoefE = 0.0;
5358 Double_t t0CoefE = 1.0;
5359 fVdriftCoef[3] = 0.0;
5361 TLine *line = new TLine();
5362 TF1 * polynome = 0x0;
5363 TF1 * polynomea = 0x0;
5364 TF1 * polynomeb = 0x0;
5368 TAxis *xpph = projPH->GetXaxis();
5369 Int_t nbins = xpph->GetNbins();
5370 Double_t lowedge = xpph->GetBinLowEdge(1);
5371 Double_t upedge = xpph->GetBinUpEdge(xpph->GetNbins());
5372 Double_t widbins = (upedge - lowedge) / nbins;
5373 Double_t limit = upedge + 0.5 * widbins;
5378 // Beginning of the signal
5379 TH1D *pentea = new TH1D("pentea","pentea",projPH->GetNbinsX(),0,(Float_t) limit);
5380 for (Int_t k = 1; k < projPH->GetNbinsX(); k++) {
5381 pentea->SetBinContent(k,(Double_t) (projPH->GetBinContent(k+1) - projPH->GetBinContent(k)));
5384 binmax = (Int_t) pentea->GetMaximumBin();
5385 if(fDebug == 2) AliInfo(Form("maximum positive bin for the positive slope %d",binmax));
5388 Double_t minnn = 0.0;
5389 Double_t maxxx = 0.0;
5391 //Determination of minnn and maxxx
5392 //case binmax = nbins -1
5394 if(binmax == (nbins-1)){
5395 minnn = pentea->GetBinCenter(binmax-2);
5396 maxxx = pentea->GetBinCenter(binmax);
5397 x[0] = pentea->GetBinCenter(binmax-2);
5398 x[1] = pentea->GetBinCenter(binmax-1);
5399 x[2] = pentea->GetBinCenter(binmax);
5400 y[0] = pentea->GetBinContent(binmax-2);
5401 y[1] = pentea->GetBinContent(binmax-1);
5402 y[2] = pentea->GetBinContent(binmax);
5403 //Calcul the polynome de Lagrange
5404 c = CalculPolynomeLagrange2(x,y);
5405 AliInfo("At the limit for beginning!");
5407 //case binmax = nbins-2
5409 if(binmax == (nbins-2)){
5410 minnn = pentea->GetBinCenter(binmax-2);
5411 maxxx = pentea->GetBinCenter(binmax+1);
5412 x[0] = pentea->GetBinCenter(binmax-2);
5413 x[1] = pentea->GetBinCenter(binmax-1);
5414 x[2] = pentea->GetBinCenter(binmax);
5415 x[3] = pentea->GetBinCenter(binmax+1);
5416 y[0] = pentea->GetBinContent(binmax-2);
5417 y[1] = pentea->GetBinContent(binmax-1);
5418 y[2] = pentea->GetBinContent(binmax);
5419 y[3] = pentea->GetBinContent(binmax+1);
5420 //Calcul the polynome de Lagrange
5421 c = CalculPolynomeLagrange3(x,y);
5423 //case binmax <= nbins-3
5425 if(binmax <= (nbins-3)){
5426 if((binmax-2) >= 1){
5427 minnn = pentea->GetBinCenter(binmax-2);
5428 maxxx = pentea->GetBinCenter(binmax+2);
5429 x[0] = pentea->GetBinCenter(binmax-2);
5430 x[1] = pentea->GetBinCenter(binmax-1);
5431 x[2] = pentea->GetBinCenter(binmax);
5432 x[3] = pentea->GetBinCenter(binmax+1);
5433 x[4] = pentea->GetBinCenter(binmax+2);
5434 y[0] = pentea->GetBinContent(binmax-2);
5435 y[1] = pentea->GetBinContent(binmax-1);
5436 y[2] = pentea->GetBinContent(binmax);
5437 y[3] = pentea->GetBinContent(binmax+1);
5438 y[4] = pentea->GetBinContent(binmax+2);
5439 //Calcul the polynome de Lagrange
5440 c = CalculPolynomeLagrange4(x,y);
5443 if((binmax-1) == 1){
5444 minnn = pentea->GetBinCenter(binmax-1);
5445 maxxx = pentea->GetBinCenter(binmax+2);
5446 x[0] = pentea->GetBinCenter(binmax-1);
5447 x[1] = pentea->GetBinCenter(binmax);
5448 x[2] = pentea->GetBinCenter(binmax+1);
5449 x[3] = pentea->GetBinCenter(binmax+2);
5450 y[0] = pentea->GetBinContent(binmax-1);
5451 y[1] = pentea->GetBinContent(binmax);
5452 y[2] = pentea->GetBinContent(binmax+1);
5453 y[3] = pentea->GetBinContent(binmax+2);
5454 //Calcul the polynome de Lagrange
5455 c = CalculPolynomeLagrange3(x,y);
5459 minnn = pentea->GetBinCenter(binmax);
5460 maxxx = pentea->GetBinCenter(binmax+2);
5461 x[0] = pentea->GetBinCenter(binmax);
5462 x[1] = pentea->GetBinCenter(binmax+1);
5463 x[2] = pentea->GetBinCenter(binmax+2);
5464 y[0] = pentea->GetBinContent(binmax);
5465 y[1] = pentea->GetBinContent(binmax+1);
5466 y[2] = pentea->GetBinContent(binmax+2);
5467 //Calcul the polynome de Lagrange
5468 c = CalculPolynomeLagrange2(x,y);
5471 //pass but should not happen
5472 if((binmax <= (nbins-3)) && (binmax < 1)){
5476 if(fDebug == 2) AliInfo(Form("For the beginning region binmax %d",binmax));
5479 polynomeb = new TF1("polb","[0]+[1]*x+[2]*x*x+[3]*x*x*x+[4]*x*x*x*x",minnn,maxxx);
5480 polynomeb->SetParameters(c[0],c[1],c[2],c[3],c[4]);
5482 AliInfo(Form("for the beginning: c[0] %f, c[1] %f, c[2] %f, c[3] %f, c[4] %f",c[0],c[1],c[2],c[3],c[4]));
5485 Double_t step = (maxxx-minnn)/10000;
5487 Double_t maxvalue = 0.0;
5488 Double_t placemaximum = minnn;
5489 for(Int_t o = 0; o < 10000; o++){
5490 if(o == 0) maxvalue = polynomeb->Eval(l);
5491 if(maxvalue < (polynomeb->Eval(l))){
5492 maxvalue = polynomeb->Eval(l);
5497 fPhd[0] = placemaximum;
5500 if(fDebug == 2) AliInfo(Form("maximum extrapolated positive bin for the positive slope %f",fPhd[0]));
5502 // Amplification region
5505 for (Int_t kbin = 1; kbin < projPH->GetNbinsX(); kbin ++) {
5506 if (((projPH->GetBinContent(kbin+1) - projPH->GetBinContent(kbin)) <= 0.0) &&
5508 (kbin > (fPhd[0]/widbins))) {
5513 if(fDebug == 2) AliInfo(Form("For the amplification region binmax %d",binmax));
5515 Double_t minn = 0.0;
5516 Double_t maxx = 0.0;
5518 //Determination of minn and maxx
5519 //case binmax = nbins
5521 if(binmax == nbins){
5522 minn = projPH->GetBinCenter(binmax-2);
5523 maxx = projPH->GetBinCenter(binmax);
5524 x[0] = projPH->GetBinCenter(binmax-2);
5525 x[1] = projPH->GetBinCenter(binmax-1);
5526 x[2] = projPH->GetBinCenter(binmax);
5527 y[0] = projPH->GetBinContent(binmax-2);
5528 y[1] = projPH->GetBinContent(binmax-1);
5529 y[2] = projPH->GetBinContent(binmax);
5530 //Calcul the polynome de Lagrange
5531 c = CalculPolynomeLagrange2(x,y);
5532 AliInfo("At the limit for the drift!");
5534 //case binmax = nbins-1
5536 if(binmax == (nbins-1)){
5537 minn = projPH->GetBinCenter(binmax-2);
5538 maxx = projPH->GetBinCenter(binmax+1);
5539 x[0] = projPH->GetBinCenter(binmax-2);
5540 x[1] = projPH->GetBinCenter(binmax-1);
5541 x[2] = projPH->GetBinCenter(binmax);
5542 x[3] = projPH->GetBinCenter(binmax+1);
5543 y[0] = projPH->GetBinContent(binmax-2);
5544 y[1] = projPH->GetBinContent(binmax-1);
5545 y[2] = projPH->GetBinContent(binmax);
5546 y[3] = projPH->GetBinContent(binmax+1);
5547 //Calcul the polynome de Lagrange
5548 c = CalculPolynomeLagrange3(x,y);
5550 //case binmax <= nbins-2
5552 if(binmax <= (nbins-2)){
5553 if((binmax-2) >= 1){
5554 minn = projPH->GetBinCenter(binmax-2);
5555 maxx = projPH->GetBinCenter(binmax+2);
5556 x[0] = projPH->GetBinCenter(binmax-2);
5557 x[1] = projPH->GetBinCenter(binmax-1);
5558 x[2] = projPH->GetBinCenter(binmax);
5559 x[3] = projPH->GetBinCenter(binmax+1);
5560 x[4] = projPH->GetBinCenter(binmax+2);
5561 y[0] = projPH->GetBinContent(binmax-2);
5562 y[1] = projPH->GetBinContent(binmax-1);
5563 y[2] = projPH->GetBinContent(binmax);
5564 y[3] = projPH->GetBinContent(binmax+1);
5565 y[4] = projPH->GetBinContent(binmax+2);
5566 //Calcul the polynome de Lagrange
5567 c = CalculPolynomeLagrange4(x,y);
5570 if((binmax-1) == 1){
5571 minn = projPH->GetBinCenter(binmax-1);
5572 maxx = projPH->GetBinCenter(binmax+2);
5573 x[0] = projPH->GetBinCenter(binmax-1);
5574 x[1] = projPH->GetBinCenter(binmax);
5575 x[2] = projPH->GetBinCenter(binmax+1);
5576 x[3] = projPH->GetBinCenter(binmax+2);
5577 y[0] = projPH->GetBinContent(binmax-1);
5578 y[1] = projPH->GetBinContent(binmax);
5579 y[2] = projPH->GetBinContent(binmax+1);
5580 y[3] = projPH->GetBinContent(binmax+2);
5581 //Calcul the polynome de Lagrange
5582 c = CalculPolynomeLagrange3(x,y);
5586 minn = projPH->GetBinCenter(binmax);
5587 maxx = projPH->GetBinCenter(binmax+2);
5588 x[0] = projPH->GetBinCenter(binmax);
5589 x[1] = projPH->GetBinCenter(binmax+1);
5590 x[2] = projPH->GetBinCenter(binmax+2);
5591 y[0] = projPH->GetBinContent(binmax);
5592 y[1] = projPH->GetBinContent(binmax+1);
5593 y[2] = projPH->GetBinContent(binmax+2);
5594 //Calcul the polynome de Lagrange
5595 c = CalculPolynomeLagrange2(x,y);
5598 //pass but should not happen
5599 if((binmax <= (nbins-2)) && (binmax < 1)){
5603 if(fDebug == 2) AliInfo(Form("For the amplification region binmax %d",binmax));
5606 polynomea = new TF1("pola","[0]+[1]*x+[2]*x*x+[3]*x*x*x+[4]*x*x*x*x",minn,maxx);
5607 polynomea->SetParameters(c[0],c[1],c[2],c[3],c[4]);
5609 AliInfo(Form("for the amplification: c[0] %f, c[1] %f, c[2] %f, c[3] %f, c[4] %f",c[0],c[1],c[2],c[3],c[4]));
5612 Double_t step = (maxx-minn)/1000;
5614 Double_t maxvalue = 0.0;
5615 Double_t placemaximum = minn;
5616 for(Int_t o = 0; o < 1000; o++){
5617 if(o == 0) maxvalue = polynomea->Eval(l);
5618 if(maxvalue < (polynomea->Eval(l))){
5619 maxvalue = polynomea->Eval(l);
5624 fPhd[1] = placemaximum;
5627 if(fDebug == 2) AliInfo(Form("For the amplification region extrapolated binmax %f",fPhd[1]));
5630 TH1D *pente = new TH1D("pente","pente", projPH->GetNbinsX(),0,(Float_t) limit);
5631 for (Int_t k = TMath::Min(binmax+4, projPH->GetNbinsX()); k < projPH->GetNbinsX(); k++) {
5632 pente->SetBinContent(k,(Double_t) (projPH->GetBinContent(k+1) - projPH->GetBinContent(k)));
5635 if(pente->GetEntries() > 0) binmin = (Int_t) pente->GetMinimumBin();
5641 AliInfo("Put the binmax from 1 to 2 to enable the fit");
5645 if((projPH->GetBinContent(binmin)-projPH->GetBinError(binmin)) < (projPH->GetBinContent(binmin+1))) put = kFALSE;
5646 if((projPH->GetBinContent(binmin)+projPH->GetBinError(binmin)) > (projPH->GetBinContent(binmin-1))) put = kFALSE;
5649 AliInfo(Form("binmin %d BinContent %f BinError %f",binmin
5650 ,projPH->GetBinContent(binmin)
5651 ,projPH->GetBinError(binmin)));
5652 AliInfo(Form("binmin-1 %d BinContent %f BinError %f",binmin-1
5653 ,projPH->GetBinContent(binmin-1)
5654 ,projPH->GetBinError(binmin-1)));
5655 AliInfo(Form("binmin+1 %d BinContent %f BinError %f",binmin+1
5656 ,projPH->GetBinContent(binmin+1)
5657 ,projPH->GetBinError(binmin+1)));
5662 Bool_t case1 = kFALSE;
5663 Bool_t case2 = kFALSE;
5664 Bool_t case4 = kFALSE;
5666 //Determination of min and max
5667 //case binmin <= nbins-3
5669 if((binmin <= (nbins-3)) && ((binmin-2) >= TMath::Min(binmax+4, projPH->GetNbinsX()))){
5670 min = pente->GetBinCenter(binmin-2);
5671 max = pente->GetBinCenter(binmin+2);
5672 x[0] = pente->GetBinCenter(binmin-2);
5673 x[1] = pente->GetBinCenter(binmin-1);
5674 x[2] = pente->GetBinCenter(binmin);
5675 x[3] = pente->GetBinCenter(binmin+1);
5676 x[4] = pente->GetBinCenter(binmin+2);
5677 y[0] = pente->GetBinContent(binmin-2);
5678 y[1] = pente->GetBinContent(binmin-1);
5679 y[2] = pente->GetBinContent(binmin);
5680 y[3] = pente->GetBinContent(binmin+1);
5681 y[4] = pente->GetBinContent(binmin+2);
5682 //Calcul the polynome de Lagrange
5683 c = CalculPolynomeLagrange4(x,y);
5685 if((pente->GetBinContent(binmin+2) <= pente->GetBinContent(binmin+1)) &&
5686 (pente->GetBinContent(binmin-2) <= pente->GetBinContent(binmin-1))) put = kFALSE;
5687 if(((binmin+3) <= (nbins-1)) &&
5688 (pente->GetBinContent(binmin+3) <= pente->GetBinContent(binmin+2)) &&
5689 ((binmin-3) >= TMath::Min(binmax+4, projPH->GetNbinsX())) &&
5690 (pente->GetBinContent(binmin-3) <= pente->GetBinContent(binmin-2))) put = kFALSE;
5691 if((pente->GetBinContent(binmin+2) <= pente->GetBinContent(binmin+1)) &&
5692 (pente->GetBinContent(binmin-2) > pente->GetBinContent(binmin-1))) case1 = kTRUE;
5693 if((pente->GetBinContent(binmin+2) > pente->GetBinContent(binmin+1)) &&
5694 (pente->GetBinContent(binmin-2) <= pente->GetBinContent(binmin-1))) case4 = kTRUE;
5696 //case binmin = nbins-2
5698 if(((binmin == (nbins-2)) && ((binmin-2) >= TMath::Min(binmax+4, projPH->GetNbinsX()))) ||
5700 min = pente->GetBinCenter(binmin-2);
5701 max = pente->GetBinCenter(binmin+1);
5702 x[0] = pente->GetBinCenter(binmin-2);
5703 x[1] = pente->GetBinCenter(binmin-1);
5704 x[2] = pente->GetBinCenter(binmin);
5705 x[3] = pente->GetBinCenter(binmin+1);
5706 y[0] = pente->GetBinContent(binmin-2);
5707 y[1] = pente->GetBinContent(binmin-1);
5708 y[2] = pente->GetBinContent(binmin);
5709 y[3] = pente->GetBinContent(binmin+1);
5710 //Calcul the polynome de Lagrange
5711 c = CalculPolynomeLagrange3(x,y);
5712 //richtung +: nothing
5714 if((pente->GetBinContent(binmin-2) <= pente->GetBinContent(binmin-1))) case2 = kTRUE;
5717 if(((binmin <= (nbins-3)) && ((binmin-1) == TMath::Min(binmax+4, projPH->GetNbinsX()))) ||
5719 min = pente->GetBinCenter(binmin-1);
5720 max = pente->GetBinCenter(binmin+2);
5721 x[0] = pente->GetBinCenter(binmin-1);
5722 x[1] = pente->GetBinCenter(binmin);
5723 x[2] = pente->GetBinCenter(binmin+1);
5724 x[3] = pente->GetBinCenter(binmin+2);
5725 y[0] = pente->GetBinContent(binmin-1);
5726 y[1] = pente->GetBinContent(binmin);
5727 y[2] = pente->GetBinContent(binmin+1);
5728 y[3] = pente->GetBinContent(binmin+2);
5729 //Calcul the polynome de Lagrange
5730 c = CalculPolynomeLagrange3(x,y);
5732 if((pente->GetBinContent(binmin+2) <= pente->GetBinContent(binmin+1))) case2 = kTRUE;
5735 if((binmin <= (nbins-3)) && (binmin == TMath::Min(binmax+4, projPH->GetNbinsX()))){
5736 min = pente->GetBinCenter(binmin);
5737 max = pente->GetBinCenter(binmin+2);
5738 x[0] = pente->GetBinCenter(binmin);
5739 x[1] = pente->GetBinCenter(binmin+1);
5740 x[2] = pente->GetBinCenter(binmin+2);
5741 y[0] = pente->GetBinContent(binmin);
5742 y[1] = pente->GetBinContent(binmin+1);
5743 y[2] = pente->GetBinContent(binmin+2);
5744 //Calcul the polynome de Lagrange
5745 c = CalculPolynomeLagrange2(x,y);
5747 if((pente->GetBinContent(binmin+2) <= pente->GetBinContent(binmin+1))) put = kFALSE;
5750 if(((binmin == (nbins-2)) && ((binmin-1) == TMath::Min(binmax+4, projPH->GetNbinsX()))) ||
5752 min = pente->GetBinCenter(binmin-1);
5753 max = pente->GetBinCenter(binmin+1);
5754 x[0] = pente->GetBinCenter(binmin-1);
5755 x[1] = pente->GetBinCenter(binmin);
5756 x[2] = pente->GetBinCenter(binmin+1);
5757 y[0] = pente->GetBinContent(binmin-1);
5758 y[1] = pente->GetBinContent(binmin);
5759 y[2] = pente->GetBinContent(binmin+1);
5760 //Calcul the polynome de Lagrange
5761 c = CalculPolynomeLagrange2(x,y);
5762 //richtung +: nothing
5763 //richtung -: nothing
5765 //case binmin = nbins-1
5767 if((binmin == (nbins-1)) && ((binmin-2) >= TMath::Min(binmax+4, projPH->GetNbinsX()))){
5768 min = pente->GetBinCenter(binmin-2);
5769 max = pente->GetBinCenter(binmin);
5770 x[0] = pente->GetBinCenter(binmin-2);
5771 x[1] = pente->GetBinCenter(binmin-1);
5772 x[2] = pente->GetBinCenter(binmin);
5773 y[0] = pente->GetBinContent(binmin-2);
5774 y[1] = pente->GetBinContent(binmin-1);
5775 y[2] = pente->GetBinContent(binmin);
5776 //Calcul the polynome de Lagrange
5777 c = CalculPolynomeLagrange2(x,y);
5778 AliInfo("At the limit for the drift!");
5779 //fluctuation too big!
5780 //richtung +: nothing
5782 if((pente->GetBinContent(binmin-2) <= pente->GetBinContent(binmin-1))) put = kFALSE;
5784 if((binmin == (nbins-1)) && ((binmin-2) < TMath::Min(binmax+4, projPH->GetNbinsX()))) {
5786 AliInfo("At the limit for the drift and not usable!");
5790 if((binmin == (nbins-2)) && ((binmin-1) < TMath::Min(binmax+4, projPH->GetNbinsX()))){
5792 AliInfo("For the drift...problem!");
5795 //pass but should not happen
5796 if((binmin <= (nbins-3)) && (binmin < TMath::Min(binmax+4, projPH->GetNbinsX()))){
5798 AliInfo("For the drift...problem!");
5801 if(fDebug == 2) AliInfo(Form("For the drift region binmax %d",binmin));
5804 polynome = new TF1("pol","[0]+[1]*x+[2]*x*x+[3]*x*x*x+[4]*x*x*x*x",min,max);
5805 polynome->SetParameters(c[0],c[1],c[2],c[3],c[4]);
5807 AliInfo(Form("c[0] %f, c[1] %f, c[2] %f, c[3] %f, c[4] %f",c[0],c[1],c[2],c[3],c[4]));
5809 //AliInfo(Form("GetMinimum of the function %f",polynome->GetMinimumX()));
5810 Double_t step = (max-min)/1000;
5812 Double_t minvalue = 0.0;
5813 Double_t placeminimum = min;
5814 for(Int_t o = 0; o < 1000; o++){
5815 if(o == 0) minvalue = polynome->Eval(l);
5816 if(minvalue > (polynome->Eval(l))){
5817 minvalue = polynome->Eval(l);
5822 fPhd[2] = placeminimum;
5825 if(fDebug == 2) AliInfo(Form("For the drift region extrapolated binmax %f",fPhd[2]));
5827 Float_t fPhdt0 = 0.0;
5828 if(fTakeTheMaxPH) fPhdt0 = fPhd[1];
5829 else fPhdt0 = fPhd[0];
5831 if ((fPhd[2] > fPhd[0]) &&
5832 (fPhd[2] > fPhd[1]) &&
5833 (fPhd[1] > fPhd[0]) &&
5835 fVdriftCoef[3] = (kDrWidth) / (fPhd[2]-fPhd[1]);
5836 if(fFitPHNDB == 3) fNumberFitSuccess++;
5837 if (fPhdt0 >= 0.0) {
5838 fT0Coef[3] = (fPhdt0 - fT0Shift) / widbins;
5839 if (fT0Coef[3] < -1.0) {
5840 fT0Coef[3] = fT0Coef[2];
5844 fT0Coef[3] = fT0Coef[2];
5848 fVdriftCoef[3] = -TMath::Abs(fVdriftCoef[2]);
5849 fT0Coef[3] = fT0Coef[2];
5852 if ((fDebug == 1) ||
5854 fCoefVdrift[3][idect] = fVdriftCoef[3];
5855 fCoefVdriftE[2] [idect] = vdriftCoefE;
5856 fCoefT0[3][idect] = fT0Coef[3];
5857 fCoefT0E[2][idect] = t0CoefE;
5861 TCanvas *cpentei = new TCanvas("cpentei","cpentei",50,50,600,800);
5864 line->SetLineColor(2);
5865 line->DrawLine(fPhd[0],0,fPhd[0],projPH->GetMaximum());
5866 line->DrawLine(fPhd[1],0,fPhd[1],projPH->GetMaximum());
5867 line->DrawLine(fPhd[2],0,fPhd[2],projPH->GetMaximum());
5868 AliInfo(Form("fPhd[0] (beginning of the signal): %f" ,(Float_t) fPhd[0]));
5869 AliInfo(Form("fPhd[1] (end of the amplification region): %f" ,(Float_t) fPhd[1]));
5870 AliInfo(Form("fPhd[2] (end of the drift region): %f" ,(Float_t) fPhd[2]));
5871 AliInfo(Form("fVriftCoef[3] (with only the drift region(default)): %f",(Float_t) fVdriftCoef[3]));
5872 TCanvas *cpentei2 = new TCanvas("cpentei2","cpentei2",50,50,600,800);
5875 TCanvas *cpentei3 = new TCanvas("cpentei3","cpentei3",50,50,600,800);
5888 projPH->SetDirectory(0);
5892 //_____________________________________________________________________________
5893 void AliTRDCalibraFit::FitPH(TH1* projPH, Int_t idect)
5896 // Fit methode for the drift velocity
5900 const Float_t kDrWidth = AliTRDgeometry::DrThick();
5903 TAxis *xpph = projPH->GetXaxis();
5904 Double_t upedge = xpph->GetBinUpEdge(xpph->GetNbins());
5906 TF1 *fPH = new TF1("fPH",AliTRDCalibraFit::PH,-0.05,3.2,6);
5907 fPH->SetParameter(0,0.469); // Scaling
5908 fPH->SetParameter(1,0.18); // Start
5909 fPH->SetParameter(2,0.0857325); // AR
5910 fPH->SetParameter(3,1.89); // DR
5911 fPH->SetParameter(4,0.08); // QA/QD
5912 fPH->SetParameter(5,0.0); // Baseline
5914 TLine *line = new TLine();
5916 fVdriftCoef[0] = 0.0;
5918 Double_t vdriftCoefE = 0.0;
5919 Double_t t0CoefE = 0.0;
5921 if (idect%fFitPHPeriode == 0) {
5923 AliInfo(Form("<AliTRDCalibraFit::FitPH> The detector %d will be fitted",idect));
5924 fPH->SetParameter(0,(projPH->Integral()-(projPH->GetBinContent(1)*projPH->GetNbinsX())) * 0.00028); // Scaling
5925 fPH->SetParameter(1,fPhd[0] - 0.1); // Start
5926 fPH->SetParameter(2,fPhd[1] - fPhd[0]); // AR
5927 fPH->SetParameter(3,fPhd[2] - fPhd[1]); // DR
5928 fPH->SetParameter(4,0.225); // QA/QD
5929 fPH->SetParameter(5,(Float_t) projPH->GetBinContent(1));
5932 projPH->Fit(fPH,"0M","",0.0,upedge);
5936 TCanvas *cpente = new TCanvas("cpente","cpente",50,50,600,800);
5938 projPH->Fit(fPH,"M+","",0.0,upedge);
5940 line->SetLineColor(4);
5941 line->DrawLine(fPH->GetParameter(1)
5943 ,fPH->GetParameter(1)
5944 ,projPH->GetMaximum());
5945 line->DrawLine(fPH->GetParameter(1)+fPH->GetParameter(2)
5947 ,fPH->GetParameter(1)+fPH->GetParameter(2)
5948 ,projPH->GetMaximum());
5949 line->DrawLine(fPH->GetParameter(1)+fPH->GetParameter(2)+fPH->GetParameter(3)
5951 ,fPH->GetParameter(1)+fPH->GetParameter(2)+fPH->GetParameter(3)
5952 ,projPH->GetMaximum());
5955 if (fPH->GetParameter(3) != 0) {
5956 if(fFitPHNDB == 0) fNumberFitSuccess++;
5957 fVdriftCoef[0] = kDrWidth / (fPH->GetParameter(3));
5958 vdriftCoefE = (fPH->GetParError(3)/fPH->GetParameter(3))*fVdriftCoef[0];
5959 fT0Coef[0] = fPH->GetParameter(1);
5960 t0CoefE = fPH->GetParError(1);
5963 fVdriftCoef[0] = -TMath::Abs(fVdriftCoef[2]);
5964 fT0Coef[0] = fT0Coef[2];
5967 if ((fDebug == 1) ||
5969 fCoefVdrift[0][idect] = fVdriftCoef[0];
5970 fCoefVdriftE[0][idect] = vdriftCoefE;
5971 fCoefT0[0][idect] = fT0Coef[0];
5972 fCoefT0E[0][idect] = t0CoefE;
5975 AliInfo(Form("fVdriftCoef[0]: %f",(Float_t) fVdriftCoef[0]));
5982 // Put the default value
5983 if ((fDebug <= 1) ||
5985 fCoefVdrift[0][idect] = -TMath::Abs(fVdriftCoef[2]);
5986 fCoefT0[0][idect] = -TMath::Abs(fT0Coef[2]);
5997 //_____________________________________________________________________________
5998 void AliTRDCalibraFit::FitPRF(TH1 *projPRF, Int_t idect)
6001 // Fit methode for the sigma of the pad response function
6005 Double_t prfCoefE = 0.0;
6009 projPRF->Fit("gaus","0M","",-fRangeFitPRF,fRangeFitPRF);
6014 TCanvas *cfit = new TCanvas("cfit","cfit",50,50,600,800);
6016 projPRF->Fit("gaus","M+","",-fRangeFitPRF,fRangeFitPRF);
6021 fPRFCoef[0] = projPRF->GetFunction("gaus")->GetParameter(2);
6022 prfCoefE = projPRF->GetFunction("gaus")->GetParError(2);
6024 if(fPRFCoef[0] <= 0.0) fPRFCoef[0] = -fPRFCoef[1];
6026 if(fFitPRFNDB == 0) fNumberFitSuccess++;
6029 if ((fDebug == 1) ||
6031 fCoefPRF[0][idect] = fPRFCoef[0];
6032 fCoefPRFE[0][idect] = prfCoefE;
6035 AliInfo(Form("fPRFCoef[0]: %f",(Float_t) fPRFCoef[0]));
6040 //_____________________________________________________________________________
6041 void AliTRDCalibraFit::RmsPRF(TH1 *projPRF, Int_t idect)
6044 // Fit methode for the sigma of the pad response function
6048 Double_t prfCoefE = 0.0;
6052 TCanvas *cfit = new TCanvas("cfit","cfit",50,50,600,800);
6058 fPRFCoef[2] = projPRF->GetRMS();
6060 if(fPRFCoef[2] <= 0.0) fPRFCoef[2] = -fPRFCoef[1];
6062 if(fFitPRFNDB == 2) fNumberFitSuccess++;
6065 if ((fDebug == 1) ||
6067 fCoefPRF[2][idect] = fPRFCoef[2];
6068 fCoefPRFE[1][idect] = prfCoefE;
6071 AliInfo(Form("fPRFCoef[2]: %f",(Float_t) fPRFCoef[2]));
6076 //_____________________________________________________________________________
6077 void AliTRDCalibraFit::FitMean(TH1 *projch, Int_t idect, Double_t nentries)
6080 // Only mean methode for the gain factor
6083 Double_t chargeCoefE1 = 0.0;
6084 if(nentries > 0) chargeCoefE1 = projch->GetRMS()/TMath::Sqrt(nentries);
6087 TCanvas *cpmean = new TCanvas("cpmean","cpmean",50,50,600,800);
6092 if(fFitChargeNDB == 1){
6093 CalculChargeCoefMean(fCountDet[0],(Int_t) idect,kTRUE);
6094 fNumberFitSuccess++;
6096 if ((fDebug == 1) ||
6098 fCoefCharge[1][idect]= fChargeCoef[1];
6099 fCoefChargeE[1][idect]= chargeCoefE1;
6103 //_____________________________________________________________________________
6104 void AliTRDCalibraFit::FitMeanW(TH1 *projch, Int_t idect)
6107 // mean w methode for the gain factor
6111 Int_t nybins = projch->GetNbinsX();
6113 //The weight function
6114 Double_t a = 0.00228515;
6115 Double_t b = -0.00231487;
6116 Double_t c = 0.00044298;
6117 Double_t d = -0.00379239;
6118 Double_t e = 0.00338349;
6126 //A arbitrary error for the moment
6127 Double_t chargeCoefE4 = 0.0;
6128 fChargeCoef[4] = 0.0;
6131 Double_t sumw = 0.0;
6133 Int_t sumAll = (Int_t) projch->GetEntries();
6134 Int_t sumCurrent = 0;
6135 for(Int_t k = 0; k <nybins; k++){
6136 Double_t fraction = Float_t(sumCurrent)/Float_t(sumAll);
6137 if (fraction>0.95) break;
6138 Double_t weight = a + b*fraction + c*fraction*fraction + d *fraction*fraction*fraction+
6139 e*fraction*fraction*fraction*fraction;
6140 sumw += weight*projch->GetBinContent(k+1)*projch->GetBinCenter(k+1);
6141 sum += weight*projch->GetBinContent(k+1);
6142 sumCurrent += (Int_t) projch->GetBinContent(k+1);
6143 //printf("fraction %f, weight %f, bincontent %f\n",fraction,weight,projch->GetBinContent(k+1));
6145 if(sum > 0.0) fChargeCoef[4] = (sumw/sum);
6148 AliInfo(Form("fChargeCoef[4] is %f for the dect %d",fChargeCoef[4],idect));
6149 TCanvas *cpmeanw = new TCanvas("cpmeanw","cpmeanw",50,50,600,800);
6154 if(fFitChargeNDB == 4){
6155 fNumberFitSuccess++;
6156 CalculChargeCoefMean(fCountDet[0],(Int_t) idect,kTRUE);
6158 if ((fDebug == 1) ||
6160 fCoefCharge[4][idect]= fChargeCoef[4];
6161 fCoefChargeE[3][idect]= chargeCoefE4;
6166 //_____________________________________________________________________________
6167 void AliTRDCalibraFit::FitCH(TH1 *projch, Int_t idect)
6170 // Fit methode for the gain factor
6173 fChargeCoef[0] = 0.0;
6174 Double_t chargeCoefE0 = 0.0;
6175 Double_t chisqrl = 0.0;
6176 Double_t chisqrg = 0.0;
6177 Double_t chisqr = 0.0;
6178 TF1 *fLandauGaus = new TF1("fLandauGaus",FuncLandauGaus,0,300,5);
6180 projch->Fit("landau","0",""
6181 ,(Float_t) fChargeCoef[1]/fBeginFitCharge
6182 ,projch->GetBinCenter(projch->GetNbinsX()));
6183 Double_t l3P0 = projch->GetFunction("landau")->GetParameter(0);
6184 Double_t l3P1 = projch->GetFunction("landau")->GetParameter(1);
6185 Double_t l3P2 = projch->GetFunction("landau")->GetParameter(2);
6186 chisqrl = projch->GetFunction("landau")->GetChisquare();
6188 projch->Fit("gaus","0",""
6189 ,(Float_t) fChargeCoef[1]/fBeginFitCharge
6190 ,projch->GetBinCenter(projch->GetNbinsX()));
6191 Double_t g3P0 = projch->GetFunction("gaus")->GetParameter(0);
6192 Double_t g3P2 = projch->GetFunction("gaus")->GetParameter(2);
6193 chisqrg = projch->GetFunction("gaus")->GetChisquare();
6195 fLandauGaus->SetParameters(l3P0,l3P1,l3P2,g3P0,g3P2);
6196 if ((fDebug <= 1) ||
6198 projch->Fit("fLandauGaus","0",""
6199 ,(Float_t) fChargeCoef[1]/fBeginFitCharge
6200 ,projch->GetBinCenter(projch->GetNbinsX()));
6201 chisqr = projch->GetFunction("fLandauGaus")->GetChisquare();
6205 TCanvas *cp = new TCanvas("cp","cp",50,50,600,800);
6207 projch->Fit("fLandauGaus","+",""
6208 ,(Float_t) fChargeCoef[1]/fBeginFitCharge
6209 ,projch->GetBinCenter(projch->GetNbinsX()));
6210 chisqr = projch->GetFunction("fLandauGaus")->GetChisquare();
6212 fLandauGaus->Draw("same");
6215 if ((projch->GetFunction("fLandauGaus")->GetParameter(1) > 0) &&
6216 (projch->GetFunction("fLandauGaus")->GetParError(1) <
6217 (0.05*projch->GetFunction("fLandauGaus")->GetParameter(1))) &&
6218 (chisqr < chisqrl) &&
6219 (chisqr < chisqrg)) {
6220 // Calcul of "real" coef
6221 if(fFitChargeNDB == 0){
6222 fNumberFitSuccess++;
6223 CalculChargeCoefMean(fCountDet[0],(Int_t) idect,kTRUE);
6226 CalculChargeCoefMean(fCountDet[0],(Int_t) idect,kFALSE);
6228 fChargeCoef[0] = projch->GetFunction("fLandauGaus")->GetParameter(1);
6229 chargeCoefE0 = projch->GetFunction("fLandauGaus")->GetParError(1);
6232 // Calcul of "real" coef
6233 CalculChargeCoefMean(fCountDet[0],(Int_t) idect,kFALSE);
6234 fChargeCoef[0] = -TMath::Abs(fChargeCoef[3]);
6238 AliInfo(Form("fChargeCoef[0]: %f",(Float_t) fChargeCoef[0]));
6239 AliInfo(Form("fChargeCoef[1]: %f",(Float_t) fChargeCoef[1]));
6242 if ((fDebug == 1) ||
6244 if (fChargeCoef[0] > 0.0) {
6245 fCoefCharge[0][idect]= fChargeCoef[0];
6246 fCoefChargeE[0][idect]= chargeCoefE0;
6256 //_____________________________________________________________________________
6257 void AliTRDCalibraFit::FitBisCH(TH1* projch, Int_t idect)
6260 // Fit methode for the gain factor more time consuming
6263 //Some parameters to initialise
6264 Double_t widthLandau, widthGaus, MPV, Integral;
6265 Double_t chisquarel = 0.0;
6266 Double_t chisquareg = 0.0;
6268 projch->Fit("landau","0M+",""
6269 ,(Float_t) fChargeCoef[1]/6
6270 ,projch->GetBinCenter(projch->GetNbinsX()));
6271 widthLandau = projch->GetFunction("landau")->GetParameter(2);
6272 chisquarel = projch->GetFunction("landau")->GetChisquare();
6274 projch->Fit("gaus","0M+",""
6275 ,(Float_t) fChargeCoef[1]/6
6276 ,projch->GetBinCenter(projch->GetNbinsX()));
6277 widthGaus = projch->GetFunction("gaus")->GetParameter(2);
6278 chisquareg = projch->GetFunction("gaus")->GetChisquare();
6280 MPV = (projch->GetFunction("landau")->GetParameter(1))/2;
6281 Integral = (projch->GetFunction("gaus")->Integral(0.3*fChargeCoef[1],3*fChargeCoef[1])
6282 + projch->GetFunction("landau")->Integral(0.3*fChargeCoef[1],3*fChargeCoef[1]))/2;
6284 // Setting fit range and start values
6286 //Double_t sv[4] = { l3P2, fChargeCoef[1], projch->Integral("width"), fG3P2 };
6287 //Double_t sv[4] = { fL3P2, fChargeCoef[1], fL3P0, fG3P2 };
6288 Double_t sv[4] = { widthLandau, MPV, Integral, widthGaus};
6289 Double_t pllo[4] = { 0.001, 0.001, projch->Integral()/3, 0.001};
6290 Double_t plhi[4] = { 300.0, 300.0, 30*projch->Integral(), 300.0};
6291 Double_t fp[4] = { 1.0, 1.0, 1.0, 1.0 };
6292 Double_t fpe[4] = { 1.0, 1.0, 1.0, 1.0 };
6293 fr[0] = 0.3 * fChargeCoef[1];
6294 fr[1] = 3.0 * fChargeCoef[1];
6295 fChargeCoef[2] = 0.0;
6296 Double_t chargeCoefE2 = 0.0;
6300 TF1 *fitsnr = LanGauFit(projch,&fr[0],&sv[0]
6305 Double_t projchPeak;
6306 Double_t projchFWHM;
6307 LanGauPro(fp,projchPeak,projchFWHM);
6309 if ((fp[1] > 0) && ((fpe[1] < (0.05*fp[1])) && (chisqr < chisquarel) && (chisqr < chisquareg))) {
6310 //if ((fp[1] > 0) && ((chisqr < chisquarel) && (chisqr < chisquareg))) {
6311 if(fFitChargeNDB == 2){
6312 fNumberFitSuccess++;
6313 CalculChargeCoefMean(fCountDet[0],(Int_t) idect,kTRUE);
6316 CalculChargeCoefMean(fCountDet[0],(Int_t) idect,kFALSE);
6318 fChargeCoef[2] = fp[1];
6319 chargeCoefE2 = fpe[1];
6320 //chargeCoefE2 = chisqr;
6323 CalculChargeCoefMean(fCountDet[0],(Int_t) idect,kFALSE);
6324 fChargeCoef[2] = -TMath::Abs(fChargeCoef[3]);
6328 AliInfo(Form("fChargeCoef[2]: %f",(Float_t) fChargeCoef[2]));
6329 TCanvas *cpy = new TCanvas("cpy","cpy",50,50,600,800);
6332 fitsnr->Draw("same");
6335 if ((fDebug == 1) ||
6337 if (fChargeCoef[2] > 0.0) {
6338 fCoefCharge[2][idect]= fChargeCoef[2];
6339 fCoefChargeE[2][idect]= chargeCoefE2;
6349 //_____________________________________________________________________________
6350 Double_t *AliTRDCalibraFit::CalculPolynomeLagrange2(Double_t *x, Double_t *y)
6353 // Calcul the coefficients of the polynome passant par ces trois points de degre 2
6356 Double_t *c = new Double_t[5];
6357 Double_t x0 = y[0]/((x[0]-x[1])*(x[0]-x[2]));
6358 Double_t x1 = y[1]/((x[1]-x[0])*(x[1]-x[2]));
6359 Double_t x2 = y[2]/((x[2]-x[0])*(x[2]-x[1]));
6364 c[1] = -(x0*(x[1]+x[2])+x1*(x[0]+x[2])+x2*(x[0]+x[1]));
6365 c[0] = x0*x[1]*x[2]+x1*x[0]*x[2]+x2*x[0]*x[1];
6371 //_____________________________________________________________________________
6372 Double_t *AliTRDCalibraFit::CalculPolynomeLagrange3(Double_t *x, Double_t *y)
6375 // Calcul the coefficients of the polynome passant par ces quatre points de degre 3
6378 Double_t *c = new Double_t[5];
6379 Double_t x0 = y[0]/((x[0]-x[1])*(x[0]-x[2])*(x[0]-x[3]));
6380 Double_t x1 = y[1]/((x[1]-x[0])*(x[1]-x[2])*(x[1]-x[3]));
6381 Double_t x2 = y[2]/((x[2]-x[0])*(x[2]-x[1])*(x[2]-x[3]));
6382 Double_t x3 = y[3]/((x[3]-x[0])*(x[3]-x[1])*(x[3]-x[2]));
6386 c[2] = -(x0*(x[1]+x[2]+x[3])
6387 +x1*(x[0]+x[2]+x[3])
6388 +x2*(x[0]+x[1]+x[3])
6389 +x3*(x[0]+x[1]+x[2]));
6390 c[1] = (x0*(x[1]*x[2]+x[1]*x[3]+x[2]*x[3])
6391 +x1*(x[0]*x[2]+x[0]*x[3]+x[2]*x[3])
6392 +x2*(x[0]*x[1]+x[0]*x[3]+x[1]*x[3])
6393 +x3*(x[0]*x[1]+x[0]*x[2]+x[1]*x[2]));
6395 c[0] = -(x0*x[1]*x[2]*x[3]
6398 +x3*x[0]*x[1]*x[2]);
6404 //_____________________________________________________________________________
6405 Double_t *AliTRDCalibraFit::CalculPolynomeLagrange4(Double_t *x, Double_t *y)
6408 // Calcul the coefficients of the polynome passant par ces cinqs points de degre 4
6411 Double_t *c = new Double_t[5];
6412 Double_t x0 = y[0]/((x[0]-x[1])*(x[0]-x[2])*(x[0]-x[3])*(x[0]-x[4]));
6413 Double_t x1 = y[1]/((x[1]-x[0])*(x[1]-x[2])*(x[1]-x[3])*(x[1]-x[4]));
6414 Double_t x2 = y[2]/((x[2]-x[0])*(x[2]-x[1])*(x[2]-x[3])*(x[2]-x[4]));
6415 Double_t x3 = y[3]/((x[3]-x[0])*(x[3]-x[1])*(x[3]-x[2])*(x[3]-x[4]));
6416 Double_t x4 = y[4]/((x[4]-x[0])*(x[4]-x[1])*(x[4]-x[2])*(x[4]-x[3]));
6418 c[4] = x0+x1+x2+x3+x4;
6419 c[3] = -(x0*(x[1]+x[2]+x[3]+x[4])
6420 +x1*(x[0]+x[2]+x[3]+x[4])
6421 +x2*(x[0]+x[1]+x[3]+x[4])
6422 +x3*(x[0]+x[1]+x[2]+x[4])
6423 +x4*(x[0]+x[1]+x[2]+x[3]));
6424 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])
6425 +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])
6426 +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])
6427 +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])
6428 +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]));
6430 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])
6431 +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])
6432 +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])
6433 +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])
6434 +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]));
6436 c[0] = (x0*x[1]*x[2]*x[3]*x[4]
6437 +x1*x[0]*x[2]*x[3]*x[4]
6438 +x2*x[0]*x[1]*x[3]*x[4]
6439 +x3*x[0]*x[1]*x[2]*x[4]
6440 +x4*x[0]*x[1]*x[2]*x[3]);
6446 //_____________________________________________________________________________
6447 void AliTRDCalibraFit::NormierungCharge()
6450 // Normalisation of the gain factor resulting for the fits
6453 // Calcul of the mean of choosen method by fFitChargeNDB
6455 //printf("total number of entries %d\n",fVectorFitCH->GetEntriesFast());
6456 for (Int_t k = 0; k < (Int_t) fVectorFitCH->GetEntriesFast(); k++) {
6458 Int_t detector = ((AliTRDFitCHInfo *) fVectorFitCH->At(k))->GetDetector();
6459 Float_t *coef = ((AliTRDFitCHInfo *) fVectorFitCH->At(k))->GetCoef();
6460 //printf("detector %d coef[0] %f\n",detector,coef[0]);
6461 if (GetChamber(detector) == 2) {
6464 if (GetChamber(detector) != 2) {
6467 for (Int_t j = 0; j < total; j++) {
6475 fScaleFitFactor = fScaleFitFactor / sum;
6478 fScaleFitFactor = 1.0;
6481 if ((fDebug == 3) ||
6483 if ((fFitChargeOn) &&
6484 (fCoefChargeDB[0]->GetEntries() > 0.0) &&
6485 (fCoefChargeDB[0]->GetSumOfWeights() > 0.0)) {
6486 fCoefChargeDB[0]->Scale(fCoefChargeDB[0]->GetEntries() / fCoefChargeDB[0]->GetSumOfWeights());
6488 if ((fMeanChargeOn) &&
6489 (fCoefChargeDB[1]->GetEntries() > 0.0) &&
6490 (fCoefChargeDB[1]->GetSumOfWeights() > 0.0)) {
6491 fCoefChargeDB[1]->Scale(fCoefChargeDB[1]->GetEntries() / fCoefChargeDB[1]->GetSumOfWeights());
6493 if ((fFitChargeBisOn) &&
6494 (fCoefChargeDB[2]->GetEntries() > 0.0) &&
6495 (fCoefChargeDB[2]->GetSumOfWeights() > 0.0)) {
6496 fCoefChargeDB[2]->Scale(fCoefChargeDB[2]->GetEntries() / fCoefChargeDB[2]->GetSumOfWeights());
6498 if ((fFitMeanWOn) &&
6499 (fCoefChargeDB[3]->GetEntries() > 0.0) &&
6500 (fCoefChargeDB[3]->GetSumOfWeights() > 0.0)) {
6501 fCoefChargeDB[3]->Scale(fCoefChargeDB[3]->GetEntries() / fCoefChargeDB[3]->GetSumOfWeights());
6507 //_____________________________________________________________________________
6508 TH1I *AliTRDCalibraFit::ReBin(TH1I *hist) const
6511 // Rebin of the 1D histo for the gain calibration if needed.
6512 // you have to choose fRebin, divider of fNumberBinCharge
6515 TAxis *xhist = hist->GetXaxis();
6516 TH1I *rehist = new TH1I("projrebin","",(Int_t) xhist->GetNbins()/fRebin
6517 ,xhist->GetBinLowEdge(1)
6518 ,xhist->GetBinUpEdge(xhist->GetNbins()));
6520 AliInfo(Form("fRebin: %d",fRebin));
6522 for (Int_t k = 1; k <= (Int_t) xhist->GetNbins()/fRebin; k++) {
6524 for (Int_t ji = i; ji < i+fRebin; ji++) {
6525 sum += hist->GetBinContent(ji);
6528 rehist->SetBinContent(k,sum);
6533 TCanvas *crebin = new TCanvas("crebin","",50,50,600,800);
6542 //_____________________________________________________________________________
6543 TH1F *AliTRDCalibraFit::ReBin(TH1F *hist) const
6546 // Rebin of the 1D histo for the gain calibration if needed
6547 // you have to choose fRebin divider of fNumberBinCharge
6550 TAxis *xhist = hist->GetXaxis();
6551 TH1F *rehist = new TH1F("projrebin","",(Int_t) xhist->GetNbins()/fRebin
6552 ,xhist->GetBinLowEdge(1)
6553 ,xhist->GetBinUpEdge(xhist->GetNbins()));
6555 AliInfo(Form("fRebin: %d",fRebin));
6557 for (Int_t k = 1; k <= (Int_t) xhist->GetNbins()/fRebin; k++) {
6559 for (Int_t ji = i; ji < i+fRebin; ji++) {
6560 sum += hist->GetBinContent(ji);
6563 rehist->SetBinContent(k,sum);
6568 TCanvas *crebin = new TCanvas("crebin","",50,50,600,800);
6577 //_____________________________________________________________________________
6578 TH1F *AliTRDCalibraFit::CorrectTheError(TGraphErrors *hist)
6581 // In the case of the vectors method the trees contains TGraphErrors for PH and PRF
6582 // to be able to add them after
6583 // We convert it to a TH1F to be able to applied the same fit function method
6584 // After having called this function you can not add the statistics anymore
6589 Int_t nbins = hist->GetN();
6590 Double_t *x = hist->GetX();
6591 Double_t *entries = hist->GetEX();
6592 Double_t *mean = hist->GetY();
6593 Double_t *square = hist->GetEY();
6594 fEntriesCurrent = 0;
6600 Double_t step = x[1] - x[0];
6601 Double_t minvalue = x[0] - step/2;
6602 Double_t maxvalue = x[(nbins-1)] + step/2;
6604 rehist = new TH1F("projcorrecterror","",nbins,minvalue,maxvalue);
6606 for (Int_t k = 0; k < nbins; k++) {
6607 rehist->SetBinContent(k+1,mean[k]);
6608 if (entries[k] > 0.0) {
6609 fEntriesCurrent += (Int_t) entries[k];
6610 Double_t d = TMath::Abs(square[k] - (mean[k]*mean[k]));
6611 rehist->SetBinError(k+1,TMath::Sqrt(d/entries[k]));
6614 rehist->SetBinError(k+1,0.0);
6623 //____________Some basic geometry function_____________________________________
6626 //_____________________________________________________________________________
6627 Int_t AliTRDCalibraFit::GetPlane(Int_t d) const
6630 // Reconstruct the plane number from the detector number
6633 return ((Int_t) (d % 6));
6637 //_____________________________________________________________________________
6638 Int_t AliTRDCalibraFit::GetChamber(Int_t d) const
6641 // Reconstruct the chamber number from the detector number
6645 return ((Int_t) (d % 30) / fgkNplan);
6649 //_____________________________________________________________________________
6650 Int_t AliTRDCalibraFit::GetSector(Int_t d) const
6653 // Reconstruct the sector number from the detector number
6657 return ((Int_t) (d / fg));
6662 //____________Fill and Init tree Gain, PRF, Vdrift and T0______________________
6665 //_____________________________________________________________________________
6666 void AliTRDCalibraFit::InitTreePRF()
6669 // Init the tree where the coefficients from the fit methods can be stored
6673 fPRFPad = new Float_t[2304];
6674 fPRF = new TTree("PRF","PRF");
6675 fPRF->Branch("detector",&fPRFDetector,"detector/I");
6676 fPRF->Branch("width" ,fPRFPad ,"width[2304]/F");
6678 // Set to default value for the plane 0 supposed to be the first one
6679 for (Int_t k = 0; k < 2304; k++) {
6686 //_____________________________________________________________________________
6687 void AliTRDCalibraFit::FillTreePRF(Int_t countdet)
6690 // Fill the tree with the sigma of the pad response function for the detector countdet
6693 Int_t numberofgroup = 0;
6694 fPRFDetector = countdet;
6697 if (GetChamber((Int_t)(countdet+1)) == 2) {
6698 numberofgroup = 1728;
6701 numberofgroup = 2304;
6704 // Reset to default value for the next
6705 for (Int_t k = 0; k < numberofgroup; k++) {
6706 if (GetPlane((Int_t) (countdet+1)) == 0) {
6709 if (GetPlane((Int_t) (countdet+1)) == 1) {
6712 if (GetPlane((Int_t) (countdet+1)) == 2) {
6715 if (GetPlane((Int_t) (countdet+1)) == 3) {
6718 if (GetPlane((Int_t) (countdet+1)) == 4) {
6721 if (GetPlane((Int_t) (countdet+1)) == 5) {
6730 //_____________________________________________________________________________
6731 void AliTRDCalibraFit::ConvertVectorFitCHTree()
6734 // Convert the vector stuff to a tree of 1D histos if the user
6735 // want to write it after the fill functions
6738 Int_t detector = -1;
6739 Int_t numberofgroup = 1;
6740 Float_t gainPad[2304];
6742 fGain = new TTree("Gain","Gain");
6743 fGain->Branch("detector",&detector,"detector/I");
6744 fGain->Branch("gainPad" ,gainPad ,"gainPad[2304]/F");
6746 Int_t loop = (Int_t) fVectorFitCH->GetEntriesFast();
6747 for (Int_t k = 0; k < loop; k++) {
6748 detector = ((AliTRDFitCHInfo *) fVectorFitCH->At(k))->GetDetector();
6749 if (GetChamber((Int_t) ((AliTRDFitCHInfo *) fVectorFitCH->At(k))->GetDetector()) == 2) {
6750 numberofgroup = 1728;
6753 numberofgroup = 2304;
6755 for (Int_t i = 0; i < numberofgroup; i++) {
6756 if (((AliTRDFitCHInfo *) fVectorFitCH->At(k))->GetCoef()[i] >= 0) {
6757 gainPad[i] = ((AliTRDFitCHInfo *) fVectorFitCH->At(k))->GetCoef()[i] * fScaleFitFactor;
6760 gainPad[i] = (Float_t) ((AliTRDFitCHInfo *) fVectorFitCH->At(k))->GetCoef()[i];
6764 if(numberofgroup < 2304){
6765 for(Int_t i = numberofgroup; i < 2304; i++){
6766 gainPad[i] = -100.0;
6774 //_____________________________________________________________________________
6775 void AliTRDCalibraFit::FillTreeVdrift(Int_t countdet)
6778 // Fill the tree with the drift velocities for the detector countdet
6781 Int_t numberofgroup = 0;
6782 fVdriftDetector = countdet;
6785 if (GetChamber((Int_t)(countdet+1)) == 2) {
6786 numberofgroup = 1728;
6789 numberofgroup = 2304;
6791 // Reset to default value the gain coef
6792 for (Int_t k = 0; k < numberofgroup; k++) {
6793 fVdriftPad[k] = -1.5;
6795 fVdriftDetector = -1;
6799 //_____________________________________________________________________________
6800 void AliTRDCalibraFit::InitTreePH()
6803 // Init the tree where the coefficients from the fit methods can be stored
6807 fVdriftPad = new Float_t[2304];
6808 fVdrift = new TTree("Vdrift","Vdrift");
6809 fVdrift->Branch("detector",&fVdriftDetector,"detector/I");
6810 fVdrift->Branch("vdrift" ,fVdriftPad ,"vdrift[2304]/F");
6811 // Set to default value for the plane 0 supposed to be the first one
6812 for (Int_t k = 0; k < 2304; k++) {
6813 fVdriftPad[k] = -1.5;
6815 fVdriftDetector = -1;
6819 //_____________________________________________________________________________
6820 void AliTRDCalibraFit::FillTreeT0(Int_t countdet)
6823 // Fill the tree with the t0 value for the detector countdet
6826 Int_t numberofgroup = 0;
6828 fT0Detector = countdet;
6831 if (GetChamber((Int_t) (countdet+1)) == 2) {
6832 numberofgroup = 1728;
6835 numberofgroup = 2304;
6837 // Reset to default value
6838 for (Int_t k = 0; k < numberofgroup; k++) {
6845 //_____________________________________________________________________________
6846 void AliTRDCalibraFit::InitTreeT0()
6849 // Init the tree where the coefficients from the fit methods can be stored
6853 fT0Pad = new Float_t[2304];
6854 fT0 = new TTree("T0","T0");
6855 fT0->Branch("detector",&fT0Detector,"detector/I");
6856 fT0->Branch("t0",fT0Pad,"t0[2304]/F");
6857 //Set to default value for the plane 0 supposed to be the first one
6858 for(Int_t k = 0; k < 2304; k++){
6866 //____________Private Functions________________________________________________
6869 //_____________________________________________________________________________
6870 Double_t AliTRDCalibraFit::PH(Double_t *x, Double_t *par)
6873 // Function for the fit
6876 //TF1 *fAsymmGauss = new TF1("fAsymmGauss",AsymmGauss,0,4,6);
6878 //PARAMETERS FOR FIT PH
6880 //fAsymmGauss->SetParameter(0,0.113755);
6881 //fAsymmGauss->SetParameter(1,0.350706);
6882 //fAsymmGauss->SetParameter(2,0.0604244);
6883 //fAsymmGauss->SetParameter(3,7.65596);
6884 //fAsymmGauss->SetParameter(4,1.00124);
6885 //fAsymmGauss->SetParameter(5,0.870597); // No tail cancelation
6893 Double_t dx = 0.005;
6894 Double_t xs = par[1];
6896 Double_t paras[2] = { 0.0, 0.0 };
6899 if ((xs >= par[1]) &&
6900 (xs < (par[1]+par[2]))) {
6901 //fAsymmGauss->SetParameter(0,par[0]);
6902 //fAsymmGauss->SetParameter(1,xs);
6903 //ss += fAsymmGauss->Eval(xx);
6906 ss += AsymmGauss(&xx,paras);
6908 if ((xs >= (par[1]+par[2])) &&
6909 (xs < (par[1]+par[2]+par[3]))) {
6910 //fAsymmGauss->SetParameter(0,par[0]*par[4]);
6911 //fAsymmGauss->SetParameter(1,xs);
6912 //ss += fAsymmGauss->Eval(xx);
6913 paras[0] = par[0]*par[4];
6915 ss += AsymmGauss(&xx,paras);
6924 //_____________________________________________________________________________
6925 Double_t AliTRDCalibraFit::AsymmGauss(Double_t *x, Double_t *par)
6928 // Function for the fit
6931 //par[0] = normalization
6939 Double_t par1save = par[1];
6940 //Double_t par2save = par[2];
6941 Double_t par2save = 0.0604244;
6942 //Double_t par3save = par[3];
6943 Double_t par3save = 7.65596;
6944 //Double_t par5save = par[5];
6945 Double_t par5save = 0.870597;
6946 Double_t dx = x[0] - par1save;
6948 Double_t sigma2 = par2save*par2save;
6949 Double_t sqrt2 = TMath::Sqrt(2.0);
6950 Double_t exp1 = par3save * TMath::Exp(-par3save * (dx - 0.5 * par3save * sigma2))
6951 * (1.0 - TMath::Erf((par3save * sigma2 - dx) / (sqrt2 * par2save)));
6952 Double_t exp2 = par5save * TMath::Exp(-par5save * (dx - 0.5 * par5save * sigma2))
6953 * (1.0 - TMath::Erf((par5save * sigma2 - dx) / (sqrt2 * par2save)));
6955 //return par[0]*(exp1+par[4]*exp2);
6956 return par[0] * (exp1 + 1.00124 * exp2);
6960 //_____________________________________________________________________________
6961 Double_t AliTRDCalibraFit::FuncLandauGaus(Double_t *x, Double_t *par)
6964 // Sum Landau + Gaus with identical mean
6967 Double_t valLandau = par[0] * TMath::Landau(x[0],par[1],par[2]);
6968 //Double_t valGaus = par[3] * TMath::Gaus(x[0],par[4],par[5]);
6969 Double_t valGaus = par[3] * TMath::Gaus(x[0],par[1],par[4]);
6970 Double_t val = valLandau + valGaus;
6976 //_____________________________________________________________________________
6977 Double_t AliTRDCalibraFit::LanGauFun(Double_t *x, Double_t *par)
6980 // Function for the fit
6983 // par[0]=Width (scale) parameter of Landau density
6984 // par[1]=Most Probable (MP, location) parameter of Landau density
6985 // par[2]=Total area (integral -inf to inf, normalization constant)
6986 // par[3]=Width (sigma) of convoluted Gaussian function
6988 // In the Landau distribution (represented by the CERNLIB approximation),
6989 // the maximum is located at x=-0.22278298 with the location parameter=0.
6990 // This shift is corrected within this function, so that the actual
6991 // maximum is identical to the MP parameter.
6994 // Numeric constants
6995 Double_t invsq2pi = 0.3989422804014; // (2 pi)^(-1/2)
6996 Double_t mpshift = -0.22278298; // Landau maximum location
6998 // Control constants
6999 Double_t np = 100.0; // Number of convolution steps
7000 Double_t sc = 5.0; // Convolution extends to +-sc Gaussian sigmas
7012 // MP shift correction
7013 mpc = par[1] - mpshift * par[0];
7015 // Range of convolution integral
7016 xlow = x[0] - sc * par[3];
7017 xupp = x[0] + sc * par[3];
7019 step = (xupp - xlow) / np;
7021 // Convolution integral of Landau and Gaussian by sum
7022 for (i = 1.0; i <= np/2; i++) {
7024 xx = xlow + (i-.5) * step;
7025 fland = TMath::Landau(xx,mpc,par[0]) / par[0];
7026 sum += fland * TMath::Gaus(x[0],xx,par[3]);
7028 xx = xupp - (i-.5) * step;
7029 fland = TMath::Landau(xx,mpc,par[0]) / par[0];
7030 sum += fland * TMath::Gaus(x[0],xx,par[3]);
7034 return (par[2] * step * sum * invsq2pi / par[3]);
7038 //_____________________________________________________________________________
7039 TF1 *AliTRDCalibraFit::LanGauFit(TH1 *his, Double_t *fitrange, Double_t *startvalues
7040 , Double_t *parlimitslo, Double_t *parlimitshi
7041 , Double_t *fitparams, Double_t *fiterrors
7042 , Double_t *chiSqr, Int_t *ndf)
7045 // Function for the fit
7049 Char_t funname[100];
7051 TF1 *ffitold = (TF1 *) gROOT->GetListOfFunctions()->FindObject(funname);
7056 TF1 *ffit = new TF1(funname,LanGauFun,fitrange[0],fitrange[1],4);
7057 ffit->SetParameters(startvalues);
7058 ffit->SetParNames("Width","MP","Area","GSigma");
7060 for (i = 0; i < 4; i++) {
7061 ffit->SetParLimits(i,parlimitslo[i],parlimitshi[i]);
7064 his->Fit(funname,"RB0"); // Fit within specified range, use ParLimits, do not plot
7066 ffit->GetParameters(fitparams); // Obtain fit parameters
7067 for (i = 0; i < 4; i++) {
7068 fiterrors[i] = ffit->GetParError(i); // Obtain fit parameter errors
7070 chiSqr[0] = ffit->GetChisquare(); // Obtain chi^2
7071 ndf[0] = ffit->GetNDF(); // Obtain ndf
7073 return (ffit); // Return fit function
7077 //_____________________________________________________________________________
7078 Int_t AliTRDCalibraFit::LanGauPro(Double_t *params, Double_t &maxx, Double_t &fwhm)
7081 // Function for the fit
7094 Int_t maxcalls = 10000;
7096 // Search for maximum
7097 p = params[1] - 0.1 * params[0];
7098 step = 0.05 * params[0];
7102 while ((l != lold) && (i < maxcalls)) {
7106 l = LanGauFun(&x,params);
7108 step = -step / 10.0;
7113 if (i == maxcalls) {
7119 // Search for right x location of fy
7120 p = maxx + params[0];
7126 while ( (l != lold) && (i < maxcalls) ) {
7131 l = TMath::Abs(LanGauFun(&x,params) - fy);
7144 // Search for left x location of fy
7146 p = maxx - 0.5 * params[0];
7152 while ((l != lold) && (i < maxcalls)) {
7156 l = TMath::Abs(LanGauFun(&x,params) - fy);
7158 step = -step / 10.0;
7163 if (i == maxcalls) {
7174 //_____________________________________________________________________________
7175 Double_t AliTRDCalibraFit::GausConstant(Double_t *x, Double_t *par)
7178 // Gaus with identical mean
7181 Double_t gauss = par[0] * TMath::Gaus(x[0],0.0,par[1])+par[2];