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 "./Cal/AliTRDCalROC.h"
79 #include "./Cal/AliTRDCalPad.h"
80 #include "./Cal/AliTRDCalDet.h"
83 ClassImp(AliTRDCalibraFit)
85 AliTRDCalibraFit* AliTRDCalibraFit::fgInstance = 0;
86 Bool_t AliTRDCalibraFit::fgTerminated = kFALSE;
88 //_____________singleton implementation_________________________________________________
89 AliTRDCalibraFit *AliTRDCalibraFit::Instance()
92 // Singleton implementation
95 if (fgTerminated != kFALSE) {
99 if (fgInstance == 0) {
100 fgInstance = new AliTRDCalibraFit();
107 //______________________________________________________________________________________
108 void AliTRDCalibraFit::Terminate()
111 // Singleton implementation
112 // Deletes the instance of this class
115 fgTerminated = kTRUE;
117 if (fgInstance != 0) {
124 //______________________________________________________________________________________
125 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 fGeo = new AliTRDgeometry();
177 fCalibraMode = new AliTRDCalibraMode();
180 for (Int_t i = 0; i < 3; i++) {
181 fWriteCoef[i] = kFALSE;
185 for (Int_t k = 0; k < 3; k++) {
189 for (Int_t i = 0; i < 3; i++) {
198 //______________________________________________________________________________________
199 AliTRDCalibraFit::AliTRDCalibraFit(const AliTRDCalibraFit &c)
205 ,fFitLagrPolOn(kFALSE)
206 ,fTakeTheMaxPH(kFALSE)
209 ,fBeginFitCharge(0.0)
215 ,fMeanChargeOn(kFALSE)
216 ,fFitChargeBisOn(kFALSE)
217 ,fFitChargeOn(kFALSE)
224 ,fNumberFitSuccess(0)
241 ,fScaleFitFactor(0.0)
252 //____________________________________________________________________________________
253 AliTRDCalibraFit::~AliTRDCalibraFit()
256 // AliTRDCalibraFit destructor
267 //_____________________________________________________________________________
268 void AliTRDCalibraFit::Destroy()
281 //_____________________________________________________________________________
282 void AliTRDCalibraFit::ClearTree()
307 //_____________________________________________________________________________
308 void AliTRDCalibraFit::Init()
311 // Init some default values
315 fWriteNameCoef = "TRD.coefficient.root";
319 fBeginFitCharge = 3.5;
324 // Internal variables
326 // Variables in the loop
327 for (Int_t k = 0; k < 4; k++) {
328 fChargeCoef[k] = 1.0;
329 fVdriftCoef[k] = 1.5;
332 fChargeCoef[4] = 1.0;
333 for (Int_t i = 0; i < 3; i++) {
337 // Local database to be changed
342 //____________Functions fit Online CH2d________________________________________
343 Bool_t AliTRDCalibraFit::FitCHOnline(TH2I *ch)
346 // Fit the 1D histos, projections of the 2D ch on the Xaxis, for each
347 // calibration group normalized the resulted coefficients (to 1 normally)
348 // and write the results in a tree
352 if((fFitChargeNDB == 0) && (!fFitChargeOn)){
353 AliInfo("You have choosen to write the default fit method but it is not on!");
356 if((fFitChargeNDB == 1) && (!fMeanChargeOn)){
357 AliInfo("You have choosen to write the mean method but it is not on!");
360 if((fFitChargeNDB == 2) && (!fFitChargeBisOn)){
361 AliInfo("You have choosen to write the second fit method but it is not on!");
364 if((fFitChargeNDB == 4) && (!fFitMeanWOn)){
365 AliInfo("You have choosen to write the mean w method but it is not on!");
370 // Number of Xbins (detectors or groups of pads)
371 TAxis *xch = ch->GetXaxis();
372 Int_t nbins = xch->GetNbins();
373 TAxis *yph = ch->GetYaxis();
374 Int_t nybins = yph->GetNbins();
375 if (!InitFit(nbins,0)) {
378 fStatisticMean = 0.0;
380 fNumberFitSuccess = 0;
383 // Init fCountDet and fCount
384 InitfCountDetAndfCount(0);
386 // Beginning of the loop betwwen dect1 and dect2
387 for (Int_t idect = fDect1[0]; idect < fDect2[0]; idect++) {
389 TH1I *projch = (TH1I *) ch->ProjectionY("projch",idect+1,idect+1,(Option_t *)"e");
390 projch->SetDirectory(0);
392 // Determination of fNnZ, fNnRphi, fNfragZ and fNfragRphi
393 UpdatefCountDetAndfCount(idect,0);
395 // Reconstruction of the row and pad group: rowmin, row max ...
396 ReconstructFitRowMinRowMax(idect, 0);
398 // Number of entries for this calibration group
399 Double_t nentries = 0.0;
401 for (Int_t k = 0; k < nybins; k++) {
402 nentries += ch->GetBinContent(ch->GetBin(idect+1,k+1));
403 mean += projch->GetBinCenter(k+1)*projch->GetBinContent(k+1);
411 // Rebin and statistic stuff
414 projch = ReBin((TH1I *) projch);
416 // This detector has not enough statistics or was off
417 if (nentries < fMinEntries) {
418 // Fill with the default infos
419 NotEnoughStatistic(idect,0);
427 // Statistics of the group fitted
428 AliInfo(Form("For the group number %d there are %f stats",idect,nentries));
429 fStatisticMean += nentries;
433 // Method Mean and fit
434 // idect is egal for fDebug = 0 and 2, only to fill the hist
435 fChargeCoef[1] = mean;
437 FitMean((TH1 *) projch,(Int_t) (idect-fDect1[0]),nentries);
440 FitCH((TH1 *) projch,(Int_t) (idect-fDect1[0]));
442 if(fFitChargeBisOn) {
443 FitBisCH((TH1 *) projch,(Int_t) (idect-fDect1[0]));
446 FitMeanW((TH1 *) projch,(Int_t) (idect-fDect1[0]));
449 // Visualise the detector for fDebug 3 or 4
450 // Here is the reconstruction of the pad and row group is used!
455 FillInfosFit(idect,0);
472 // 0 no plot, 1 and 4 error plot, 3 and 4 DB plot
483 if (fNumberFit > 0) {
484 AliInfo(Form("There are %d with at least one entries.",fNumberEnt));
485 AliInfo(Form("%d fits have been proceeded (sucessfully or not...).",fNumberFit));
486 AliInfo(Form("There is a mean statistic of: %d over these fitted histograms and %d successfulled fits"
487 , (Int_t) fStatisticMean/fNumberFit, fNumberFitSuccess));
488 fStatisticMean = fStatisticMean / fNumberFit;
491 AliInfo(Form("There are %d with at least one entries. There is no fit!",fNumberEnt));
495 ConvertVectorFitCHTree();
504 //____________Functions fit Online CH2d________________________________________
505 Bool_t AliTRDCalibraFit::FitCHOnline()
508 // Reconstruct a 1D histo from the vectorCH for each calibration group,
509 // fit the histo, normalized the resulted coefficients (to 1 normally)
510 // and write the results in a tree
514 if((fFitChargeNDB == 0) && (!fFitChargeOn)){
515 AliInfo("You have choosen to write the default fit method but it is not on!");
518 if((fFitChargeNDB == 1) && (!fMeanChargeOn)){
519 AliInfo("You have choosen to write the mean method but it is not on!");
522 if((fFitChargeNDB == 2) && (!fFitChargeBisOn)){
523 AliInfo("You have choosen to write the second fit method but it is not on!");
526 if((fFitChargeNDB == 4) && (!fFitMeanWOn)){
527 AliInfo("You have choosen to write the mean w method but it is not on!");
533 if (!fCalibraVector) {
534 AliError("You have first to set the calibravector before using this function!");
538 // Number of Xbins (detectors or groups of pads)
542 fStatisticMean = 0.0;
544 fNumberFitSuccess = 0;
547 // Init fCountDet and fCount
548 InitfCountDetAndfCount(0);
550 // Beginning of the loop between dect1 and dect2
551 for (Int_t idect = fDect1[0]; idect < fDect2[0]; idect++) {
553 // Search if the group is in the VectorCH
554 Int_t place = fCalibraVector->SearchInVector(idect,0);
561 projch = fCalibraVector->ConvertVectorCTHisto(place,(const char *) name);
562 projch->SetDirectory(0);
565 // Determination of fNnZ, fNnRphi, fNfragZ and fNfragRphi
566 UpdatefCountDetAndfCount(idect,0);
568 // Reconstruction of the row and pad group: rowmin, row max ...
569 ReconstructFitRowMinRowMax(idect,0);
571 // Number of entries and mean
572 Double_t nentries = 0.0;
575 for (Int_t k = 0; k < fCalibraVector->GetNumberBinCharge(); k++) {
576 nentries += projch->GetBinContent(k+1);
577 mean += projch->GetBinCenter(k+1)*projch->GetBinContent(k+1);
585 // Rebin and statistic stuff
589 projch = ReBin((TH1F *) projch);
592 // This detector has not enough statistics or was not found in VectorCH
595 (nentries < fMinEntries))) {
597 // Fill with the default infos
598 NotEnoughStatistic(idect,0);
609 // Statistic of the histos fitted
610 AliInfo(Form("For the group number %d there are %f stats",idect,nentries));
611 fStatisticMean += nentries;
615 // Method Mean and fit
616 // idect is egal for fDebug = 0 and 2, only to fill the hist
617 fChargeCoef[1] = mean;
619 FitMean((TH1 *) projch,(Int_t) (idect-fDect1[0]),nentries);
622 FitCH((TH1 *) projch,(Int_t) (idect-fDect1[0]));
624 if(fFitChargeBisOn) {
625 FitBisCH((TH1 *) projch,(Int_t) (idect-fDect1[0]));
628 FitMeanW((TH1 *) projch,(Int_t) (idect-fDect1[0]));
631 // Visualise the detector for fDebug 3 or 4
632 // Here is the reconstruction of the pad and row group is used!
638 FillInfosFit(idect,0);
655 // 0 no plot, 1 and 4 error plot, 3 and 4 DB plot
666 if (fNumberFit > 0) {
667 AliInfo(Form("There are %d with at least one entries.",fNumberEnt));
668 AliInfo(Form("%d fits have been proceeded (sucessfully or not...).",fNumberFit));
669 AliInfo(Form("There is a mean statistic of: %d over these fitted histograms and %d successfulled fits"
670 ,(Int_t) fStatisticMean/fNumberFit, fNumberFitSuccess));
671 fStatisticMean = fStatisticMean / fNumberFit;
674 AliInfo(Form("There are %d with at least one entries. There is no fit!",fNumberEnt));
678 ConvertVectorFitCHTree();
687 //____________Functions fit Online CH2d________________________________________
688 Bool_t AliTRDCalibraFit::FitCHOnline(TTree *tree)
691 // Look if the calibration group can be found in the tree, if yes take the
692 // histo, fit it, normalized the resulted coefficients (to 1 normally) and
693 // write the results in a tree
697 if((fFitChargeNDB == 0) && (!fFitChargeOn)){
698 AliInfo("You have choosen to write the default fit method but it is not on!");
701 if((fFitChargeNDB == 1) && (!fMeanChargeOn)){
702 AliInfo("You have choosen to write the mean method but it is not on!");
705 if((fFitChargeNDB == 2) && (!fFitChargeBisOn)){
706 AliInfo("You have choosen to write the second fit method but it is not on!");
709 if((fFitChargeNDB == 4) && (!fFitMeanWOn)){
710 AliInfo("You have choosen to write the mean w method but it is not on!");
715 // Number of Xbins (detectors or groups of pads)
719 fStatisticMean = 0.0;
721 fNumberFitSuccess = 0;
725 fCalibraVector = new AliTRDCalibraVector();
728 // Init fCountDet and fCount
729 InitfCountDetAndfCount(0);
731 tree->SetBranchAddress("histo",&projch);
732 TObjArray *vectorplace = fCalibraVector->ConvertTreeVector(tree);
734 // Beginning of the loop between dect1 and dect2
735 for (Int_t idect = fDect1[0]; idect < fDect2[0]; idect++) {
737 //Search if the group is in the VectorCH
738 Int_t place = fCalibraVector->SearchInTreeVector(vectorplace,idect);
743 tree->GetEntry(place);
746 // Determination of fNnZ, fNnRphi, fNfragZ and fNfragRphi
747 UpdatefCountDetAndfCount(idect,0);
749 // Reconstruction of the row and pad group: rowmin, row max ...
750 ReconstructFitRowMinRowMax(idect,0);
752 // Number of entries and mean
753 Double_t nentries = 0.0;
756 for (Int_t k = 0; k < projch->GetXaxis()->GetNbins(); k++) {
757 nentries += projch->GetBinContent(k+1);
758 mean += projch->GetBinCenter(k+1)*projch->GetBinContent(k+1);
767 // Rebin and statistic stuff
771 projch = ReBin((TH1F *) projch);
774 // This detector has not enough statistics or was not found in VectorCH
777 (nentries < fMinEntries))) {
779 // Fill with the default infos
780 NotEnoughStatistic(idect,0);
786 // Statistics of the group fitted
787 AliInfo(Form("For the group number %d there are %f stats",idect,nentries));
789 fStatisticMean += nentries;
791 // Method Mean and fit
792 // idect is egal for fDebug = 0 and 2, only to fill the hist
793 fChargeCoef[1] = mean;
795 FitMean((TH1 *) projch,(Int_t) (idect-fDect1[0]),nentries);
798 FitCH((TH1 *) projch,(Int_t) (idect-fDect1[0]));
800 if(fFitChargeBisOn) {
801 FitBisCH((TH1 *) projch,(Int_t) (idect-fDect1[0]));
804 FitMeanW((TH1 *) projch,(Int_t) (idect-fDect1[0]));
807 // Visualise the detector for fDebug 3 or 4
808 // Here is the reconstruction of the pad and row group is used!
814 FillInfosFit(idect,0);
826 // 0 no plot, 1 and 4 error plot, 3 and 4 DB plot
837 if (fNumberFit > 0) {
838 AliInfo(Form("There are %d with at least one entries.",fNumberEnt));
839 AliInfo(Form("%d fits have been proceeded (sucessfully or not...).",fNumberFit));
840 AliInfo(Form("There is a mean statistic of: %d over these fitted histograms and %d successfulled fits"
841 ,(Int_t) fStatisticMean/fNumberFit,fNumberFitSuccess));
842 fStatisticMean = fStatisticMean / fNumberFit;
845 AliInfo(Form("There are %d with at least one entries. There is no fit!",fNumberEnt));
849 ConvertVectorFitCHTree();
859 //________________functions fit Online PH2d____________________________________
860 Bool_t AliTRDCalibraFit::FitPHOnline(TProfile2D *ph)
863 // Take the 1D profiles (average pulse height), projections of the 2D PH
864 // on the Xaxis, for each calibration group
865 // Fit or use the slope of the average pulse height to reconstruct the
866 // drift velocity write the results in a tree
867 // A first calibration of T0 is also made using the same method (slope method)
871 if((fFitPHNDB == 0) && (!fFitPHOn)){
872 AliInfo("You have choosen to write the fit method but it is not on!");
875 if((fFitPHNDB == 1) && (!fFitPol2On)){
876 AliInfo("You have choosen to write the Pol2 method but it is not on!");
879 if((fFitPHNDB == 3) && (!fFitLagrPolOn)){
880 AliInfo("You have choosen to write the LagrPol2 method but it is not on!");
884 // Number of Xbins (detectors or groups of pads)
885 TAxis *xph = ph->GetXaxis();
886 TAxis *yph = ph->GetYaxis();
887 Int_t nbins = xph->GetNbins();
888 Int_t nybins = yph->GetNbins();
889 if (!InitFit(nbins,1)) {
892 fStatisticMean = 0.0;
894 fNumberFitSuccess = 0;
897 // Init fCountDet and fCount
898 InitfCountDetAndfCount(1);
900 // Beginning of the loop
901 for (Int_t idect = fDect1[1]; idect < fDect2[1]; idect++) {
903 TH1D *projph = (TH1D *) ph->ProjectionY("projph",idect+1,idect+1,(Option_t *) "e");
904 projph->SetDirectory(0);
906 // Number of entries for this calibration group
907 Double_t nentries = 0;
908 for (Int_t k = 0; k < nybins; k++) {
909 nentries += ph->GetBinEntries(ph->GetBin(idect+1,k+1));
915 // Determination of fNnZ, fNnRphi, fNfragZ and fNfragRphi
916 UpdatefCountDetAndfCount(idect,1);
918 // Reconstruction of the row and pad group: rowmin, row max ...
919 ReconstructFitRowMinRowMax(idect,1);
921 // Rebin and statistic stuff
922 // This detector has not enough statistics or was off
923 if (nentries < fMinEntries) {
925 // Fill with the default values
926 NotEnoughStatistic(idect,1);
937 // Statistics of the histos fitted
938 AliInfo(Form("For the group number %d there are %f stats",idect,nentries));
940 fStatisticMean += nentries;
942 // Calcul of "real" coef
943 CalculVdriftCoefMean(fCountDet[1],(Int_t) (idect - fDect1[1]));
944 CalculT0CoefMean(fCountDet[1],(Int_t) (idect - fDect1[1]));
946 // Method Mean and fit
947 // idect is egal for fDebug = 0 and 2, only to fill the hist
949 FitPente((TH1 *) projph,(Int_t) (idect - fDect1[1]));
952 FitLagrangePoly((TH1 *) projph,(Int_t) (idect - fDect1[1]));
955 FitPH((TH1 *) projph,(Int_t) (idect - fDect1[1]));
958 // Visualise the detector for fDebug 3 or 4
959 // Here is the reconstruction of the pad and row group is used!
964 // Fill the tree if end of a detector or only the pointer to the branch!!!
965 FillInfosFit(idect,1);
975 // 0 no plot, 1 and 4 error plot, 3 and 4 DB plot
988 if (fNumberFit > 0) {
989 AliInfo(Form("There are %d with at least one entries.",fNumberEnt));
990 AliInfo(Form("%d fits have been proceeded (sucessfully or not...).",fNumberFit));
991 AliInfo(Form("There is a mean statistic of: %d over these fitted histograms and %d successfulled fits"
992 ,(Int_t) fStatisticMean/fNumberFit,fNumberFitSuccess));
993 fStatisticMean = fStatisticMean / fNumberFit;
996 AliInfo(Form("There are %d with at least one entries. There is no fit!",fNumberEnt));
1008 //____________Functions fit Online PH2d________________________________________
1009 Bool_t AliTRDCalibraFit::FitPHOnline()
1012 // Reconstruct the average pulse height from the vectorPH for each
1013 // calibration group
1014 // Fit or use the slope of the average pulse height to reconstruct the
1015 // drift velocity write the results in a tree
1016 // A first calibration of T0 is also made using the same method (slope method)
1020 if((fFitPHNDB == 0) && (!fFitPHOn)){
1021 AliInfo("You have choosen to write the fit method but it is not on!");
1024 if((fFitPHNDB == 1) && (!fFitPol2On)){
1025 AliInfo("You have choosen to write the Pol2 method but it is not on!");
1028 if((fFitPHNDB == 3) && (!fFitLagrPolOn)){
1029 AliInfo("You have choosen to write the LagrPol2 method but it is not on!");
1034 if (!fCalibraVector) {
1035 AliError("You have first to set the calibravector before using this function!");
1040 // Number of Xbins (detectors or groups of pads)
1041 if (!InitFit(0,1)) {
1044 fStatisticMean = 0.0;
1046 fNumberFitSuccess = 0;
1049 // Init fCountDet and fCount
1050 InitfCountDetAndfCount(1);
1052 // Beginning of the loop
1053 for (Int_t idect = fDect1[1]; idect < fDect2[1]; idect++) {
1055 // Search if the group is in the VectorCH
1056 Int_t place = fCalibraVector->SearchInVector(idect,1);
1065 projph = CorrectTheError((TGraphErrors *) (fCalibraVector->ConvertVectorPHisto(place,(const char *) name)));
1066 projph->SetDirectory(0);
1069 // Determination of fNnZ, fNnRphi, fNfragZ and fNfragRphi
1070 UpdatefCountDetAndfCount(idect,1);
1072 // Reconstruction of the row and pad group: rowmin, row max ...
1073 ReconstructFitRowMinRowMax(idect,1);
1075 // Rebin and statistic stuff
1076 // This detector has not enough statistics or was off
1077 if ((place == -1) ||
1079 (fEntriesCurrent < fMinEntries))) {
1081 // Fill with the default values
1082 NotEnoughStatistic(idect,1);
1093 // Statistic of the histos fitted
1094 AliInfo(Form("For the group number %d there are %d stats",idect,fEntriesCurrent));
1096 fStatisticMean += fEntriesCurrent;
1098 // Calcul of "real" coef
1099 CalculVdriftCoefMean(fCountDet[1],(Int_t) (idect - fDect1[1]));
1100 CalculT0CoefMean(fCountDet[1],(Int_t) (idect - fDect1[1]));
1102 // Method Mean and fit
1103 // idect is egal for fDebug = 0 and 2, only to fill the hist
1105 FitPente((TH1 *) projph,(Int_t) (idect - fDect1[1]));
1108 FitLagrangePoly((TH1 *) projph,(Int_t) (idect - fDect1[1]));
1111 FitPH((TH1 *) projph,(Int_t) (idect - fDect1[1]));
1114 // Visualise the detector for fDebug 3 or 4
1115 // Here is the reconstruction of the pad and row group is used!
1121 // Fill the tree if end of a detector or only the pointer to the branch!!!
1122 FillInfosFit(idect,1);
1133 // 0 no plot, 1 and 4 error plot, 3 and 4 DB plot
1134 if ((fDebug == 1) ||
1139 if ((fDebug == 4) ||
1146 if (fNumberFit > 0) {
1147 AliInfo(Form("There are %d with at least one entries.",fNumberEnt));
1148 AliInfo(Form("%d fits have been proceeded (sucessfully or not...).",fNumberFit));
1149 AliInfo(Form("There is a mean statistic of: %d over these fitted histograms and %d successfulled fits"
1150 ,(Int_t) fStatisticMean/fNumberFit,fNumberFitSuccess));
1151 fStatisticMean = fStatisticMean / fNumberFit;
1154 AliInfo(Form("There are %d with at least one entries. There is no fit!",fNumberEnt));
1157 // Write the things!
1158 if (fWriteCoef[1]) {
1166 //____________Functions fit Online PH2d________________________________________
1167 Bool_t AliTRDCalibraFit::FitPHOnline(TTree *tree)
1170 // Look if the calibration group can be found in the tree, if yes take the
1171 // histo, fit it, and write the results in a tree
1172 // A first calibration of T0 is also made using the same method (slope method)
1176 if ((fFitPHNDB == 0) && (!fFitPHOn)){
1177 AliInfo("You have choosen to write the fit method but it is not on!");
1180 if ((fFitPHNDB == 1) && (!fFitPol2On)){
1181 AliInfo("You have choosen to write the Pol2 method but it is not on!");
1184 if ((fFitPHNDB == 3) && (!fFitLagrPolOn)){
1185 AliInfo("You have choosen to write the LagrPol2 method but it is not on!");
1189 // Number of Xbins (detectors or groups of pads)
1190 if (!InitFit(0,1)) {
1193 fStatisticMean = 0.0;
1195 fNumberFitSuccess = 0;
1199 fCalibraVector = new AliTRDCalibraVector();
1201 // Init fCountDet and fCount
1202 InitfCountDetAndfCount(1);
1203 TGraphErrors *projphtree = 0x0;
1204 tree->SetBranchAddress("histo",&projphtree);
1205 TObjArray *vectorplace = fCalibraVector->ConvertTreeVector(tree);
1207 // Beginning of the loop
1208 for (Int_t idect = fDect1[1]; idect < fDect2[1]; idect++) {
1210 // Search if the group is in the VectorCH
1211 Int_t place = fCalibraVector->SearchInTreeVector(vectorplace,idect);
1219 tree->GetEntry(place);
1220 projph = CorrectTheError(projphtree);
1223 // Determination of fNnZ, fNnRphi, fNfragZ and fNfragRphi
1224 UpdatefCountDetAndfCount(idect,1);
1226 // Reconstruction of the row and pad group: rowmin, row max ...
1227 ReconstructFitRowMinRowMax(idect,1);
1229 // Rebin and statistic stuff
1230 // This detector has not enough statistics or was off
1233 (fEntriesCurrent < fMinEntries))) {
1235 // Fill with the default values
1236 NotEnoughStatistic(idect,1);
1247 // Statistics of the histos fitted
1248 AliInfo(Form("For the group number %d there are %d stats",idect,fEntriesCurrent));
1250 fStatisticMean += fEntriesCurrent;
1252 // Calcul of "real" coef
1253 CalculVdriftCoefMean(fCountDet[1],(Int_t) (idect - fDect1[1]));
1254 CalculT0CoefMean(fCountDet[1],(Int_t) (idect - fDect1[1]));
1256 // Method Mean and fit
1257 // idect is egal for fDebug = 0 and 2, only to fill the hist
1259 FitPente((TH1 *) projph,(Int_t) (idect - fDect1[1]));
1262 FitLagrangePoly((TH1 *) projph,(Int_t) (idect - fDect1[1]));
1265 FitPH((TH1 *) projph,(Int_t) (idect - fDect1[1]));
1268 // Visualise the detector for fDebug 3 or 4
1269 // Here is the reconstruction of the pad and row group is used!
1275 // Fill the tree if end of a detector or only the pointer to the branch!!!
1276 FillInfosFit(idect,1);
1286 // 0 no plot, 1 and 4 error plot, 3 and 4 DB plot
1287 if ((fDebug == 1) ||
1292 if ((fDebug == 4) ||
1299 if (fNumberFit > 0) {
1300 AliInfo(Form("There are %d with at least one entries.",fNumberEnt));
1301 AliInfo(Form("%d fits have been proceeded (sucessfully or not...).",fNumberFit));
1302 AliInfo(Form("There is a mean statistic of: %d over these fitted histograms and %d successfulled fits"
1303 ,(Int_t) fStatisticMean/fNumberFit,fNumberFitSuccess));
1304 fStatisticMean = fStatisticMean / fNumberFit;
1307 AliInfo(Form("There are %d with at least one entries. There is no fit!",fNumberEnt));
1310 // Write the things!
1311 if (fWriteCoef[1]) {
1319 //____________Functions fit Online PRF2d_______________________________________
1320 Bool_t AliTRDCalibraFit::FitPRFOnline(TProfile2D *prf)
1323 // Take the 1D profiles (pad response function), projections of the 2D PRF
1324 // on the Xaxis, for each calibration group
1325 // Fit with a gaussian to reconstruct the sigma of the pad response function
1326 // write the results in a tree
1330 if ((fFitPRFNDB == 2) && (!fRMSPRFOn)){
1331 AliInfo("You have choosen to write the RMS method but it is not on!");
1334 if ((fFitPRFNDB == 0) && (!fFitPRFOn)){
1335 AliInfo("You have choosen to write the fit method but it is not on!");
1339 // Number of Xbins (detectors or groups of pads)
1340 TAxis *xprf = prf->GetXaxis();
1341 TAxis *yprf = prf->GetYaxis();
1342 Int_t nybins = yprf->GetNbins();
1343 Int_t nbins = xprf->GetNbins();
1344 if (!InitFit(nbins,2)) {
1347 fStatisticMean = 0.0;
1349 fNumberFitSuccess = 0;
1352 // Init fCountDet and fCount
1353 InitfCountDetAndfCount(2);
1355 // Beginning of the loop
1356 for (Int_t idect = fDect1[2]; idect < fDect2[2]; idect++) {
1358 TH1D *projprf = (TH1D *) prf->ProjectionY("projprf",idect+1,idect+1,(Option_t *) "e");
1359 projprf->SetDirectory(0);
1361 // Number of entries for this calibration group
1362 Double_t nentries = 0;
1363 for (Int_t k = 0; k < nybins; k++) {
1364 nentries += prf->GetBinEntries(prf->GetBin(idect+1,k+1));
1366 if(nentries > 0) fNumberEnt++;
1368 // Determination of fNnZ, fNnRphi, fNfragZ and fNfragRphi
1369 UpdatefCountDetAndfCount(idect,2);
1371 // Reconstruction of the row and pad group: rowmin, row max ...
1372 ReconstructFitRowMinRowMax(idect,2);
1374 // Rebin and statistic stuff
1375 // This detector has not enough statistics or was off
1376 if (nentries < fMinEntries) {
1378 // Fill with the default values
1379 NotEnoughStatistic(idect,2);
1390 // Statistics of the histos fitted
1391 AliInfo(Form("For the group number %d there are %f stats",idect,nentries));
1393 fStatisticMean += nentries;
1395 // Calcul of "real" coef
1396 if ((fDebug == 1) ||
1398 CalculPRFCoefMean(fCountDet[2],(Int_t) (idect - fDect1[2]));
1401 // Method Mean and fit
1402 // idect is egal for fDebug = 0 and 2, only to fill the hist
1404 FitPRF((TH1 *) projprf,(Int_t) (idect - fDect1[2]));
1407 RmsPRF((TH1 *) projprf,(Int_t) (idect - fDect1[2]));
1411 // Visualise the detector for fDebug 3 or 4
1412 // Here is the reconstruction of the pad and row group is used!
1417 // Fill the tree if end of a detector or only the pointer to the branch!!!
1418 FillInfosFit(idect,2);
1428 // No plot, 1 and 4 error plot, 3 and 4 DB plot
1429 if ((fDebug == 1) ||
1433 if ((fDebug == 4) ||
1439 if (fNumberFit > 0) {
1440 AliInfo(Form("There are %d with at least one entries.",fNumberEnt));
1441 AliInfo(Form("%d fits have been proceeded (sucessfully or not...).",fNumberFit));
1442 AliInfo(Form("There is a mean statistic of: %d over these fitted histograms and %d successfulled fits"
1443 ,(Int_t) fStatisticMean/fNumberFit,fNumberFitSuccess));
1444 fStatisticMean = fStatisticMean / fNumberFit;
1447 AliInfo(Form("There are %d with at least one entries. There is no fit!",fNumberEnt));
1450 // Write the things!
1451 if (fWriteCoef[2]) {
1459 //____________Functions fit Online PRF2d_______________________________________
1460 Bool_t AliTRDCalibraFit::FitPRFOnline()
1463 // Reconstruct the 1D histo (pad response function) from the vectorPRD for
1464 // each calibration group
1465 // Fit with a gaussian to reconstruct the sigma of the pad response function
1466 // write the results in a tree
1470 if ((fFitPRFNDB == 2) && (!fRMSPRFOn)){
1471 AliInfo("You have choosen to write the RMS method but it is not on!");
1474 if ((fFitPRFNDB == 0) && (!fFitPRFOn)){
1475 AliInfo("You have choosen to write the fit method but it is not on!");
1480 if (!fCalibraVector) {
1481 AliError("You have first to set the calibravector before using this function!");
1485 // Number of Xbins (detectors or groups of pads)
1486 if (!InitFit(0,2)) {
1489 fStatisticMean = 0.0;
1491 fNumberFitSuccess = 0;
1494 // Init fCountDet and fCount
1495 InitfCountDetAndfCount(2);
1497 // Beginning of the loop
1498 for (Int_t idect = fDect1[2]; idect < fDect2[2]; idect++) {
1500 // Search if the group is in the VectorCH
1501 Int_t place = fCalibraVector->SearchInVector(idect,2);
1504 TH1F *projprf = 0x0;
1505 TString name("PRF");
1510 projprf = CorrectTheError((TGraphErrors *) (fCalibraVector->ConvertVectorPHisto(place,(const char *)name)));
1511 projprf->SetDirectory(0);
1514 // Determination of fNnZ, fNnRphi, fNfragZ and fNfragRphi
1515 UpdatefCountDetAndfCount(idect,2);
1517 // Reconstruction of the row and pad group: rowmin, row max ...
1518 ReconstructFitRowMinRowMax(idect,2);
1520 // Rebin and statistic stuff
1521 // This detector has not enough statistics or was off
1522 if ((place == -1) ||
1524 (fEntriesCurrent < fMinEntries))) {
1526 // Fill with the default values
1527 NotEnoughStatistic(idect,2);
1538 // Statistic of the histos fitted
1539 AliInfo(Form("For the group number %d there are %d stats", idect,fEntriesCurrent));
1541 fStatisticMean += fEntriesCurrent;
1543 // Calcul of "real" coef
1544 if ((fDebug == 1) ||
1546 CalculPRFCoefMean(fCountDet[2],(Int_t) (idect-fDect1[2]));
1549 // Method Mean and fit
1550 // idect is egal for fDebug = 0 and 2, only to fill the hist
1552 FitPRF((TH1 *) projprf,(Int_t) (idect-fDect1[2]));
1555 RmsPRF((TH1 *) projprf,(Int_t) (idect-fDect1[2]));
1558 // Visualise the detector for fDebug 3 or 4
1559 // Here is the reconstruction of the pad and row group is used!
1563 // Fill the tree if end of a detector or only the pointer to the branch!!!
1564 FillInfosFit(idect,2);
1574 // No plot, 1 and 4 error plot, 3 and 4 DB plot
1575 if ((fDebug == 1) ||
1579 if ((fDebug == 4) ||
1585 if (fNumberFit > 0) {
1586 AliInfo(Form("There are %d with at least one entries.",fNumberEnt));
1587 AliInfo(Form("%d fits have been proceeded (sucessfully or not...).",fNumberFit));
1588 AliInfo(Form("There is a mean statistic of: %d over these fitted histograms and %d successfulled fits"
1589 ,(Int_t) fStatisticMean/fNumberFit,fNumberFitSuccess));
1592 AliInfo(Form("There are %d with at least one entries. There is no fit!",fNumberEnt));
1595 // Write the things!
1596 if (fWriteCoef[2]) {
1604 //____________Functions fit Online PRF2d_______________________________________
1605 Bool_t AliTRDCalibraFit::FitPRFOnline(TTree *tree)
1608 // Look if the calibration group can be found in the tree, if yes take
1609 // the histo, fit it, and write the results in a tree
1613 if ((fFitPRFNDB == 2) && (!fRMSPRFOn)){
1614 AliInfo("You have choosen to write the RMS method but it is not on!");
1617 if ((fFitPRFNDB == 0) && (!fFitPRFOn)){
1618 AliInfo("You have choosen to write the fit method but it is not on!");
1622 // Number of Xbins (detectors or groups of pads)
1623 if (!InitFit(0,2)) {
1626 fStatisticMean = 0.0;
1628 fNumberFitSuccess = 0;
1632 fCalibraVector = new AliTRDCalibraVector();
1634 // Init fCountDet and fCount
1635 InitfCountDetAndfCount(2);
1636 TGraphErrors *projprftree = 0x0;
1637 tree->SetBranchAddress("histo",&projprftree);
1638 TObjArray *vectorplace = fCalibraVector->ConvertTreeVector(tree);
1640 // Beginning of the loop
1641 for (Int_t idect = fDect1[2]; idect < fDect2[2]; idect++) {
1643 // Search if the group is in the VectorCH
1644 Int_t place = fCalibraVector->SearchInTreeVector(vectorplace,idect);
1647 TH1F *projprf = 0x0;
1652 tree->GetEntry(place);
1653 projprf = CorrectTheError(projprftree);
1656 // Determination of fNnZ, fNnRphi, fNfragZ and fNfragRphi
1657 UpdatefCountDetAndfCount(idect,2);
1659 // Reconstruction of the row and pad group: rowmin, row max ...
1660 ReconstructFitRowMinRowMax(idect,2);
1662 // Rebin and statistic stuff
1663 // This detector has not enough statistics or was off
1664 if ((place == -1) ||
1666 (fEntriesCurrent < fMinEntries))) {
1668 // Fill with the default values
1669 NotEnoughStatistic(idect,2);
1680 // Statistics of the histos fitted
1681 AliInfo(Form("For the group number %d there are %d stats",idect,fEntriesCurrent));
1683 fStatisticMean += fEntriesCurrent;
1685 // Calcul of "real" coef
1686 if ((fDebug == 1) ||
1688 CalculPRFCoefMean(fCountDet[2],(Int_t) (idect - fDect1[2]));
1691 // Method Mean and fit
1692 // idect is egal for fDebug = 0 and 2, only to fill the hist
1694 FitPRF((TH1 *) projprf,(Int_t) (idect - fDect1[2]));
1697 RmsPRF((TH1 *) projprf,(Int_t) (idect - fDect1[2]));
1700 // Visualise the detector for fDebug 3 or 4
1701 // Here is the reconstruction of the pad and row group is used!
1705 // Fill the tree if end of a detector or only the pointer to the branch!!!
1706 FillInfosFit(idect,2);
1716 // No plot, 1 and 4 error plot, 3 and 4 DB plot
1717 if ((fDebug == 1) ||
1721 if ((fDebug == 4) ||
1727 if (fNumberFit > 0) {
1728 AliInfo(Form("There are %d with at least one entries.",fNumberEnt));
1729 AliInfo(Form("%d fits have been proceeded (sucessfully or not...).",fNumberFit));
1730 AliInfo(Form("There is a mean statistic of: %d over these fitted histograms and %d successfulled fits"
1731 ,(Int_t) fStatisticMean/fNumberFit,fNumberFitSuccess));
1732 fStatisticMean = fStatisticMean / fNumberFit;
1735 AliInfo(Form("There are %d with at least one entries. There is no fit!",fNumberEnt));
1738 // Write the things!
1739 if (fWriteCoef[2]) {
1747 //____________Functions for seeing if the pad is really okey___________________
1749 //_____________________________________________________________________________
1750 Bool_t AliTRDCalibraFit::SetModeCalibrationFromTObject(TObject *object, Int_t i)
1753 // Set fNz[i] and fNrphi[i] of the AliTRDCalibraFit::Instance()
1754 // corresponding to the given TObject
1757 const char *nametitle = object->GetTitle();
1760 const Char_t *patternz0 = "Nz0";
1761 const Char_t *patternz1 = "Nz1";
1762 const Char_t *patternz2 = "Nz2";
1763 const Char_t *patternz3 = "Nz3";
1764 const Char_t *patternz4 = "Nz4";
1765 const Char_t *patternrphi0 = "Nrphi0";
1766 const Char_t *patternrphi1 = "Nrphi1";
1767 const Char_t *patternrphi2 = "Nrphi2";
1768 const Char_t *patternrphi3 = "Nrphi3";
1769 const Char_t *patternrphi4 = "Nrphi4";
1770 const Char_t *patternrphi5 = "Nrphi5";
1771 const Char_t *patternrphi6 = "Nrphi6";
1774 UShort_t testrphi = 0;
1777 if (strstr(nametitle,patternz0)) {
1779 fCalibraMode->SetNz(i, 0);
1781 if (strstr(nametitle,patternz1)) {
1783 fCalibraMode->SetNz(i ,1);
1785 if (strstr(nametitle,patternz2)) {
1787 fCalibraMode->SetNz(i ,2);
1789 if (strstr(nametitle,patternz3)) {
1791 fCalibraMode->SetNz(i ,3);
1793 if (strstr(nametitle,patternz4)) {
1795 fCalibraMode->SetNz(i ,4);
1799 if (strstr(nametitle,patternrphi0)) {
1801 fCalibraMode->SetNrphi(i ,0);
1803 if (strstr(nametitle,patternrphi1)) {
1805 fCalibraMode->SetNrphi(i, 1);
1807 if (strstr(nametitle,patternrphi2)) {
1809 fCalibraMode->SetNrphi(i, 2);
1811 if (strstr(nametitle,patternrphi3)) {
1813 fCalibraMode->SetNrphi(i, 3);
1815 if (strstr(nametitle,patternrphi4)) {
1817 fCalibraMode->SetNrphi(i, 4);
1819 if (strstr(nametitle,patternrphi5)) {
1821 fCalibraMode->SetNrphi(i, 5);
1823 if (strstr(nametitle,patternrphi6)) {
1825 fCalibraMode->SetNrphi(i, 6);
1828 // Look if all is okey
1834 fCalibraMode->SetNrphi(i ,0);
1835 fCalibraMode->SetNz(i ,0);
1841 //_____________________________________________________________________________
1842 AliTRDCalDet *AliTRDCalibraFit::CreateDetObjectTree(TTree *tree, Int_t i)
1845 // It creates the AliTRDCalDet object from the tree of the coefficient
1846 // for the calibration i (i != 2)
1847 // It takes the mean value of the coefficients per detector
1848 // This object has to be written in the database
1851 // Create the DetObject
1852 AliTRDCalDet *object = 0x0;
1854 object = new AliTRDCalDet("ChamberGainFactor","GainFactor (detector value)");
1857 object = new AliTRDCalDet("ChamberVdrift","TRD drift velocities (detector value)");
1860 object = new AliTRDCalDet("ChamberT0","T0 (detector value)");
1864 Int_t detector = -1;
1865 Float_t values[2304];
1866 tree->SetBranchAddress("detector",&detector);
1868 tree->SetBranchAddress("gainPad",values);
1871 tree->SetBranchAddress("vdrift" ,values);
1874 tree->SetBranchAddress("t0" ,values);
1877 // For calculating the mean
1880 Int_t numberofentries = tree->GetEntries();
1882 if (numberofentries != 540) {
1883 AliInfo("The tree is not complete");
1886 for (Int_t det = 0; det < numberofentries; ++det) {
1887 tree->GetEntry(det);
1888 if (GetChamber(detector) == 2) {
1896 for (Int_t k = 0; k < nto; k++) {
1897 mean += TMath::Abs(values[k]) / nto;
1901 for (Int_t k = 0; k < nto; k++) {
1902 if(k == 0) mean = values[k];
1903 if(mean > values[k]) mean = values[k];
1906 object->SetValue(detector,mean);
1913 //_____________________________________________________________________________
1914 TObject *AliTRDCalibraFit::CreatePadObjectTree(TTree *tree, Int_t i
1915 , AliTRDCalDet *detobject)
1918 // It Creates the AliTRDCalPad object from the tree of the
1919 // coefficient for the calibration i (i != 2)
1920 // You need first to create the object for the detectors,
1921 // where the mean value is put.
1922 // This object has to be written in the database
1925 // Create the DetObject
1926 AliTRDCalPad *object = 0x0;
1928 object = new AliTRDCalPad("GainFactor","GainFactor (local variations)");
1931 object = new AliTRDCalPad("LocalVdrift","TRD drift velocities (local variations)");
1934 object = new AliTRDCalPad("LocalT0","T0 (local variations)");
1938 Int_t detector = -1;
1939 Float_t values[2304];
1940 tree->SetBranchAddress("detector",&detector);
1942 tree->SetBranchAddress("gainPad",values);
1945 tree->SetBranchAddress("vdrift" ,values);
1948 tree->SetBranchAddress("t0" ,values);
1953 Int_t numberofentries = tree->GetEntries();
1955 if (numberofentries != 540) {
1956 AliInfo("The tree is not complete");
1959 for (Int_t det = 0; det < numberofentries; ++det) {
1960 tree->GetEntry(det);
1961 AliTRDCalROC *calROC = object->GetCalROC(detector);
1962 mean = detobject->GetValue(detector);
1963 if ((mean == 0) && (i != 3)) {
1966 Int_t rowMax = calROC->GetNrows();
1967 Int_t colMax = calROC->GetNcols();
1968 for (Int_t row = 0; row < rowMax; ++row) {
1969 for (Int_t col = 0; col < colMax; ++col) {
1970 if(i != 3) calROC->SetValue(col,row,TMath::Abs(values[(Int_t) (col*rowMax+row)])/mean);
1971 else calROC->SetValue(col,row,values[(Int_t) (col*rowMax+row)]-mean);
1981 //_____________________________________________________________________________
1982 TObject *AliTRDCalibraFit::CreatePadObjectTree(TTree *tree)
1985 // It Creates the AliTRDCalPad object from the tree of the
1986 // coefficient for the calibration PRF (i = 2)
1987 // This object has to be written in the database
1990 // Create the DetObject
1991 AliTRDCalPad *object = new AliTRDCalPad("PRFWidth","PRFWidth");
1994 Int_t detector = -1;
1995 Float_t values[2304];
1996 tree->SetBranchAddress("detector",&detector);
1997 tree->SetBranchAddress("width" ,values);
2000 Int_t numberofentries = tree->GetEntries();
2002 if (numberofentries != 540) {
2003 AliInfo("The tree is not complete");
2006 for (Int_t det = 0; det < numberofentries; ++det) {
2007 tree->GetEntry(det);
2008 AliTRDCalROC *calROC = object->GetCalROC(detector);
2009 Int_t rowMax = calROC->GetNrows();
2010 Int_t colMax = calROC->GetNcols();
2011 for (Int_t row = 0; row < rowMax; ++row) {
2012 for (Int_t col = 0; col < colMax; ++col) {
2013 calROC->SetValue(col,row,TMath::Abs(values[(Int_t) (col*rowMax+row)]));
2022 //_____________________________________________________________________________
2023 void AliTRDCalibraFit::SetPeriodeFitPH(Int_t periodeFitPH)
2026 // Set FitPH if 1 then each detector will be fitted
2029 if (periodeFitPH > 0) {
2030 fFitPHPeriode = periodeFitPH;
2033 AliInfo("periodeFitPH must be higher than 0!");
2038 //_____________________________________________________________________________
2039 void AliTRDCalibraFit::SetFitPRFNDB(Int_t fitPRFNDB)
2042 // TO choose the method that you write into the database
2045 if ((fitPRFNDB >= 3) || (fitPRFNDB == 1)) {
2046 AliInfo("fitPRFNDB is not a correct number!");
2049 fFitPRFNDB = fitPRFNDB;
2054 //_____________________________________________________________________________
2055 void AliTRDCalibraFit::SetFitChargeNDB(Int_t fitChargeNDB)
2058 // To choose the method that you write into the database
2060 if ((fitChargeNDB >= 5) || (fitChargeNDB == 3)) {
2061 AliInfo("fitChargeNDB is not a correct number!");
2064 fFitChargeNDB = fitChargeNDB;
2069 //_____________________________________________________________________________
2070 void AliTRDCalibraFit::SetFitPHNDB(Int_t fitPHNDB)
2073 // To choose the method that you write into the database
2076 if ((fitPHNDB >= 4) || (fitPHNDB == 2)) {
2077 AliInfo("fitPHNDB is not a correct number!");
2080 fFitPHNDB = fitPHNDB;
2085 //_____________________________________________________________________________
2086 void AliTRDCalibraFit::SetBeginFitCharge(Float_t beginFitCharge)
2089 // The fit of the deposited charge distribution begins at
2090 // histo->Mean()/beginFitCharge
2091 // You can here set beginFitCharge
2094 if (beginFitCharge > 0) {
2095 fBeginFitCharge = beginFitCharge;
2098 AliInfo("beginFitCharge must be strict positif!");
2103 //_____________________________________________________________________________
2104 void AliTRDCalibraFit::SetT0Shift(Float_t t0Shift)
2107 // The t0 calculated with the maximum positif slope is shift from t0Shift
2108 // You can here set t0Shift
2115 AliInfo("t0Shift must be strict positif!");
2120 //_____________________________________________________________________________
2121 void AliTRDCalibraFit::SetRangeFitPRF(Float_t rangeFitPRF)
2124 // The fit of the PRF is from -rangeFitPRF to rangeFitPRF
2125 // You can here set rangeFitPRF
2128 if ((rangeFitPRF > 0) &&
2129 (rangeFitPRF <= 1.5)) {
2130 fRangeFitPRF = rangeFitPRF;
2133 AliInfo("rangeFitPRF must be between 0 and 1.0");
2138 //_____________________________________________________________________________
2139 void AliTRDCalibraFit::SetRebin(Short_t rebin)
2142 // Rebin with rebin time less bins the Ch histo
2143 // You can set here rebin that should divide the number of bins of CH histo
2148 AliInfo("You have to be sure that fRebin divides fNumberBinCharge used!");
2151 AliInfo("You have to choose a positiv value!");
2156 //____________Pad Calibration Public___________________________________________
2158 //____________Protected Functions______________________________________________
2159 //____________Create the 2D histo to be filled online__________________________
2161 //____________Fit______________________________________________________________
2162 //____________Create histos if fDebug == 1 or fDebug >= 3______________________
2164 //_____________________________________________________________________________
2165 void AliTRDCalibraFit::InitArrayFitPH()
2168 // Initialise fCoefVdrift[3] and fCoefVdriftE[2] to the right dimension
2171 Int_t nbins = fDect2[1]-fDect1[1];
2173 fCoefVdrift[2] = new Double_t[nbins];
2175 // Init the pointer to nbins
2177 fCoefVdrift[0] = new Double_t[nbins];
2178 fCoefVdriftE[0] = new Double_t[nbins];
2179 for(Int_t k = 0; k < nbins; k++){
2180 fCoefVdriftE[0][k] = 0.0;
2186 fCoefVdrift[1] = new Double_t[nbins];
2187 fCoefVdriftE[1] = new Double_t[nbins];
2188 for(Int_t k = 0; k < nbins; k++){
2189 fCoefVdriftE[1][k] = 0.0;
2193 fCoefVdrift[3] = new Double_t[nbins];
2194 fCoefVdriftE[2] = new Double_t[nbins];
2195 for(Int_t k = 0; k < nbins; k++){
2196 fCoefVdriftE[2][k] = 0.0;
2202 //_____________________________________________________________________________
2203 void AliTRDCalibraFit::InitArrayFitT0()
2206 // Initialise fCoefT0[3] and fCoefT0E[2] to the right dimension
2209 Int_t nbins = fDect2[1]-fDect1[1];
2211 fCoefT0[2] = new Double_t[nbins];
2213 // Init the pointer to nbins
2215 fCoefT0[0] = new Double_t[nbins];
2216 fCoefT0E[0] = new Double_t[nbins];
2217 for(Int_t k = 0; k < nbins; k++){
2218 fCoefT0E[0][k] = 0.0;
2222 fCoefT0[1] = new Double_t[nbins];
2223 fCoefT0E[1] = new Double_t[nbins];
2224 for(Int_t k = 0; k < nbins; k++){
2225 fCoefT0E[1][k] = 0.0;
2229 fCoefT0[3] = new Double_t[nbins];
2230 fCoefT0E[2] = new Double_t[nbins];
2231 for(Int_t k = 0; k < nbins; k++){
2232 fCoefT0E[2][k] = 0.0;
2238 //_____________________________________________________________________________
2239 void AliTRDCalibraFit::InitArrayFitCH()
2242 // Initialise fCoefCharge[4] and fCoefChargeE[3] to the right dimension
2245 Int_t nbins = fDect2[0]-fDect1[0];
2247 //Init the pointer to nbins
2249 fCoefCharge[1] = new Double_t[nbins];
2250 fCoefChargeE[1] = new Double_t[nbins];
2251 for(Int_t k = 0; k < nbins; k++){
2252 fCoefChargeE[1][k] = 0.0;
2256 fCoefCharge[4] = new Double_t[nbins];
2257 fCoefChargeE[3] = new Double_t[nbins];
2258 for(Int_t k = 0; k < nbins; k++){
2259 fCoefChargeE[3][k] = 0.0;
2263 fCoefCharge[0] = new Double_t[nbins];
2264 fCoefChargeE[0] = new Double_t[nbins];
2265 for(Int_t k = 0; k < nbins; k++){
2266 fCoefChargeE[0][k] = 0.0;
2270 if(fFitChargeBisOn){
2271 fCoefCharge[2] = new Double_t[nbins];
2272 fCoefChargeE[2] = new Double_t[nbins];
2273 for(Int_t k = 0; k < nbins; k++){
2274 fCoefChargeE[2][k] = 0.0;
2278 fCoefCharge[3] = new Double_t[nbins];
2282 //_____________________________________________________________________________
2283 void AliTRDCalibraFit::InitArrayFitPRF()
2286 // Initialise fCoefPRF[2] and fCoefPRFE to the right dimension
2289 Int_t nbins = fDect2[2]-fDect1[2];
2290 fCoefPRF[1] = new Double_t[nbins];
2292 //Init the pointer to nbins
2294 fCoefPRF[0] = new Double_t[nbins];
2295 fCoefPRFE[0] = new Double_t[nbins];
2296 for(Int_t k = 0; k < nbins; k++){
2297 fCoefPRFE[0][k] = 0.0;
2301 fCoefPRF[2] = new Double_t[nbins];
2302 fCoefPRFE[1] = new Double_t[nbins];
2303 for(Int_t k = 0; k < nbins; k++){
2304 fCoefPRFE[1][k] = 0.0;
2309 //_____________________________________________________________________________
2310 void AliTRDCalibraFit::CreateFitHistoPRFDB(Int_t rowMax, Int_t colMax)
2313 // Create the histos for fDebug = 3 and fDebug = 4 (Fit functions)
2316 fCoefPRFDB[0] = new TH2F("coefPRF0","",rowMax,0,rowMax,colMax,0,colMax);
2317 fCoefPRFDB[0]->SetStats(0);
2318 fCoefPRFDB[0]->SetXTitle("row Number");
2319 fCoefPRFDB[0]->SetYTitle("col Number");
2320 fCoefPRFDB[0]->SetZTitle("PRF width [pad width units]");
2321 fCoefPRFDB[0]->SetFillColor(6);
2322 fCoefPRFDB[0]->SetLineColor(6);
2325 fCoefPRFDB[1] = new TH2F("coefPRF1","",rowMax,0,rowMax,colMax,0,colMax);
2326 fCoefPRFDB[1]->SetStats(0);
2327 fCoefPRFDB[1]->SetXTitle("row Number");
2328 fCoefPRFDB[1]->SetYTitle("col Number");
2329 fCoefPRFDB[1]->SetZTitle("PRF width [pad width units]");
2330 fCoefPRFDB[1]->SetFillColor(1);
2331 fCoefPRFDB[1]->SetLineColor(1);
2335 //_____________________________________________________________________________
2336 void AliTRDCalibraFit::CreateFitHistoCHDB(Int_t rowMax, Int_t colMax)
2339 // Create the histos for fDebug = 3 and fDebug = 4 (Fit functions)
2343 fCoefChargeDB[0] = new TH2F("coefchargedb0","",rowMax,0,rowMax,colMax,0,colMax);
2344 fCoefChargeDB[0]->SetStats(0);
2345 fCoefChargeDB[0]->SetXTitle("row Number");
2346 fCoefChargeDB[0]->SetYTitle("col Number");
2347 fCoefChargeDB[0]->SetZTitle("f_{g} Fit method");
2348 fCoefChargeDB[0]->SetFillColor(6);
2349 fCoefChargeDB[0]->SetLineColor(6);
2351 if(fFitChargeBisOn){
2352 fCoefChargeDB[2] = new TH2F("coefchargedb2","",rowMax,0,rowMax,colMax,0,colMax);
2353 fCoefChargeDB[2]->SetStats(0);
2354 fCoefChargeDB[2]->SetXTitle("row Number");
2355 fCoefChargeDB[2]->SetYTitle("col Number");
2356 fCoefChargeDB[2]->SetZTitle("f_{g} Fitbis method");
2357 fCoefChargeDB[2]->SetFillColor(8);
2358 fCoefChargeDB[2]->SetLineColor(8);
2361 fCoefChargeDB[3] = new TH2F("coefchargedb3","",rowMax,0,rowMax,colMax,0,colMax);
2362 fCoefChargeDB[3]->SetStats(0);
2363 fCoefChargeDB[3]->SetXTitle("row Number");
2364 fCoefChargeDB[3]->SetYTitle("col Number");
2365 fCoefChargeDB[3]->SetFillColor(1);
2366 fCoefChargeDB[3]->SetLineColor(1);
2370 fCoefChargeDB[1] = new TH2F("coefchargedb1","",rowMax,0,rowMax,colMax,0,colMax);
2371 fCoefChargeDB[1]->SetStats(0);
2372 fCoefChargeDB[1]->SetXTitle("row Number");
2373 fCoefChargeDB[1]->SetYTitle("col Number");
2374 fCoefChargeDB[1]->SetZTitle("f_{g} Mean method");
2375 fCoefChargeDB[1]->SetFillColor(2);
2376 fCoefChargeDB[1]->SetLineColor(2);
2381 //_____________________________________________________________________________
2382 void AliTRDCalibraFit::CreateFitHistoPHDB(Int_t rowMax, Int_t colMax)
2385 // Create the histos for fDebug = 3 and fDebug = 4 (Fit functions)
2389 fCoefVdriftDB[0] = new TH2F("coefvdriftdb0","",rowMax,0,rowMax,colMax,0,colMax);
2390 fCoefVdriftDB[0]->SetStats(0);
2391 fCoefVdriftDB[0]->SetXTitle("row Number");
2392 fCoefVdriftDB[0]->SetYTitle("col Number");
2393 fCoefVdriftDB[0]->SetZTitle("v_{drift} Fit method");
2394 fCoefVdriftDB[0]->SetFillColor(6);
2395 fCoefVdriftDB[0]->SetLineColor(6);
2399 fCoefVdriftDB[1] = new TH2F("coefvdriftdb1","",rowMax,0,rowMax,colMax,0,colMax);
2400 fCoefVdriftDB[1]->SetStats(0);
2401 fCoefVdriftDB[1]->SetXTitle("row Number");
2402 fCoefVdriftDB[1]->SetYTitle("col Number");
2403 fCoefVdriftDB[1]->SetZTitle("v_{drift} slope method");
2404 fCoefVdriftDB[1]->SetFillColor(2);
2405 fCoefVdriftDB[1]->SetLineColor(2);
2408 fCoefVdriftDB[2] = new TH2F("coefvdriftdb1","",rowMax,0,rowMax,colMax,0,colMax);
2409 fCoefVdriftDB[2]->SetStats(0);
2410 fCoefVdriftDB[2]->SetXTitle("row Number");
2411 fCoefVdriftDB[2]->SetYTitle("col Number");
2412 fCoefVdriftDB[2]->SetZTitle("v_{drift} slope method");
2413 fCoefVdriftDB[2]->SetFillColor(1);
2414 fCoefVdriftDB[2]->SetLineColor(1);
2419 //_____________________________________________________________________________
2420 void AliTRDCalibraFit::CreateFitHistoT0DB(Int_t rowMax, Int_t colMax)
2423 // Create the histos for fDebug = 3 and fDebug = 4 (Fit functions)
2427 fCoefT0DB[0] = new TH2F("coefT0db0","",rowMax,0,rowMax,colMax,0,colMax);
2428 fCoefT0DB[0]->SetStats(0);
2429 fCoefT0DB[0]->SetXTitle("row Number");
2430 fCoefT0DB[0]->SetYTitle("col Number");
2431 fCoefT0DB[0]->SetZTitle("t0 Fit method");
2432 fCoefT0DB[0]->SetFillColor(6);
2433 fCoefT0DB[0]->SetLineColor(6);
2436 fCoefT0DB[1] = new TH2F("coefT0db1","",rowMax,0,rowMax,colMax,0,colMax);
2437 fCoefT0DB[1]->SetStats(0);
2438 fCoefT0DB[1]->SetXTitle("row Number");
2439 fCoefT0DB[1]->SetYTitle("col Number");
2440 fCoefT0DB[1]->SetZTitle("t0 slope method");
2441 fCoefT0DB[1]->SetFillColor(2);
2442 fCoefT0DB[1]->SetLineColor(2);
2445 fCoefT0DB[2] = new TH2F("coefT0db1","",rowMax,0,rowMax,colMax,0,colMax);
2446 fCoefT0DB[2]->SetStats(0);
2447 fCoefT0DB[2]->SetXTitle("row Number");
2448 fCoefT0DB[2]->SetYTitle("col Number");
2449 fCoefT0DB[2]->SetZTitle("t0 slope method");
2450 fCoefT0DB[2]->SetFillColor(1);
2451 fCoefT0DB[2]->SetLineColor(1);
2456 //_____________________________________________________________________________
2457 Bool_t AliTRDCalibraFit::FillVectorFitCH(Int_t countdet)
2460 // For the Fit functions fill the vector FitCH special for the gain calibration
2463 AliTRDFitCHInfo *fitCHInfo = new AliTRDFitCHInfo();
2466 if (GetChamber(countdet) == 2) {
2473 //printf("For the detector %d , ntotal %d and fCoefCH[0] %f\n",countdet,ntotal,fCoefCH[0]);
2474 Float_t *coef = new Float_t[ntotal];
2475 for (Int_t i = 0; i < ntotal; i++) {
2476 coef[i] = fCoefCH[i];
2479 Int_t detector = countdet;
2481 fitCHInfo->SetCoef(coef);
2482 fitCHInfo->SetDetector(detector);
2483 fVectorFitCH->Add((TObject *) fitCHInfo);
2489 //____________Functions for initialising the AliTRDCalibraFit in the code_________
2490 Bool_t AliTRDCalibraFit::InitFit(Int_t nbins, Int_t i)
2493 // Init the calibration mode (Nz, Nrphi), the histograms for
2494 // debugging the fit methods if fDebug > 0,
2497 gStyle->SetPalette(1);
2498 gStyle->SetOptStat(1111);
2499 gStyle->SetPadBorderMode(0);
2500 gStyle->SetCanvasColor(10);
2501 gStyle->SetPadLeftMargin(0.13);
2502 gStyle->SetPadRightMargin(0.01);
2504 // Mode groups of pads: the total number of bins!
2505 Int_t numberofbinsexpected = 0;
2506 fCalibraMode->ModePadCalibration(2,i);
2507 fCalibraMode->ModePadFragmentation(0,2,0,i);
2508 fCalibraMode->SetDetChamb2(i);
2510 AliInfo(Form("For the chamber 2: %d",fCalibraMode->GetDetChamb2(i)));
2512 numberofbinsexpected += 6 * 18 * fCalibraMode->GetDetChamb2(i);
2513 fCalibraMode->ModePadCalibration(0,i);
2514 fCalibraMode->ModePadFragmentation(0,0,0,i);
2515 fCalibraMode->SetDetChamb0(i);
2517 AliInfo(Form("For the other chamber 0: %d",fCalibraMode->GetDetChamb0(i)));
2519 numberofbinsexpected += 6 * 4 * 18 * fCalibraMode->GetDetChamb0(i);
2521 // Quick verification that we have the good pad calibration mode if 2D histos!
2523 if (numberofbinsexpected != nbins) {
2524 AliInfo("It doesn't correspond to the mode of pad group calibration!");
2529 // Security for fDebug 3 and 4
2530 if ((fDebug >= 3) &&
2534 AliInfo("This detector doesn't exit!");
2538 // Determine fDet1 and fDet2
2542 fDect1[i] = fFitVoir;
2543 fDect2[i] = fDect1[i] +1;
2547 fDect2[i] = numberofbinsexpected;
2550 fCalibraMode->CalculXBins(AliTRDgeometry::GetDetector(fDet[0],fDet[1],fDet[2]),i);
2551 fDect1[i] = fCalibraMode->GetXbins(i);
2552 fCalibraMode->CalculXBins((AliTRDgeometry::GetDetector(fDet[0],fDet[1],fDet[2])+1),i);
2553 fDect2[i] = fCalibraMode->GetXbins(i);
2556 // Create the histos for debugging
2561 // Init the VectorFitCH
2562 fVectorFitCH = new TObjArray();
2563 fCoefCH = new Float_t[2304];
2564 for (Int_t k = 0; k < 2304; k++) {
2567 fScaleFitFactor = 0.0;
2569 // Number of Xbins(detectors or groups of pads) if Vector2d
2570 // Quick verification that we are not out of range!
2571 if (fCalibraVector) {
2573 (fCalibraVector->GetVectorCH()->GetEntriesFast() > 0) &&
2574 ((Int_t) fCalibraVector->GetPlaCH()->GetEntriesFast() > 0)) {
2575 if ((Int_t) fCalibraVector->GetVectorCH()->GetEntriesFast() > numberofbinsexpected) {
2576 AliInfo("ch doesn't correspond to the mode of pad group calibration!");
2579 if ((Int_t) fCalibraVector->GetVectorCH()->GetEntriesFast() !=
2580 (Int_t) fCalibraVector->GetPlaCH()->GetEntriesFast()) {
2581 AliInfo("VectorCH doesn't correspond to PlaCH!");
2588 // Debugging: Create the histos
2591 // fDebug == 0 nothing
2598 // fDebug == 2 and fFitVoir no histo
2600 if (fFitVoir < numberofbinsexpected) {
2601 AliInfo(Form("We will see the fit of the object %d",fFitVoir));
2604 AliInfo("fFitVoir is out of range of the histo!");
2609 // fDebug == 3 or 4 and fDet
2611 if ((fCalibraMode->GetNz(0) == 0) && (fCalibraMode->GetNrphi(0) == 0)) {
2612 AliInfo("Do you really want to see one detector without pad groups?");
2616 AliInfo(Form("You will see the detector: iPlane %d, iChamb %d, iSect %d"
2617 ,fDet[0],fDet[1],fDet[2]));
2618 // A little geometry:
2619 Int_t rowMax = fGeo->GetRowMax(fDet[0],fDet[1],fDet[2]);
2620 Int_t colMax = fGeo->GetColMax(fDet[0]);
2621 // Create the histos to visualise
2622 CreateFitHistoCHDB(rowMax,colMax);
2634 // Number of Xbins (detectors or groups of pads) if vector2d
2635 // Quick verification that we are not out of range!
2636 if (fCalibraVector) {
2638 (fCalibraVector->GetVectorPH()->GetEntriesFast() > 0) &&
2639 ((Int_t) fCalibraVector->GetPlaPH()->GetEntriesFast() > 0)) {
2640 if ((Int_t) fCalibraVector->GetVectorPH()->GetEntriesFast() > numberofbinsexpected) {
2641 AliInfo("ph doesn't correspond to the mode of pad group calibration!");
2644 if ((Int_t) fCalibraVector->GetVectorPH()->GetEntriesFast() !=
2645 (Int_t) fCalibraVector->GetPlaPH()->GetEntriesFast()) {
2646 AliInfo("VectorPH doesn't correspond to PlaPH!");
2657 // Debugging: Create the histos
2660 // fDebug == 0 nothing
2664 // Create the histos replique de ph
2669 // fDebug == 2 and fFitVoir no histo
2671 if (fFitVoir < numberofbinsexpected) {
2672 AliInfo(Form("We will see the fit of the object %d",fFitVoir));
2675 AliInfo("fFitVoir is out of range of the histo!");
2680 // fDebug == 3 or 4 and fDet
2682 if ((fCalibraMode->GetNz(1) == 0) &&
2683 (fCalibraMode->GetNrphi(1) == 0)) {
2684 AliInfo("Do you really want to see one detector without pad groups?");
2688 AliInfo(Form("You will see the detector: iPlane %d, iChamb %d, iSect %d"
2689 ,fDet[0],fDet[1],fDet[2]));
2690 // A little geometry:
2691 Int_t rowMax = fGeo->GetRowMax(fDet[0],fDet[1],fDet[2]);
2692 Int_t colMax = fGeo->GetColMax(fDet[0]);
2693 // Create the histos to visualise
2694 CreateFitHistoPHDB(rowMax,colMax);
2695 CreateFitHistoT0DB(rowMax,colMax);
2708 // Number of Xbins(detectors or groups of pads) if vector2d
2709 if (fCalibraVector){
2711 (fCalibraVector->GetVectorPRF()->GetEntriesFast() > 0) &&
2712 (fCalibraVector->GetPlaPRF()->GetEntriesFast() > 0)) {
2713 // Quick verification that we are not out of range!
2714 if ((Int_t) fCalibraVector->GetVectorPRF()->GetEntriesFast() > numberofbinsexpected) {
2715 AliInfo("ch doesn't correspond to the mode of pad group calibration!");
2718 if ((Int_t) fCalibraVector->GetVectorPRF()->GetEntriesFast() !=
2719 (Int_t) fCalibraVector->GetPlaPRF()->GetEntriesFast()) {
2720 AliInfo("VectorPRF doesn't correspond to PlaCH!");
2730 // Debugging: Create the histos
2733 // fDebug == 0 nothing
2737 // Create the histos replique de ch
2741 // fDebug == 2 and fFitVoir no histo
2743 if (fFitVoir < numberofbinsexpected) {
2744 AliInfo(Form("We will see the fit of the object %d",fFitVoir));
2747 AliInfo("fFitVoir is out of range of the histo!");
2752 // fDebug == 3 or 4 and fDet
2754 if ((fCalibraMode->GetNz(2) == 0) &&
2755 (fCalibraMode->GetNrphi(2) == 0)) {
2756 AliInfo("Do you really want to see one detector without pad groups?");
2760 AliInfo(Form("You will see the detector: iPlane %d, iChamb %d, iSect %d"
2761 ,fDet[0],fDet[1],fDet[2]));
2762 // A little geometry:
2763 Int_t rowMax = fGeo->GetRowMax(fDet[0],fDet[1],fDet[2]);
2764 Int_t colMax = fGeo->GetColMax(fDet[0]);
2765 // Create the histos to visualise
2766 CreateFitHistoPRFDB(rowMax,colMax);
2779 //____________Functions for initialising the AliTRDCalibraFit in the code_________
2780 void AliTRDCalibraFit::InitfCountDetAndfCount(Int_t i)
2783 // Init the current detector where we are fCountDet and the
2784 // next fCount for the functions Fit...
2787 // Loop on the Xbins of ch!!
2788 fCountDet[i] = -1; // Current detector
2789 fCount[i] = 0; // To find the next detector
2794 // Set countdet to the detector
2795 fCountDet[i] = AliTRDgeometry::GetDetector(fDet[0],fDet[1],fDet[2]);
2797 // Determination of fNnZ, fNnRphi, fNfragZ and fNfragRphi
2798 fCalibraMode->ModePadCalibration(fDet[1],i);
2799 fCalibraMode->ModePadFragmentation(fDet[0],fDet[1],fDet[2],i);
2801 // Set counter to write at the end of the detector
2802 fCount[i] = fDect1[i] + fCalibraMode->GetNfragZ(i)*fCalibraMode->GetNfragRphi(i);
2808 //____________Functions for initialising the AliTRDCalibraFit in the code_________
2809 void AliTRDCalibraFit::UpdatefCountDetAndfCount(Int_t idect, Int_t i)
2812 // See if we are in a new detector and update the
2813 // variables fNfragZ and fNfragRphi if yes
2816 // Determination of fNnZ, fNnRphi, fNfragZ and fNfragRphi
2817 // If fDebug == 1 or 0
2818 if ((fDebug == 0) ||
2821 if (fCount[i] == idect) {
2823 // On en est au detector
2826 // Determination of fNnZ, fNnRphi, fNfragZ and fNfragRphi
2827 fCalibraMode->ModePadCalibration((Int_t) GetChamber(fCountDet[i]),i);
2828 fCalibraMode->ModePadFragmentation((Int_t) GetPlane(fCountDet[i])
2829 ,(Int_t) GetChamber(fCountDet[i])
2830 ,(Int_t) GetSector(fCountDet[i]),i);
2832 // Set for the next detector
2833 fCount[i] += fCalibraMode->GetNfragZ(i)*fCalibraMode->GetNfragRphi(i);
2841 //____________Functions for initialising the AliTRDCalibraFit in the code_________
2842 void AliTRDCalibraFit::ReconstructFitRowMinRowMax(Int_t idect, Int_t i)
2845 // Reconstruct the min pad row, max pad row, min pad col and
2846 // max pad col of the calibration group for the Fit functions
2850 fCalibraMode->ReconstructionRowPadGroup((Int_t) (idect-(fCount[i]-(fCalibraMode->GetNfragZ(i)
2851 *fCalibraMode->GetNfragRphi(i)))),i);
2854 fCalibraMode->ReconstructionRowPadGroup((Int_t) (idect-fDect1[i]),i);
2859 //____________Functions for initialising the AliTRDCalibraFit in the code_________
2860 Bool_t AliTRDCalibraFit::NotEnoughStatistic(Int_t idect, Int_t i)
2863 // For the case where there are not enough entries in the histograms
2864 // of the calibration group, the value present in the choosen database
2865 // will be put. A negativ sign enables to know that a fit was not possible.
2869 AliTRDcalibDB *cal = AliTRDcalibDB::Instance();
2871 AliInfo("Could not get calibDB");
2876 AliInfo(Form("The element %d in this detector %d has not enough statistic to be fitted"
2877 ,idect-(fCount[i]-(fCalibraMode->GetNfragZ(i)*fCalibraMode->GetNfragRphi(i))),fCountDet[i]));
2880 AliInfo("The element has not enough statistic to be fitted");
2883 if ((i == 0) && (fDebug != 2)) {
2885 // Calcul the coef from the database choosen
2886 CalculChargeCoefMean(fCountDet[0],(Int_t) (idect-fDect1[0]),kFALSE);
2888 // Fill the coefCH[2304] with negative value to say: not fitted
2889 AliInfo(Form("The row min %d, the row max %d, the colmin %d and the col max %d"
2890 ,fCalibraMode->GetRowMin(0)
2891 ,fCalibraMode->GetRowMax(0)
2892 ,fCalibraMode->GetColMin(0)
2893 ,fCalibraMode->GetColMax(0)));
2894 for (Int_t k = fCalibraMode->GetRowMin(0); k < fCalibraMode->GetRowMax(0); k++) {
2895 for (Int_t j = fCalibraMode->GetColMin(0); j < fCalibraMode->GetColMax(0); j++) {
2896 if (GetChamber(fCountDet[0]) == 2) {
2897 fCoefCH[(Int_t)(j*12+k)] = -TMath::Abs(fChargeCoef[3]);
2899 if (GetChamber(fCountDet[0]) != 2) {
2900 fCoefCH[(Int_t)(j*16+k)] = -TMath::Abs(fChargeCoef[3]);
2905 // Put the default value negative
2906 if ((fDebug == 1) ||
2909 if (fFitChargeBisOn) {
2910 fCoefCharge[2][idect-fDect1[0]]=-TMath::Abs(fChargeCoef[3]);
2912 if (fMeanChargeOn) {
2913 fCoefCharge[1][idect-fDect1[0]]=-TMath::Abs(fChargeCoef[3]);
2916 fCoefCharge[0][idect-fDect1[0]]=-TMath::Abs(fChargeCoef[3]);
2921 // End of one detector
2922 if ((idect == (fCount[0]-1))) {
2923 FillVectorFitCH((Int_t) fCountDet[0]);
2925 for (Int_t k = 0; k < 2304; k++) {
2932 if ((i == 1) && (fDebug != 2)) {
2934 CalculVdriftCoefMean(fCountDet[1],(Int_t) (idect-fDect1[1]));
2935 CalculT0CoefMean(fCountDet[1],(Int_t) (idect-fDect1[1]));
2937 // Put the default value (time0 can be negativ, so we stay with + )
2938 if ((fDebug == 1) ||
2942 fCoefVdrift[0][(idect-fDect1[1])] = -fVdriftCoef[2];
2943 fCoefT0[0][(idect-fDect1[1])] = fT0Coef[2];
2947 fCoefVdrift[1][(idect-fDect1[1])] = -fVdriftCoef[2];
2948 fCoefT0[1][(idect-fDect1[1])] = fT0Coef[2];
2951 fCoefVdrift[3][(idect-fDect1[1])] = -fVdriftCoef[2];
2952 fCoefT0[3][(idect-fDect1[1])] = fT0Coef[2];
2957 // Put the default value
2960 fVdriftCoef[0] = fVdriftCoef[2];
2961 fT0Coef[0] = fT0Coef[2];
2964 fVdriftCoef[1] = fVdriftCoef[2];
2965 fT0Coef[1] = fT0Coef[2];
2968 fVdriftCoef[3] = fVdriftCoef[2];
2969 fT0Coef[3] = fT0Coef[2];
2975 // Fill the tree if end of a detector.
2976 // The pointer to the branch stays with the default value negative!!!
2978 // Pointer to the branch
2979 for (Int_t k = fCalibraMode->GetRowMin(1); k < fCalibraMode->GetRowMax(1); k++) {
2980 for (Int_t j = fCalibraMode->GetColMin(1); j < fCalibraMode->GetColMax(1); j++) {
2981 if (GetChamber(fCountDet[1]) == 2) {
2982 fVdriftPad[(Int_t)(j*12+k)] = -TMath::Abs(fVdriftCoef[2]);
2984 if (GetChamber(fCountDet[1]) != 2) {
2985 fVdriftPad[(Int_t)(j*16+k)] = -TMath::Abs(fVdriftCoef[2]);
2990 // End of one detector
2991 if ((idect == (fCount[1]-1)) && (fDebug != 2)) {
2992 FillTreeVdrift((Int_t) fCountDet[1]);
2996 // Fill the tree if end of a detector.
2997 // The pointer to the branch stays with the default value positive!!!
2998 // Pointer to the branch
2999 for (Int_t k = fCalibraMode->GetRowMin(1); k < fCalibraMode->GetRowMax(1); k++) {
3000 for (Int_t j = fCalibraMode->GetColMin(1); j < fCalibraMode->GetColMax(1); j++) {
3001 if (GetChamber(fCountDet[1]) == 2) {
3002 fT0Pad[(Int_t)(j*12+k)] = fT0Coef[2];
3004 if (GetChamber(fCountDet[1]) != 2) {
3005 fT0Pad[(Int_t)(j*16+k)] = fT0Coef[2];
3010 // End of one detector
3011 if ((idect == (fCount[1]-1)) && (fDebug != 2)) {
3012 FillTreeT0((Int_t) fCountDet[1]);
3017 if ((i == 2) && (fDebug != 2)) {
3019 CalculPRFCoefMean(fCountDet[2],(Int_t) (idect-fDect1[2]));
3021 if ((fDebug == 1) ||
3024 fCoefPRF[0][(idect-fDect1[2])] = -fPRFCoef[1];
3027 fCoefPRF[2][(idect-fDect1[2])] = -fPRFCoef[1];
3033 fPRFCoef[0] = fPRFCoef[1];
3036 fPRFCoef[2] = fPRFCoef[1];
3041 // Fill the tree if end of a detector.
3042 // The pointer to the branch stays with the default value 1.5!!!
3043 // Pointer to the branch
3044 for (Int_t k = fCalibraMode->GetRowMin(2); k < fCalibraMode->GetRowMax(2); k++) {
3045 for (Int_t j = fCalibraMode->GetColMin(2); j < fCalibraMode->GetColMax(2); j++) {
3046 if((fGeo->GetColMax(GetPlane(fCountDet[2])) != (j+1)) && (j != 0)){
3047 if (GetChamber(fCountDet[2]) == 2) {
3048 fPRFPad[(Int_t)(j*12+k)] = -fPRFCoef[1];
3050 if (GetChamber(fCountDet[2]) != 2) {
3051 fPRFPad[(Int_t)(j*16+k)] = -fPRFCoef[1];
3056 if (GetChamber(fCountDet[2]) == 2) {
3057 fPRFPad[(Int_t)(j*12+k)] = -((Float_t) cal->GetPRFWidth(fCountDet[2],j,k));
3059 if (GetChamber(fCountDet[2]) != 2) {
3060 fPRFPad[(Int_t)(j*16+k)] = -((Float_t) cal->GetPRFWidth(fCountDet[2],j,k));
3064 if (GetChamber(fCountDet[2]) == 2) {
3065 fPRFPad[(Int_t)(j*12+k)] = -((Float_t) GetPRFDefault(GetPlane(fCountDet[2])));
3067 if (GetChamber(fCountDet[2]) != 2) {
3068 fPRFPad[(Int_t)(j*16+k)] = -((Float_t) GetPRFDefault(GetPlane(fCountDet[2])));
3075 // End of one detector
3076 if ((idect == (fCount[2]-1)) && (fDebug != 2)) {
3077 FillTreePRF((Int_t) fCountDet[2]);
3086 //____________Functions for initialising the AliTRDCalibraFit in the code_________
3087 Bool_t AliTRDCalibraFit::FillInfosFit(Int_t idect, Int_t i)
3090 // Fill the coefficients found with the fits or other
3091 // methods from the Fit functions
3095 AliTRDcalibDB *cal = AliTRDcalibDB::Instance();
3097 AliInfo("Could not get calibDB");
3101 if ((i == 0) && (fDebug != 2)) {
3102 // Fill the coefCH[2304] with fChargeCoef[0]
3103 // that would be negativ only if the fit failed totally
3104 //printf("for fCountDet %d we have %f\n",fCountDet[0],fChargeCoef[fFitChargeNDB]);
3105 //printf("RowMin %d RowMax %d ColMin %d ColMax %d\n",fCalibraMode->GetRowMin(0),fCalibraMode->GetRowMax(0),fCalibraMode->GetColMin(0),fCalibraMode->GetColMax(0));
3106 for (Int_t k = fCalibraMode->GetRowMin(0); k < fCalibraMode->GetRowMax(0); k++) {
3107 for (Int_t j = fCalibraMode->GetColMin(0); j < fCalibraMode->GetColMax(0); j++) {
3108 if (GetChamber(fCountDet[0]) == 2) {
3109 fCoefCH[(Int_t)(j*12+k)] = fChargeCoef[fFitChargeNDB];
3111 if (GetChamber(fCountDet[0]) != 2) {
3112 fCoefCH[(Int_t)(j*16+k)] = fChargeCoef[fFitChargeNDB];
3116 // End of one detector
3117 if ((idect == (fCount[0]-1))) {
3118 FillVectorFitCH((Int_t) fCountDet[0]);
3120 for (Int_t k = 0; k < 2304; k++) {
3126 if ((i == 1) && (fDebug != 2)) {
3129 // Pointer to the branch: fVdriftCoef[1] will ne negativ only if the fit failed totally
3130 for (Int_t k = fCalibraMode->GetRowMin(1); k < fCalibraMode->GetRowMax(1); k++) {
3131 for (Int_t j = fCalibraMode->GetColMin(1); j < fCalibraMode->GetColMax(1); j++) {
3132 if (GetChamber(fCountDet[1]) == 2) {
3133 fVdriftPad[(Int_t)(j*12+k)]=fVdriftCoef[fFitPHNDB];
3135 if (GetChamber(fCountDet[1]) != 2) {
3136 fVdriftPad[(Int_t)(j*16+k)]=fVdriftCoef[fFitPHNDB];
3140 // End of one detector
3141 if ((idect == (fCount[1]-1)) && (fDebug != 2)) {
3142 FillTreeVdrift((Int_t) fCountDet[1]);
3146 // Pointer to the branch: fT0Coef[1] will ne negativ only if the fit failed totally
3147 for (Int_t k = fCalibraMode->GetRowMin(1); k < fCalibraMode->GetRowMax(1); k++) {
3148 for (Int_t j = fCalibraMode->GetColMin(1); j < fCalibraMode->GetColMax(1); j++) {
3149 if (GetChamber(fCountDet[1]) == 2) {
3150 fT0Pad[(Int_t)(j*12+k)]=fT0Coef[fFitPHNDB];
3152 if (GetChamber(fCountDet[1]) != 2) {
3153 fT0Pad[(Int_t)(j*16+k)]=fT0Coef[fFitPHNDB];
3157 // End of one detector
3158 if ((idect == (fCount[1]-1)) && (fDebug != 2)) {
3159 FillTreeT0((Int_t) fCountDet[1]);
3164 if ((i == 2) && (fDebug != 2)) {
3165 // Pointer to the branch
3166 for (Int_t k = fCalibraMode->GetRowMin(2); k < fCalibraMode->GetRowMax(2); k++) {
3167 for (Int_t j = fCalibraMode->GetColMin(2); j < fCalibraMode->GetColMax(2); j++) {
3168 if ((fGeo->GetColMax(GetPlane(fCountDet[2])) != (j+1)) && (j != 0)) {
3169 if (GetChamber(fCountDet[2]) == 2) {
3170 fPRFPad[(Int_t)(j*12+k)] = fPRFCoef[fFitPRFNDB];
3172 if (GetChamber(fCountDet[2]) != 2) {
3173 fPRFPad[(Int_t)(j*16+k)] = fPRFCoef[fFitPRFNDB];
3178 if (GetChamber(fCountDet[2]) == 2) {
3179 fPRFPad[(Int_t)(j*12+k)] = (Float_t) cal->GetPRFWidth(fCountDet[2],j,k);
3181 if (GetChamber(fCountDet[2]) != 2) {
3182 fPRFPad[(Int_t)(j*16+k)] = (Float_t) cal->GetPRFWidth(fCountDet[2],j,k);
3186 if (GetChamber(fCountDet[2]) == 2) {
3187 fPRFPad[(Int_t)(j*12+k)] = (Float_t) GetPRFDefault(GetPlane(fCountDet[2]));
3189 if (GetChamber(fCountDet[2]) != 2) {
3190 fPRFPad[(Int_t)(j*16+k)] = (Float_t) GetPRFDefault(GetPlane(fCountDet[2]));
3196 // End of one detector
3197 if ((idect == (fCount[2]-1)) && (fDebug != 2)) {
3198 FillTreePRF((Int_t) fCountDet[2]);
3206 //____________Functions for initialising the AliTRDCalibraFit in the code_________
3207 Bool_t AliTRDCalibraFit::WriteFitInfos(Int_t i)
3210 // In the case the user wants to write a file with a tree of the found
3211 // coefficients for the calibration before putting them in the database
3214 TFile *fout = TFile::Open(fWriteNameCoef,"UPDATE");
3215 // Check if the file could be opened
3216 if (!fout || !fout->IsOpen()) {
3217 AliInfo("No File found!");
3221 if ((i == 0) && (fDebug != 2)) {
3223 if ((fDebug == 4) ||
3228 fout->WriteTObject(fGain,fGain->GetName(),(Option_t *) "writedelete");
3231 if ((i == 1) && (fDebug != 2)) {
3233 if ((fDebug == 4) ||
3238 fout->WriteTObject(fVdrift,fVdrift->GetName(),(Option_t *) "writedelete");
3240 if ((fDebug == 4) ||
3245 fout->WriteTObject(fT0,fT0->GetName(),(Option_t *) "writedelete");
3248 if ((i == 2) && (fDebug != 2)) {
3250 if ((fDebug == 4) ||
3255 fout->WriteTObject(fPRF,fPRF->GetName(),(Option_t *) "writedelete");
3265 //____________Fill Coef DB in case of visualisation of one detector____________
3268 //_____________________________________________________________________________
3269 void AliTRDCalibraFit::FillCoefVdriftDB()
3272 // Fill the histos for fDebug = 3 and fDebug = 4 to visualise the detector
3275 for (Int_t row = fCalibraMode->GetRowMin(1); row < fCalibraMode->GetRowMax(1); row++) {
3276 for (Int_t col = fCalibraMode->GetColMin(1); col < fCalibraMode->GetColMax(1); col++) {
3278 fCoefVdriftDB[1]->SetBinContent(row+1,col+1,TMath::Abs(fVdriftCoef[1]));
3281 fCoefVdriftDB[0]->SetBinContent(row+1,col+1,TMath::Abs(fVdriftCoef[0]));
3283 if (fFitLagrPolOn ) {
3284 fCoefVdriftDB[2]->SetBinContent(row+1,col+1,TMath::Abs(fVdriftCoef[3]));
3291 //_____________________________________________________________________________
3292 void AliTRDCalibraFit::FillCoefT0DB()
3295 // Fill the histos for fDebug = 3 and fDebug = 4 to visualise the detector
3298 for (Int_t row = fCalibraMode->GetRowMin(1); row < fCalibraMode->GetRowMax(1); row++) {
3299 for (Int_t col = fCalibraMode->GetColMin(1); col < fCalibraMode->GetColMax(1); col++) {
3301 fCoefT0DB[1]->SetBinContent(row+1,col+1,TMath::Abs(fT0Coef[1]));
3304 fCoefT0DB[0]->SetBinContent(row+1,col+1,TMath::Abs(fT0Coef[0]));
3306 if (fFitLagrPolOn) {
3307 fCoefT0DB[2]->SetBinContent(row+1,col+1,TMath::Abs(fT0Coef[3]));
3314 //_____________________________________________________________________________
3315 void AliTRDCalibraFit::FillCoefChargeDB()
3318 // Fill the histos for fDebug = 3 and fDebug = 4 to visualise the detector
3321 for (Int_t row = fCalibraMode->GetRowMin(0); row < fCalibraMode->GetRowMax(0); row++) {
3322 for (Int_t col = fCalibraMode->GetColMin(0); col < fCalibraMode->GetColMax(0); col++) {
3323 if (fMeanChargeOn) {
3324 fCoefChargeDB[1]->SetBinContent(row+1,col+1,TMath::Abs(fChargeCoef[1]));
3326 if (fFitChargeBisOn) {
3327 fCoefChargeDB[2]->SetBinContent(row+1,col+1,TMath::Abs(fChargeCoef[2]));
3330 fCoefChargeDB[0]->SetBinContent(row+1,col+1,TMath::Abs(fChargeCoef[0]));
3333 fCoefChargeDB[3]->SetBinContent(row+1,col+1,TMath::Abs(fChargeCoef[4]));
3340 //_____________________________________________________________________________
3341 void AliTRDCalibraFit::FillCoefPRFDB()
3344 // Fill the histos for fDebug = 3 and fDebug = 4 to visualise the detector
3347 for (Int_t row = fCalibraMode->GetRowMin(2); row < fCalibraMode->GetRowMax(2); row++) {
3348 for (Int_t col = fCalibraMode->GetColMin(2); col < fCalibraMode->GetColMax(2); col++) {
3349 fCoefPRFDB[0]->SetBinContent(row+1,col+1,fPRFCoef[0]);
3354 for (Int_t row = fCalibraMode->GetRowMin(2); row < fCalibraMode->GetRowMax(2); row++) {
3355 for (Int_t col = fCalibraMode->GetColMin(2); col < fCalibraMode->GetColMax(2); col++) {
3356 fCoefPRFDB[1]->SetBinContent(row+1,col+1,fPRFCoef[2]);
3364 //____________Plot histos CoefPRF....__________________________________________
3367 //_____________________________________________________________________________
3368 void AliTRDCalibraFit::PlotWriteCH()
3371 // Scale the coefficients to one, create the graph errors and write them if wanted
3374 //TObjArray of the grapherrors and so on
3375 TObjArray *listofgraphs = new TObjArray();
3377 Int_t nbins = fDect2[0]-fDect1[0];
3380 // We will check fScaleFitFactor for the fFitChargeNDB, otherwise we calculate and normalise to 1
3381 // It can be that fScaleFitFactor is different from scale if we have taken a no default database as reference
3386 counter[0] = 0; //how many groups are fitted for 0
3387 counter[1] = 0; //how many groups are with mean for 1
3388 counter[2] = 0; //how many groups are fitted for 2
3389 counter[3] = 0; //how many groups are fitted for 4
3391 Double_t scale = 1.0;
3394 // Is -1 if no fit or mean, is 1 if fit or mean
3395 Double_t *xValuesFitted = new Double_t[nbins];
3396 Double_t *xValuesFittedMean = new Double_t[nbins];
3397 Double_t *xValuesFittedBis = new Double_t[nbins];
3398 Double_t *xValuesFittedMeanW = new Double_t[nbins];
3399 for(Int_t k = 0; k < nbins; k ++){
3400 xValuesFitted[k] = -1;
3401 xValuesFittedMean[k] = -1;
3402 xValuesFittedMeanW[k] = -1;
3403 xValuesFittedBis[k] = -1;
3408 for(Int_t l = 0; l < nbins; l++){
3409 if(fCoefCharge[0][l] > 0){
3410 sum += fCoefCharge[0][l];
3411 xValuesFitted[counter[0]]= l;
3416 if(sum > 0.0) scale = counter[0]/sum;
3417 if(fFitChargeNDB == 0){
3418 if(scale != fScaleFitFactor){
3419 AliInfo(Form("The normalisation is different from a nomalisation to one."));
3420 AliInfo(Form("For one we have %f and here %f",scale,fScaleFitFactor));
3422 AliInfo(Form("It is not normal because we didn't choose a reference database!"));
3425 scale = fScaleFitFactor;
3427 for(Int_t l = 0; l < nbins; l++){
3428 if(fCoefCharge[0][l] > 0){
3429 fCoefCharge[0][l]=fCoefCharge[0][l]*scale;
3430 fCoefChargeE[0][l]=fCoefChargeE[0][l]*scale;
3437 for(Int_t l = 0; l < nbins; l++){
3438 if(fCoefCharge[4][l] > 0){
3439 sum += fCoefCharge[4][l];
3440 xValuesFittedMeanW[counter[3]]= l;
3445 if(sum > 0.0) scale = counter[3]/sum;
3446 if(fFitChargeNDB == 4){
3447 if(scale != fScaleFitFactor){
3448 AliInfo(Form("The normalisation is different from a nomalisation to one."));
3449 AliInfo(Form("For one we have %f and here %f",scale,fScaleFitFactor));
3451 AliInfo(Form("It is not normal because we didn't choose a reference database!"));
3454 scale = fScaleFitFactor;
3456 for(Int_t l = 0; l < nbins; l++){
3457 if(fCoefCharge[4][l] > 0){
3458 fCoefCharge[4][l]=fCoefCharge[4][l]*scale;
3459 fCoefChargeE[3][l]=fCoefChargeE[3][l]*scale;
3466 for(Int_t l = 0; l < nbins; l++){
3467 if(fCoefCharge[1][l] > 0){
3468 sum += fCoefCharge[1][l];
3469 xValuesFittedMean[counter[1]]= l;
3474 if(sum > 0.0) scale = counter[1]/sum;
3475 if(fFitChargeNDB == 1){
3476 if(scale != fScaleFitFactor){
3477 AliInfo(Form("The normalisation is different from a nomalisation to one."));
3478 AliInfo(Form("For one we have %f and here %f",scale,fScaleFitFactor));
3480 AliInfo(Form("It is not normal because we didn't choose a reference database!"));
3483 scale = fScaleFitFactor;
3485 for(Int_t l = 0; l < nbins; l++){
3486 if(fCoefCharge[1][l] > 0){
3487 fCoefCharge[1][l]=fCoefCharge[1][l]*scale;
3488 fCoefChargeE[1][l]=fCoefChargeE[1][l]*scale;
3493 if(fFitChargeBisOn){
3495 for(Int_t l = 0; l < nbins; l++){
3496 if(fCoefCharge[2][l] > 0){
3497 sum += fCoefCharge[2][l];
3498 xValuesFittedBis[counter[2]]= l;
3503 if(sum > 0.0) scale = counter[2]/sum;
3504 if(fFitChargeNDB == 0){
3505 if(scale != fScaleFitFactor){
3506 AliInfo(Form("The normalisation is different from a nomalisation to one."));
3507 AliInfo(Form("For one we have %f and here %f",scale,fScaleFitFactor));
3509 AliInfo(Form("It is not normal because we didn't choose a reference database!"));
3512 scale = fScaleFitFactor;
3514 for(Int_t l = 0; l < nbins; l++){
3515 if(fCoefCharge[2][l] > 0){
3516 fCoefCharge[2][l]=fCoefCharge[2][l]*scale;
3517 fCoefChargeE[2][l]=fCoefChargeE[2][l]*scale;
3522 // Create the X and Xerror
3523 Double_t *xValues = new Double_t[nbins];
3524 Double_t *xValuesE = new Double_t[nbins];
3525 for(Int_t k = 0; k < nbins; k ++){
3530 // Create the graph erros and plot them
3531 TCanvas *cch1 = new TCanvas("cch1","",50,50,600,800);
3533 TLegend *legch1 = new TLegend(0.4,0.6,0.89,0.89);
3535 TGraph *graphCharge3 = new TGraph(nbins,xValues,fCoefCharge[3]);
3536 graphCharge3->SetName("coefcharge3");
3537 graphCharge3->SetTitle("");
3538 graphCharge3->GetXaxis()->SetTitle("Det/Pad groups");
3539 graphCharge3->GetYaxis()->SetTitle("gain factor");
3540 graphCharge3->SetLineColor(4);
3541 graphCharge3->SetMarkerStyle(25);
3542 graphCharge3->SetMarkerColor(4);
3543 listofgraphs->Add((TObject *)graphCharge3);
3544 legch1->AddEntry(graphCharge3,"f_{g} simulated","p");
3545 graphCharge3->Draw("AP");
3548 TGraphErrors *graphCharge0 = new TGraphErrors(nbins,xValues,fCoefCharge[0],xValuesE,fCoefChargeE[0]);
3549 graphCharge0->SetName("coefcharge0");
3550 graphCharge0->SetTitle("");
3551 graphCharge0->GetXaxis()->SetTitle("Det/Pad groups");
3552 graphCharge0->GetYaxis()->SetTitle("gain factor");
3553 graphCharge0->SetMarkerColor(6);
3554 graphCharge0->SetLineColor(6);
3555 graphCharge0->SetMarkerStyle(26);
3556 listofgraphs->Add((TObject *)graphCharge0);
3557 legch1->AddEntry(graphCharge0,"f_{g} fit","p");
3558 graphCharge0->Draw("P");
3561 TGraphErrors *graphCharge4 = new TGraphErrors(nbins,xValues,fCoefCharge[4],xValuesE,fCoefChargeE[3]);
3562 graphCharge4->SetName("coefcharge4");
3563 graphCharge4->SetTitle("");
3564 graphCharge4->GetXaxis()->SetTitle("Det/Pad groups");
3565 graphCharge4->GetYaxis()->SetTitle("gain factor");
3566 graphCharge4->SetMarkerColor(1);
3567 graphCharge4->SetLineColor(1);
3568 graphCharge4->SetMarkerStyle(30);
3569 listofgraphs->Add((TObject *)graphCharge4);
3570 legch1->AddEntry(graphCharge4,"f_{g} Mean W","p");
3571 graphCharge4->Draw("P");
3573 if (fMeanChargeOn) {
3574 TGraphErrors *graphCharge1 = new TGraphErrors(nbins,xValues,fCoefCharge[1],xValuesE,fCoefChargeE[1]);
3575 graphCharge1->SetName("coefcharge1");
3576 graphCharge1->GetXaxis()->SetTitle("Det/Pad groups");
3577 graphCharge1->GetYaxis()->SetTitle("gain factor");
3578 graphCharge1->SetTitle("");
3579 graphCharge1->SetMarkerColor(2);
3580 graphCharge1->SetLineColor(2);
3581 graphCharge1->SetMarkerStyle(24);
3582 legch1->AddEntry(graphCharge1,"f_{g} mean","p");
3583 graphCharge1->Draw("P");
3584 listofgraphs->Add((TObject *)graphCharge1);
3586 if (fFitChargeBisOn ) {
3587 TGraphErrors *graphCharge2 = new TGraphErrors(nbins,xValues,fCoefCharge[2],xValuesE,fCoefChargeE[2]);
3588 graphCharge2->SetName("coefcharge2");
3589 graphCharge2->SetTitle("");
3590 graphCharge2->GetXaxis()->SetTitle("Det/Pad groups");
3591 graphCharge2->GetYaxis()->SetTitle("gain factor");
3592 graphCharge2->SetMarkerColor(8);
3593 graphCharge2->SetLineColor(8);
3594 graphCharge2->SetMarkerStyle(25);
3595 legch1->AddEntry(graphCharge2,"f_{g} fitbis","p");
3596 graphCharge2->Draw("P");
3597 listofgraphs->Add((TObject *)graphCharge2);
3599 legch1->Draw("same");
3601 //Create the arrays and the graphs for the delta
3603 TCanvas *cch2 = new TCanvas("cch2","",50,50,600,800);
3605 TLegend *legch3 = new TLegend(0.4,0.6,0.89,0.89);
3606 TLegend *legch2 = new TLegend(0.4,0.6,0.89,0.89);
3610 Double_t *yValuesDelta = new Double_t[counter[0]];
3611 for(Int_t k = 0; k < counter[0]; k++){
3612 if (fCoefCharge[3][(Int_t)(xValuesFitted[k])] > 0.0) {
3613 yValuesDelta[k] = (fCoefCharge[0][(Int_t)xValuesFitted[k]]-fCoefCharge[3][(Int_t)xValuesFitted[k]])
3614 / fCoefCharge[3][(Int_t)xValuesFitted[k]];
3617 yValuesDelta[k] = 0.0;
3620 TGraph *graphDeltaCharge0 = new TGraph(counter[0],&xValuesFitted[0],yValuesDelta);
3621 graphDeltaCharge0->SetName("deltacharge0");
3622 graphDeltaCharge0->GetXaxis()->SetTitle("Det/Pad groups");
3623 graphDeltaCharge0->GetYaxis()->SetTitle("#Deltag/g_{sim}");
3624 graphDeltaCharge0->SetMarkerColor(6);
3625 graphDeltaCharge0->SetTitle("");
3626 graphDeltaCharge0->SetLineColor(6);
3627 graphDeltaCharge0->SetMarkerStyle(26);
3628 listofgraphs->Add((TObject *)graphDeltaCharge0);
3629 legch3->AddEntry(graphDeltaCharge0,"fit","p");
3630 graphDeltaCharge0->Draw("AP");
3633 TH1I *histoErrorCharge0 = new TH1I("errorcharge0","",100 ,-0.10,0.10);
3634 histoErrorCharge0->SetXTitle("#Deltag/g_{sim}");
3635 histoErrorCharge0->SetYTitle("counts");
3636 histoErrorCharge0->SetLineColor(6);
3637 histoErrorCharge0->SetLineStyle(1);
3638 histoErrorCharge0->SetStats(0);
3639 Double_t maxvalue = 0.0;
3640 for(Int_t k = 0; k < counter[0]; k++){
3641 histoErrorCharge0->Fill(yValuesDelta[k]);
3642 if(k == 0) maxvalue = TMath::Abs(yValuesDelta[k]);
3643 if(maxvalue < (TMath::Abs(yValuesDelta[k]))) maxvalue = TMath::Abs(yValuesDelta[k]);
3645 AliInfo(Form("The maximum deviation found dor the fit method is %f",maxvalue));
3646 legch2->AddEntry(histoErrorCharge0,"f_{g} fit","l");
3647 histoErrorCharge0->Draw();
3648 listofgraphs->Add((TObject *)histoErrorCharge0);
3654 Double_t *yValuesDelta = new Double_t[counter[3]];
3655 for(Int_t k = 0; k < counter[3]; k++){
3656 if (fCoefCharge[3][(Int_t)(xValuesFittedMeanW[k])] > 0.0) {
3657 yValuesDelta[k] = (fCoefCharge[4][(Int_t)xValuesFittedMeanW[k]]-fCoefCharge[3][(Int_t)xValuesFittedMeanW[k]])
3658 / fCoefCharge[3][(Int_t)xValuesFittedMeanW[k]];
3661 yValuesDelta[k] = 0.0;
3664 TGraph *graphDeltaCharge4 = new TGraph(counter[3],&xValuesFittedMeanW[0],yValuesDelta);
3665 graphDeltaCharge4->SetName("deltacharge4");
3666 graphDeltaCharge4->GetXaxis()->SetTitle("Det/Pad groups");
3667 graphDeltaCharge4->GetYaxis()->SetTitle("#Deltag/g_{sim}");
3668 graphDeltaCharge4->SetMarkerColor(1);
3669 graphDeltaCharge4->SetTitle("");
3670 graphDeltaCharge4->SetLineColor(1);
3671 graphDeltaCharge4->SetMarkerStyle(30);
3672 listofgraphs->Add((TObject *)graphDeltaCharge4);
3673 legch3->AddEntry(graphDeltaCharge4,"Mean W","p");
3675 graphDeltaCharge4->Draw("AP");
3678 graphDeltaCharge4->Draw("P");
3682 TH1I *histoErrorCharge4 = new TH1I("errorcharge4","",100 ,-0.10,0.10);
3683 histoErrorCharge4->SetXTitle("#Deltag/g_{sim}");
3684 histoErrorCharge4->SetYTitle("counts");
3685 histoErrorCharge4->SetLineColor(1);
3686 histoErrorCharge4->SetLineStyle(1);
3687 histoErrorCharge4->SetStats(0);
3688 Double_t maxvalue = 0.0;
3689 for(Int_t k = 0; k < counter[3]; k++){
3690 histoErrorCharge4->Fill(yValuesDelta[k]);
3691 if(k == 0) maxvalue = yValuesDelta[k];
3692 if(maxvalue < (TMath::Abs(yValuesDelta[k]))) maxvalue = TMath::Abs(yValuesDelta[k]);
3694 AliInfo(Form("The maximum deviation found for the meanW method is %f",maxvalue));
3695 legch2->AddEntry(histoErrorCharge4,"f_{g} Mean W","l");
3697 histoErrorCharge4->Draw();
3700 histoErrorCharge4->Draw("same");
3702 listofgraphs->Add((TObject *)histoErrorCharge4);
3706 if (fMeanChargeOn) {
3708 Double_t *yValuesDeltaMean = new Double_t[counter[1]];
3709 for (Int_t k = 0; k < counter[1]; k++){
3710 if (fCoefCharge[3][(Int_t)xValuesFittedMean[k]] > 0.0) {
3711 yValuesDeltaMean[k] = (fCoefCharge[1][(Int_t)xValuesFittedMean[k]]-fCoefCharge[3][(Int_t)xValuesFittedMean[k]])
3712 / fCoefCharge[3][(Int_t)xValuesFittedMean[k]];
3715 yValuesDeltaMean[k] = 0.0;
3718 TGraph *graphDeltaCharge1 = new TGraph(counter[1],&xValuesFittedMean[0],yValuesDeltaMean);
3719 graphDeltaCharge1->SetName("deltacharge1");
3720 graphDeltaCharge1->GetXaxis()->SetTitle("Det/Pad groups");
3721 graphDeltaCharge1->GetYaxis()->SetTitle("#Deltag/g_{sim}");
3722 graphDeltaCharge1->SetMarkerColor(2);
3723 graphDeltaCharge1->SetMarkerStyle(24);
3724 graphDeltaCharge1->SetLineColor(2);
3725 graphDeltaCharge1->SetTitle("");
3726 legch3->AddEntry(graphDeltaCharge1,"mean","p");
3728 graphDeltaCharge1->Draw("AP");
3731 graphDeltaCharge1->Draw("P");
3733 listofgraphs->Add((TObject *)graphDeltaCharge1);
3736 TH1I *histoErrorCharge1 = new TH1I("errorcharge1","",100 ,-0.10,0.10);
3737 histoErrorCharge1->SetXTitle("#Deltag/g_{sim}");
3738 histoErrorCharge1->SetYTitle("counts");
3739 histoErrorCharge1->SetLineColor(2);
3740 histoErrorCharge1->SetLineStyle(2);
3741 histoErrorCharge1->SetStats(0);
3742 Double_t maxvalue = 0.0;
3743 for(Int_t k = 0; k < counter[1]; k++){
3744 histoErrorCharge1->Fill(yValuesDeltaMean[k]);
3745 if(k == 0) maxvalue = TMath::Abs(yValuesDeltaMean[k]);
3746 if(maxvalue < (TMath::Abs(yValuesDeltaMean[k]))) maxvalue = TMath::Abs(yValuesDeltaMean[k]);
3748 AliInfo(Form("The maximum deviation found for the mean method is %f",maxvalue));
3749 legch2->AddEntry(histoErrorCharge1,"f_{g} mean","l");
3751 histoErrorCharge1->Draw();
3754 histoErrorCharge1->Draw("same");
3756 listofgraphs->Add((TObject *)histoErrorCharge1);
3760 if (fFitChargeBisOn) {
3762 Double_t *yValuesDeltaBis = new Double_t[counter[2]];
3763 for(Int_t k = 0; k < counter[2]; k++){
3764 if (fCoefCharge[3][(Int_t)xValuesFittedBis[k]] > 0.0) {
3765 yValuesDeltaBis[k] = (fCoefCharge[2][(Int_t)xValuesFittedBis[k]]-fCoefCharge[3][(Int_t)xValuesFittedBis[k]])
3766 / fCoefCharge[3][(Int_t)xValuesFittedBis[k]];
3769 yValuesDeltaBis[k] = 0.0;
3772 TGraph *graphDeltaCharge2 = new TGraph(counter[2],&xValuesFittedBis[0],yValuesDeltaBis);
3773 graphDeltaCharge2->SetName("deltacharge2");
3774 graphDeltaCharge2->GetXaxis()->SetTitle("Det/Pad groups");
3775 graphDeltaCharge2->GetYaxis()->SetTitle("#Deltag/g_{sim}");
3776 graphDeltaCharge2->SetMarkerColor(8);
3777 graphDeltaCharge2->SetLineColor(8);
3778 graphDeltaCharge2->SetMarkerStyle(25);
3779 legch3->AddEntry(graphDeltaCharge2,"fit","p");
3780 graphDeltaCharge2->SetTitle("");
3782 graphDeltaCharge2->Draw("AP");
3785 graphDeltaCharge2->Draw("P");
3787 listofgraphs->Add((TObject *)graphDeltaCharge2);
3790 TH1I *histoErrorCharge2 = new TH1I("errorcharge2","",100 ,-0.10, 0.10);
3791 histoErrorCharge2->SetXTitle("#Deltag/g_{sim}");
3792 histoErrorCharge2->SetYTitle("counts");
3793 histoErrorCharge2->SetLineColor(8);
3794 histoErrorCharge2->SetLineStyle(5);
3795 histoErrorCharge2->SetLineWidth(3);
3796 histoErrorCharge2->SetStats(0);
3797 Double_t maxvalue = 0.0;
3798 for(Int_t k = 0; k < counter[2]; k++){
3799 histoErrorCharge2->Fill(yValuesDeltaBis[k]);
3800 if(k == 0) maxvalue = TMath::Abs(yValuesDeltaBis[k]);
3801 if(maxvalue < (TMath::Abs(yValuesDeltaBis[k]))) maxvalue = TMath::Abs(yValuesDeltaBis[k]);
3803 AliInfo(Form("The maximum deviation found for the fit bis method is %f",maxvalue));
3804 legch2->AddEntry(histoErrorCharge2,"f_{g} fitbis","l");
3806 histoErrorCharge2->Draw();
3809 histoErrorCharge2->Draw("same");
3811 listofgraphs->Add((TObject *)histoErrorCharge2);
3812 //it doesn't matter anymore but...
3817 legch3->Draw("same");
3819 legch2->Draw("same");
3823 TFile *fout = TFile::Open(fWriteNameCoef,"UPDATE");
3824 // Check if the file could be opened
3825 if (!fout || !fout->IsOpen()) {
3826 AliInfo("No File found!");
3830 for(Int_t k = 0; k < listofgraphs->GetEntriesFast(); k++){
3831 fout->WriteTObject((TObject *) listofgraphs->At(k),((TObject *)listofgraphs->At(k))->GetName()
3832 ,(Option_t *) "OverWrite");
3840 //_____________________________________________________________________________
3841 void AliTRDCalibraFit::PlotWritePH()
3844 // create the graph errors and write them if wanted
3847 //TObjArray of the grapherrors and so on
3848 TObjArray *listofgraphs = new TObjArray();
3850 Int_t nbins = fDect2[1]-fDect1[1];
3852 //See the number of fitted for delta
3860 Double_t *xValuesFitted = new Double_t[nbins];
3861 Double_t *xValuesFittedPH = new Double_t[nbins];
3862 Double_t *xValuesFittedLP = new Double_t[nbins];
3863 for(Int_t k = 0; k < nbins; k ++){
3864 xValuesFitted[k] = -1;
3865 xValuesFittedPH[k] = -1;
3866 xValuesFittedLP[k] = -1;
3870 for(Int_t l = 0; l < nbins; l++){
3871 if(fCoefVdrift[1][l] > 0){
3872 xValuesFitted[counter[1]]=l;
3878 for(Int_t l = 0; l < nbins; l++){
3879 if(fCoefVdrift[3][l] > 0){
3880 xValuesFittedLP[counter[2]]=l;
3886 for(Int_t l = 0; l < nbins; l++){
3887 if(fCoefVdrift[0][l] > 0){
3888 xValuesFittedPH[counter[0]]= l;
3894 //Create the X and Xerror
3895 Double_t *xValues = new Double_t[nbins];
3896 Double_t *xValuesE = new Double_t[nbins];
3897 for(Int_t k = 0; k < nbins; k ++){
3902 //Create the graph erros and plot them
3903 TCanvas *cph1 = new TCanvas("cph1","",50,50,600,800);
3906 TGraph *graphVdrift2 = new TGraph(nbins,xValues,fCoefVdrift[2]);
3907 graphVdrift2->SetName("coefvdrift2");
3908 graphVdrift2->SetTitle("");
3909 graphVdrift2->GetXaxis()->SetTitle("Det/Pad groups");
3910 graphVdrift2->GetYaxis()->SetTitle("Vdrift [cm/#mus]");
3911 graphVdrift2->SetLineColor(4);
3912 listofgraphs->Add((TObject *)graphVdrift2);
3913 TLegend *legph1 = new TLegend(0.4,0.6,0.89,0.89);
3914 legph1->AddEntry(graphVdrift2,"Vdrift simulated","l");
3915 graphVdrift2->Draw("AL");
3918 TGraphErrors *graphVdrift1 = new TGraphErrors(nbins,xValues,fCoefVdrift[1],xValuesE,fCoefVdriftE[1]);
3919 graphVdrift1->SetName("coefvdrift1");
3920 graphVdrift1->SetTitle("");
3921 graphVdrift1->GetXaxis()->SetTitle("Det/Pad groups");
3922 graphVdrift1->GetYaxis()->SetTitle("Vdrift [cm/#mus]");
3923 graphVdrift1->SetMarkerColor(6);
3924 graphVdrift1->SetLineColor(6);
3925 graphVdrift1->SetMarkerStyle(26);
3926 listofgraphs->Add((TObject *)graphVdrift1);
3927 legph1->AddEntry(graphVdrift1,"Vdrift fit","p");
3928 graphVdrift1->Draw("P");
3931 TGraphErrors *graphVdrift0 = new TGraphErrors(nbins,xValues,fCoefVdrift[0],xValuesE,fCoefVdriftE[0]);
3932 graphVdrift0->SetName("coefVdrift0");
3933 graphVdrift0->GetXaxis()->SetTitle("Det/Pad groups");
3934 graphVdrift0->GetYaxis()->SetTitle("Vdrift [cm/#mus]");
3935 graphVdrift0->SetTitle("");
3936 graphVdrift0->SetMarkerColor(2);
3937 graphVdrift0->SetLineColor(2);
3938 graphVdrift0->SetMarkerStyle(24);
3939 legph1->AddEntry(graphVdrift0,"v_{fit PH}","p");
3940 graphVdrift0->Draw("P");
3941 listofgraphs->Add((TObject *)graphVdrift0);
3943 if (fFitLagrPolOn) {
3944 TGraphErrors *graphVdrift3 = new TGraphErrors(nbins,xValues,fCoefVdrift[3],xValuesE,fCoefVdriftE[2]);
3945 graphVdrift3->SetName("coefVdrift3");
3946 graphVdrift3->GetXaxis()->SetTitle("Det/Pad groups");
3947 graphVdrift3->GetYaxis()->SetTitle("Vdrift [cm/#mus]");
3948 graphVdrift3->SetTitle("");
3949 graphVdrift3->SetMarkerColor(1);
3950 graphVdrift3->SetLineColor(1);
3951 graphVdrift3->SetMarkerStyle(28);
3952 legph1->AddEntry(graphVdrift3,"v_{LagrPol}","p");
3953 graphVdrift3->Draw("P");
3954 listofgraphs->Add((TObject *)graphVdrift3);
3956 legph1->Draw("same");
3958 //Create the arrays and the graphs for the delta
3959 TCanvas *cph2 = new TCanvas("cph2","",50,50,600,800);
3961 TLegend *legph3 = new TLegend(0.4,0.6,0.89,0.89);
3962 TLegend *legph2 = new TLegend(0.4,0.6,0.89,0.89);
3967 Double_t *yValuesDelta = new Double_t[counter[1]];
3968 for (Int_t k = 0; k < counter[1]; k++){
3969 if (fCoefVdrift[2][(Int_t)(xValuesFitted[k])] > 0.0) {
3970 yValuesDelta[k] = (fCoefVdrift[1][(Int_t)xValuesFitted[k]]-fCoefVdrift[2][(Int_t)xValuesFitted[k]])
3971 / fCoefVdrift[2][(Int_t)xValuesFitted[k]];
3974 yValuesDelta[k] = 0.0;
3977 TGraph *graphDeltaVdrift1 = new TGraph(counter[1],&xValuesFitted[0],yValuesDelta);
3978 graphDeltaVdrift1->SetName("deltavdrift1");
3979 graphDeltaVdrift1->GetXaxis()->SetTitle("Det/Pad groups");
3980 graphDeltaVdrift1->GetYaxis()->SetTitle("#Deltav/v_{sim}");
3981 graphDeltaVdrift1->SetMarkerColor(6);
3982 graphDeltaVdrift1->SetTitle("");
3983 graphDeltaVdrift1->SetLineColor(6);
3984 graphDeltaVdrift1->SetMarkerStyle(26);
3985 listofgraphs->Add((TObject *)graphDeltaVdrift1);
3986 legph3->AddEntry(graphDeltaVdrift1,"v_{slope method}","p");
3987 graphDeltaVdrift1->Draw("AP");
3990 TH1I *histoErrorVdrift1 = new TH1I("errorvdrift1","",100 ,-0.2,0.2);
3991 histoErrorVdrift1->SetXTitle("#Deltav/v_{sim}");
3992 histoErrorVdrift1->SetYTitle("counts");
3993 histoErrorVdrift1->SetLineColor(6);
3994 histoErrorVdrift1->SetLineStyle(1);
3995 histoErrorVdrift1->SetStats(0);
3996 Double_t maxvalue = 0.0;
3997 for(Int_t k = 0; k < counter[1]; k++){
3998 histoErrorVdrift1->Fill(yValuesDelta[k]);
3999 if(k == 0) maxvalue = yValuesDelta[k];
4000 if(maxvalue < (TMath::Abs(yValuesDelta[k]))) maxvalue = TMath::Abs(yValuesDelta[k]);
4002 AliInfo(Form("The maximum deviation found for the Pol2 method is %f",maxvalue));
4003 legph2->AddEntry(histoErrorVdrift1,"v_{slope method}","l");
4004 histoErrorVdrift1->Draw();
4005 listofgraphs->Add((TObject *)histoErrorVdrift1);
4011 Double_t *yValuesDeltaPH = new Double_t[counter[0]];
4012 for(Int_t k = 0; k < counter[0]; k++){
4013 if(fCoefVdrift[2][(Int_t)xValuesFittedPH[k]] > 0.0) {
4014 yValuesDeltaPH[k] = (fCoefVdrift[0][(Int_t)xValuesFittedPH[k]]-fCoefVdrift[2][(Int_t)xValuesFittedPH[k]])/fCoefVdrift[2][(Int_t)xValuesFittedPH[k]];
4016 else yValuesDeltaPH[k] = 0.0;
4018 TGraph *graphDeltaVdrift0 = new TGraph(counter[0],&xValuesFittedPH[0],yValuesDeltaPH);
4019 graphDeltaVdrift0->SetName("deltavdrift0");
4020 graphDeltaVdrift0->GetXaxis()->SetTitle("Det/Pad groups");
4021 graphDeltaVdrift0->GetYaxis()->SetTitle("#Deltav/v_{sim}");
4022 graphDeltaVdrift0->SetMarkerColor(2);
4023 graphDeltaVdrift0->SetMarkerStyle(24);
4024 graphDeltaVdrift0->SetLineColor(2);
4025 graphDeltaVdrift0->SetTitle("");
4026 legph3->AddEntry(graphDeltaVdrift0,"v_{fit PH}","p");
4028 graphDeltaVdrift0->Draw("P");
4031 graphDeltaVdrift0->Draw("AP");
4033 listofgraphs->Add((TObject *)graphDeltaVdrift0);
4035 TH1I *histoErrorVdrift0 = new TH1I("errorvdrift0","",100 ,-0.2,0.2);
4036 histoErrorVdrift0->SetXTitle("#Deltav/v_{sim}");
4037 histoErrorVdrift0->SetYTitle("counts");
4038 histoErrorVdrift0->SetLineColor(2);
4039 histoErrorVdrift0->SetLineStyle(2);
4040 histoErrorVdrift0->SetStats(0);
4041 Double_t maxvalue = 0.0;
4042 for(Int_t k = 0; k < counter[0]; k++){
4043 histoErrorVdrift0->Fill(yValuesDeltaPH[k]);
4044 if(k == 0) maxvalue = yValuesDeltaPH[k];
4045 if(maxvalue < (TMath::Abs(yValuesDeltaPH[k]))) maxvalue = TMath::Abs(yValuesDeltaPH[k]);
4047 AliInfo(Form("The maximum deviation found for the fit method is %f",maxvalue));
4048 legph2->AddEntry(histoErrorVdrift0,"v_{fit PH}","l");
4050 histoErrorVdrift0->Draw("same");
4053 histoErrorVdrift0->Draw();
4055 listofgraphs->Add((TObject *)histoErrorVdrift0);
4059 if (fFitLagrPolOn) {
4061 Double_t *yValuesDeltaPH = new Double_t[counter[2]];
4062 for (Int_t k = 0; k < counter[2]; k++){
4063 if (fCoefVdrift[2][(Int_t)xValuesFittedLP[k]] > 0.0) {
4064 yValuesDeltaPH[k] = (fCoefVdrift[3][(Int_t)xValuesFittedLP[k]]-fCoefVdrift[2][(Int_t)xValuesFittedLP[k]])
4065 / fCoefVdrift[2][(Int_t)xValuesFittedLP[k]];
4068 yValuesDeltaPH[k] = 0.0;
4071 TGraph *graphDeltaVdrift3 = new TGraph(counter[2],&xValuesFittedLP[0],yValuesDeltaPH);
4072 graphDeltaVdrift3->SetName("deltavdrift3");
4073 graphDeltaVdrift3->GetXaxis()->SetTitle("Det/Pad groups");
4074 graphDeltaVdrift3->GetYaxis()->SetTitle("#Deltav/v_{sim}");
4075 graphDeltaVdrift3->SetMarkerColor(1);
4076 graphDeltaVdrift3->SetMarkerStyle(28);
4077 graphDeltaVdrift3->SetLineColor(1);
4078 graphDeltaVdrift3->SetTitle("");
4079 legph3->AddEntry(graphDeltaVdrift3,"v_{LagrPol}","p");
4081 graphDeltaVdrift3->Draw("P");
4084 graphDeltaVdrift3->Draw("AP");
4086 listofgraphs->Add((TObject *)graphDeltaVdrift3);
4088 TH1I *histoErrorVdrift3 = new TH1I("errorvdrift3","",100 ,-0.2,0.2);
4089 histoErrorVdrift3->SetXTitle("#Deltav/v_{sim}");
4090 histoErrorVdrift3->SetYTitle("counts");
4091 histoErrorVdrift3->SetLineColor(1);
4092 histoErrorVdrift3->SetLineStyle(1);
4093 histoErrorVdrift3->SetStats(0);
4094 Double_t maxvalue = 0.0;
4095 for(Int_t k = 0; k < counter[2]; k++){
4096 histoErrorVdrift3->Fill(yValuesDeltaPH[k]);
4097 if(k == 0) maxvalue = yValuesDeltaPH[k];
4098 if(maxvalue < (TMath::Abs(yValuesDeltaPH[k]))) maxvalue = TMath::Abs(yValuesDeltaPH[k]);
4100 AliInfo(Form("The maximum deviation found for the LagrPol method is %f",maxvalue));
4101 legph2->AddEntry(histoErrorVdrift3,"v_{LagrPol}","l");
4103 histoErrorVdrift3->Draw("same");
4106 histoErrorVdrift3->Draw();
4108 listofgraphs->Add((TObject *)histoErrorVdrift3);
4112 legph3->Draw("same");
4114 legph2->Draw("same");
4118 TFile *fout = TFile::Open(fWriteNameCoef,"UPDATE");
4119 // Check if the file could be opened
4120 if (!fout || !fout->IsOpen()) {
4121 AliInfo("No File found!");
4125 for(Int_t k = 0; k < listofgraphs->GetEntriesFast(); k++){
4126 fout->WriteTObject((TObject *) listofgraphs->At(k),((TObject *)listofgraphs->At(k))->GetName()
4127 ,(Option_t *) "OverWrite");
4135 //_____________________________________________________________________________
4136 void AliTRDCalibraFit::PlotWriteT0()
4139 // create the graph errors and write them if wanted
4142 //TObjArray of the grapherrors and so on
4143 TObjArray *listofgraphs = new TObjArray();
4145 Int_t nbins = fDect2[1]-fDect1[1];
4147 //See the number of fitted for delta: here T0 can be negative, we don't use the sign but the error
4148 //and the grapherrors of the coefficients contained the no fitted with error 0.0
4156 Double_t *xValuesFitted = new Double_t[nbins];
4157 Double_t *xValuesFittedPH = new Double_t[nbins];
4158 Double_t *xValuesFittedLP = new Double_t[nbins];
4159 for(Int_t k = 0; k < nbins; k ++){
4160 xValuesFitted[k] = -1;
4161 xValuesFittedPH[k] = -1;
4162 xValuesFittedLP[k] = -1;
4166 for(Int_t l = 0; l < nbins; l++){
4167 if(fCoefT0E[1][l] != 0.0){
4168 xValuesFitted[counter[1]]=l;
4175 for(Int_t l = 0; l < nbins; l++){
4176 if(fCoefT0E[0][l] != 0.0){
4177 xValuesFittedPH[counter[0]]= l;
4184 for(Int_t l = 0; l < nbins; l++){
4185 if(fCoefT0E[2][l] == 1.0){
4186 xValuesFittedLP[counter[2]]= l;
4192 //Create the X and Xerror
4193 Double_t *xValues = new Double_t[nbins];
4194 Double_t *xValuesE = new Double_t[nbins];
4195 for(Int_t k = 0; k < nbins; k ++){
4200 //Create the graph erros and plot them
4201 TCanvas *ct01 = new TCanvas("ct01","",50,50,600,800);
4203 TLegend *legt01 = new TLegend(0.4,0.6,0.89,0.89);
4205 TGraph *graphT02 = new TGraph(nbins,xValues,fCoefT0[2]);
4206 graphT02->SetName("coeft02");
4207 graphT02->SetTitle("");
4208 graphT02->GetXaxis()->SetTitle("Det/Pad groups");
4209 graphT02->GetYaxis()->SetTitle("T0 [time bins]");
4210 graphT02->SetLineColor(4);
4211 listofgraphs->Add((TObject *)graphT02);
4212 legt01->AddEntry(graphT02,"T0 simulated","l");
4213 graphT02->Draw("AL");
4216 TGraphErrors *graphT01 = new TGraphErrors(nbins,xValues,fCoefT0[1],xValuesE,fCoefT0E[1]);
4217 graphT01->SetName("coeft01");
4218 graphT01->SetTitle("");
4219 graphT01->GetXaxis()->SetTitle("Det/Pad groups");
4220 graphT01->GetYaxis()->SetTitle("T0 [time bins]");
4221 graphT01->SetMarkerColor(6);
4222 graphT01->SetLineColor(6);
4223 graphT01->SetMarkerStyle(26);
4224 listofgraphs->Add((TObject *)graphT01);
4225 legt01->AddEntry(graphT01,"T0 slope method","p");
4226 graphT01->Draw("P");
4229 TGraphErrors *graphT00 = new TGraphErrors(nbins,xValues,fCoefT0[0],xValuesE,fCoefT0E[0]);
4230 graphT00->SetName("coeft00");
4231 graphT00->GetXaxis()->SetTitle("Det/Pad groups");
4232 graphT00->GetYaxis()->SetTitle("T0 [time bins]");
4233 graphT00->SetTitle("");
4234 graphT00->SetMarkerColor(2);
4235 graphT00->SetLineColor(2);
4236 graphT00->SetMarkerStyle(24);
4237 legt01->AddEntry(graphT00,"T0 fit","p");
4238 graphT00->Draw("P");
4239 listofgraphs->Add((TObject *)graphT00);
4241 if (fFitLagrPolOn) {
4242 TGraphErrors *graphT03 = new TGraphErrors(nbins,xValues,fCoefT0[3],xValuesE,xValuesE);
4243 graphT03->SetName("coeft03");
4244 graphT03->GetXaxis()->SetTitle("Det/Pad groups");
4245 graphT03->GetYaxis()->SetTitle("T0 [time bins]");
4246 graphT03->SetTitle("");
4247 graphT03->SetMarkerColor(1);
4248 graphT03->SetLineColor(1);
4249 graphT03->SetMarkerStyle(28);
4250 legt01->AddEntry(graphT03,"T0 LagrPol","p");
4251 graphT03->Draw("P");
4252 listofgraphs->Add((TObject *)graphT03);
4254 legt01->Draw("same");
4256 //Create the arrays and the graphs for the delta
4257 TCanvas *ct02 = new TCanvas("ct02","",50,50,600,800);
4259 TLegend *legt03 = new TLegend(0.4,0.6,0.89,0.89);
4260 TLegend *legt02 = new TLegend(0.4,0.6,0.89,0.89);
4265 Double_t *yValuesDelta = new Double_t[counter[1]];
4266 for(Int_t k = 0; k < counter[1]; k++){
4267 yValuesDelta[k] = (fCoefT0[1][(Int_t)xValuesFitted[k]]-fCoefT0[2][(Int_t)xValuesFitted[k]]);
4269 TGraph *graphDeltaT01 = new TGraph(counter[1],&xValuesFitted[0],yValuesDelta);
4270 graphDeltaT01->SetName("deltat01");
4271 graphDeltaT01->GetXaxis()->SetTitle("Det/Pad groups");
4272 graphDeltaT01->GetYaxis()->SetTitle("#Deltat0 [time bins]");
4273 graphDeltaT01->SetMarkerColor(6);
4274 graphDeltaT01->SetTitle("");
4275 graphDeltaT01->SetLineColor(6);
4276 graphDeltaT01->SetMarkerStyle(26);
4277 listofgraphs->Add((TObject *)graphDeltaT01);
4278 legt03->AddEntry(graphDeltaT01,"T0_{slope method}","p");
4279 graphDeltaT01->Draw("AP");
4282 TH1I *histoErrorT01 = new TH1I("errort01","",100 ,-0.2,0.2);
4283 histoErrorT01->SetXTitle("#Deltat0 [time bins]");
4284 histoErrorT01->SetYTitle("counts");
4285 histoErrorT01->SetLineColor(6);
4286 histoErrorT01->SetLineStyle(1);
4287 histoErrorT01->SetStats(0);
4288 Double_t maxvalue = 0.0;
4289 for(Int_t k = 0; k < counter[1]; k++){
4290 histoErrorT01->Fill(yValuesDelta[k]);
4291 if(k == 0) maxvalue = yValuesDelta[k];
4292 if(maxvalue < (TMath::Abs(yValuesDelta[k]))) maxvalue = (TMath::Abs(yValuesDelta[k]));
4294 AliInfo(Form("The maximum deviation found for the Pol2 method is %f",maxvalue));
4295 legt02->AddEntry(histoErrorT01,"T0_{slope method}","l");
4296 histoErrorT01->Draw();
4297 listofgraphs->Add((TObject *)histoErrorT01);
4302 Double_t *yValuesDeltaPH = new Double_t[counter[0]];
4303 for(Int_t k = 0; k < counter[0]; k++){
4304 yValuesDeltaPH[k] = (fCoefT0[0][(Int_t)xValuesFittedPH[k]]-fCoefT0[2][(Int_t)xValuesFittedPH[k]]);
4306 TGraph *graphDeltaT00 = new TGraph(counter[0],&xValuesFittedPH[0],yValuesDeltaPH);
4307 graphDeltaT00->SetName("deltat00");
4308 graphDeltaT00->GetXaxis()->SetTitle("Det/Pad groups");
4309 graphDeltaT00->GetYaxis()->SetTitle("#Deltat0 [time bins]");
4310 graphDeltaT00->SetMarkerColor(2);
4311 graphDeltaT00->SetMarkerStyle(24);
4312 graphDeltaT00->SetLineColor(2);
4313 graphDeltaT00->SetTitle("");
4314 legt03->AddEntry(graphDeltaT00,"T0_{fit PH}","p");
4316 graphDeltaT00->Draw("P");
4319 graphDeltaT00->Draw("AP");
4321 listofgraphs->Add((TObject *)graphDeltaT00);
4323 TH1I *histoErrorT00 = new TH1I("errort00","",100 ,-0.2,0.2);
4324 histoErrorT00->SetXTitle("#Deltat0 [time bins]");
4325 histoErrorT00->SetYTitle("counts");
4326 histoErrorT00->SetLineColor(2);
4327 histoErrorT00->SetLineStyle(2);
4328 histoErrorT00->SetStats(0);
4329 Double_t maxvalue = 0.0;
4330 for(Int_t k = 0; k < counter[0]; k++){
4331 histoErrorT00->Fill(yValuesDeltaPH[k]);
4332 if(k == 0) maxvalue = yValuesDeltaPH[k];
4333 if(maxvalue < (TMath::Abs(yValuesDeltaPH[k]))) maxvalue = (TMath::Abs(yValuesDeltaPH[k]));
4335 AliInfo(Form("The maximum deviation found for the fit method is %f",maxvalue));
4336 legt02->AddEntry(histoErrorT00,"T0_{fit PH}","l");
4338 histoErrorT00->Draw("same");
4341 histoErrorT00->Draw();
4343 listofgraphs->Add((TObject *)histoErrorT00);
4347 if (fFitLagrPolOn) {
4349 Double_t *yValuesDeltaPH = new Double_t[counter[2]];
4350 for(Int_t k = 0; k < counter[2]; k++){
4351 yValuesDeltaPH[k] = (fCoefT0[3][(Int_t)xValuesFittedLP[k]]-fCoefT0[2][(Int_t)xValuesFittedLP[k]]);
4353 TGraph *graphDeltaT03 = new TGraph(counter[2],&xValuesFittedLP[0],yValuesDeltaPH);
4354 graphDeltaT03->SetName("deltat03");
4355 graphDeltaT03->GetXaxis()->SetTitle("Det/Pad groups");
4356 graphDeltaT03->GetYaxis()->SetTitle("#Deltat0 [time bins]");
4357 graphDeltaT03->SetMarkerColor(1);
4358 graphDeltaT03->SetMarkerStyle(28);
4359 graphDeltaT03->SetLineColor(1);
4360 graphDeltaT03->SetTitle("");
4361 legt03->AddEntry(graphDeltaT03,"T0_{LagrPol}","p");
4363 graphDeltaT03->Draw("P");
4366 graphDeltaT03->Draw("AP");
4368 listofgraphs->Add((TObject *)graphDeltaT03);
4370 TH1I *histoErrorT03 = new TH1I("errort03","",100 ,-0.2,0.2);
4371 histoErrorT03->SetXTitle("#Deltat0 [time bins]");
4372 histoErrorT03->SetYTitle("counts");
4373 histoErrorT03->SetLineColor(1);
4374 histoErrorT03->SetLineStyle(1);
4375 histoErrorT03->SetStats(0);
4376 Double_t maxvalue = 0.0;
4377 for(Int_t k = 0; k < counter[2]; k++){
4378 histoErrorT03->Fill(yValuesDeltaPH[k]);
4379 if(k == 0) maxvalue = yValuesDeltaPH[k];
4380 if(maxvalue < (TMath::Abs(yValuesDeltaPH[k]))) maxvalue = (TMath::Abs(yValuesDeltaPH[k]));
4382 AliInfo(Form("The maximum deviation found for the LagrPol method is %f",maxvalue));
4383 legt02->AddEntry(histoErrorT03,"T0_{LagrPol}","l");
4385 histoErrorT03->Draw("same");
4388 histoErrorT03->Draw();
4390 listofgraphs->Add((TObject *)histoErrorT03);
4395 legt03->Draw("same");
4397 legt02->Draw("same");
4401 TFile *fout = TFile::Open(fWriteNameCoef,"UPDATE");
4402 // Check if the file could be opened
4403 if (!fout || !fout->IsOpen()) {
4404 AliInfo("No File found!");
4408 for(Int_t k = 0; k < listofgraphs->GetEntriesFast(); k++){
4409 fout->WriteTObject((TObject *) listofgraphs->At(k),((TObject *)listofgraphs->At(k))->GetName()
4410 ,(Option_t *) "OverWrite");
4418 //_____________________________________________________________________________
4419 void AliTRDCalibraFit::PlotWritePRF()
4422 // create the graph errors and write them if wanted
4425 //TObjArray of the grapherrors and so on
4426 TObjArray *listofgraphs = new TObjArray();
4428 Int_t nbins = fDect2[2]-fDect1[2];
4430 //See the number of fitted for delta
4437 Double_t *xValuesFitted = new Double_t[nbins];
4438 for(Int_t k = 0; k < nbins; k ++){
4439 xValuesFitted[k] = -1;
4441 Double_t *xValuesRMS = new Double_t[nbins];
4442 for(Int_t k = 0; k < nbins; k ++){
4447 for(Int_t l = 0; l < nbins; l++){
4448 if(fCoefPRF[0][l] > 0){
4449 xValuesFitted[counter[0]]=l;
4455 for(Int_t l = 0; l < nbins; l++){
4456 if(fCoefPRF[2][l] > 0){
4457 xValuesRMS[counter[1]]=l;
4464 //Create the X and Xerror
4465 Double_t *xValues = new Double_t[nbins];
4466 Double_t *xValuesE = new Double_t[nbins];
4467 for(Int_t k = 0; k < nbins; k ++){
4472 //Create the graph erros and plot them
4473 TCanvas *cprf1 = new TCanvas("cprf1","",50,50,600,800);
4475 TLegend *legprf1 = new TLegend(0.4,0.6,0.89,0.89);
4477 TGraph *graphPRF1 = new TGraph(nbins,xValues,fCoefPRF[1]);
4478 graphPRF1->SetName("coefprf1");
4479 graphPRF1->SetTitle("");
4480 graphPRF1->GetXaxis()->SetTitle("Det/Pad groups");
4481 graphPRF1->GetYaxis()->SetTitle("PRF width [p.u]");
4482 graphPRF1->SetLineColor(4);
4483 graphPRF1->SetMarkerColor(4);
4484 graphPRF1->SetMarkerStyle(25);
4485 graphPRF1->SetMarkerSize(0.7);
4486 listofgraphs->Add((TObject *)graphPRF1);
4487 legprf1->AddEntry(graphPRF1,"PRF width simulated","p");
4488 graphPRF1->Draw("AP");
4491 TGraphErrors *graphPRF0 = new TGraphErrors(nbins,xValues,fCoefPRF[0],xValuesE,fCoefPRFE[0]);
4492 graphPRF0->SetName("coefprf0");
4493 graphPRF0->SetTitle("");
4494 graphPRF0->GetXaxis()->SetTitle("Det/Pad groups");
4495 graphPRF0->GetYaxis()->SetTitle("PRF Width [p.u]");
4496 graphPRF0->SetMarkerColor(6);
4497 graphPRF0->SetLineColor(6);
4498 graphPRF0->SetMarkerStyle(26);
4499 listofgraphs->Add((TObject *)graphPRF0);
4500 legprf1->AddEntry(graphPRF0,"PRF fit","p");
4501 graphPRF0->Draw("P");
4504 TGraphErrors *graphPRF2 = new TGraphErrors(nbins,xValues,fCoefPRF[2],xValuesE,fCoefPRFE[1]);
4505 graphPRF2->SetName("coefprf2");
4506 graphPRF2->SetTitle("");
4507 graphPRF2->GetXaxis()->SetTitle("Det/Pad groups");
4508 graphPRF2->GetYaxis()->SetTitle("PRF Width [p.u]");
4509 graphPRF2->SetMarkerColor(1);
4510 graphPRF2->SetLineColor(1);
4511 graphPRF2->SetMarkerStyle(28);
4512 listofgraphs->Add((TObject *)graphPRF2);
4513 legprf1->AddEntry(graphPRF2,"PRF Rms","p");
4514 graphPRF2->Draw("P");
4516 legprf1->Draw("same");
4519 //Create the arrays and the graphs for the delta
4520 TCanvas *cprf2 = new TCanvas("cprf2","",50,50,600,800);
4523 TLegend *legprf3 = new TLegend(0.4,0.6,0.89,0.89);
4524 TLegend *legprf2 = new TLegend(0.4,0.6,0.89,0.89);
4528 Double_t *yValuesDelta = new Double_t[counter[0]];
4529 for(Int_t k = 0; k < counter[0]; k++){
4530 if(fCoefPRF[1][(Int_t)xValuesFitted[k]] > 0.0){
4531 yValuesDelta[k] = (fCoefPRF[0][(Int_t)xValuesFitted[k]]-fCoefPRF[1][(Int_t)xValuesFitted[k]])
4532 / (fCoefPRF[1][(Int_t)xValuesFitted[k]]);
4535 TGraph *graphDeltaPRF0 = new TGraph(counter[0],&xValuesFitted[0],yValuesDelta);
4536 graphDeltaPRF0->SetName("deltaprf0");
4537 graphDeltaPRF0->GetXaxis()->SetTitle("Det/Pad groups");
4538 graphDeltaPRF0->GetYaxis()->SetTitle("#Delta#sigma/#sigma_{sim}");
4539 graphDeltaPRF0->SetMarkerColor(6);
4540 graphDeltaPRF0->SetTitle("");
4541 graphDeltaPRF0->SetLineColor(6);
4542 graphDeltaPRF0->SetMarkerStyle(26);
4543 listofgraphs->Add((TObject *)graphDeltaPRF0);
4544 legprf3->AddEntry(graphDeltaPRF0,"#sigma_{fit}","p");
4545 graphDeltaPRF0->Draw("AP");
4548 TH1I *histoErrorPRF0 = new TH1I("errorprf10","",100 ,-0.1,0.2);
4549 histoErrorPRF0->SetXTitle("#Delta#sigma/#sigma_{sim}");
4550 histoErrorPRF0->SetYTitle("counts");
4551 histoErrorPRF0->SetLineColor(6);
4552 histoErrorPRF0->SetLineStyle(1);
4553 histoErrorPRF0->SetStats(0);
4554 Double_t maxvalue = 0.0;
4555 for(Int_t k = 0; k < counter[0]; k++){
4556 histoErrorPRF0->Fill(yValuesDelta[k]);
4557 if(k == 0) maxvalue = yValuesDelta[k];
4558 if(maxvalue < (TMath::Abs(yValuesDelta[k]))) maxvalue = (TMath::Abs(yValuesDelta[k]));
4560 AliInfo(Form("The maximum deviation for the fit method is %f",maxvalue));
4561 legprf2->AddEntry(histoErrorPRF0,"#sigma_{fit}","l");
4562 histoErrorPRF0->Draw();
4563 listofgraphs->Add((TObject *)histoErrorPRF0);
4569 Double_t *yValuesDelta = new Double_t[counter[1]];
4570 for(Int_t k = 0; k < counter[1]; k++){
4571 if(fCoefPRF[1][(Int_t)xValuesRMS[k]] > 0.0){
4572 yValuesDelta[k] = (fCoefPRF[2][(Int_t)xValuesRMS[k]]-fCoefPRF[1][(Int_t)xValuesRMS[k]])
4573 / (fCoefPRF[1][(Int_t)xValuesRMS[k]]);
4576 TGraph *graphDeltaPRF2 = new TGraph(counter[1],&xValuesRMS[0],yValuesDelta);
4577 graphDeltaPRF2->SetName("deltaprf2");
4578 graphDeltaPRF2->GetXaxis()->SetTitle("Det/Pad groups");
4579 graphDeltaPRF2->GetYaxis()->SetTitle("#Delta#sigma/#sigma_{sim}");
4580 graphDeltaPRF2->SetMarkerColor(1);
4581 graphDeltaPRF2->SetTitle("");
4582 graphDeltaPRF2->SetLineColor(1);
4583 graphDeltaPRF2->SetMarkerStyle(28);
4584 listofgraphs->Add((TObject *)graphDeltaPRF2);
4585 legprf3->AddEntry(graphDeltaPRF2,"#sigma_{rms}","p");
4587 graphDeltaPRF2->Draw("P");
4590 graphDeltaPRF2->Draw("AP");
4594 TH1I *histoErrorPRF2 = new TH1I("errorprf12","",100 ,-0.1,0.2);
4595 histoErrorPRF2->SetXTitle("#Delta#sigma/#sigma_{sim}");
4596 histoErrorPRF2->SetYTitle("counts");
4597 histoErrorPRF2->SetLineColor(1);
4598 histoErrorPRF2->SetLineStyle(1);
4599 histoErrorPRF2->SetStats(0);
4600 Double_t maxvalue = 0.0;
4601 for(Int_t k = 0; k < counter[1]; k++){
4602 histoErrorPRF2->Fill(yValuesDelta[k]);
4603 if(k == 0) maxvalue = yValuesDelta[k];
4604 if(maxvalue < TMath::Abs(yValuesDelta[k])) maxvalue = TMath::Abs(yValuesDelta[k]);
4606 AliInfo(Form("The maximum deviation for the rms is %f",maxvalue));
4607 legprf2->AddEntry(histoErrorPRF2,"#sigma_{rms}","l");
4609 histoErrorPRF2->Draw("same");
4612 histoErrorPRF2->Draw();
4614 listofgraphs->Add((TObject *)histoErrorPRF2);
4619 legprf3->Draw("same");
4621 legprf2->Draw("same");
4626 TFile *fout = TFile::Open(fWriteNameCoef,"UPDATE");
4627 // Check if the file could be opened
4628 if (!fout || !fout->IsOpen()) {
4629 AliInfo("No File found!");
4633 for(Int_t k = 0; k < listofgraphs->GetEntriesFast(); k++){
4634 fout->WriteTObject((TObject *) listofgraphs->At(k),((TObject *)listofgraphs->At(k))->GetName()
4635 ,(Option_t *) "OverWrite");
4644 //____________Plot histos DB___________________________________________________
4647 //_____________________________________________________________________________
4648 void AliTRDCalibraFit::PlotCHDB()
4651 // Plot the histos for fDebug = 3 and fDebug = 4 to visualise the detector
4654 TCanvas *cchdb = new TCanvas("cchdb","",50,50,600,800);
4656 if(fFitChargeOn) nb++;
4657 if(fFitChargeBisOn) nb++;
4658 if(fMeanChargeOn) nb++;
4659 if(fFitMeanWOn) nb++;
4661 cchdb->Divide(nb,1);
4665 fCoefChargeDB[1]->Draw("LEGO");
4670 fCoefChargeDB[0]->Draw("LEGO");
4675 fCoefChargeDB[3]->Draw("LEGO");
4678 if(fFitChargeBisOn){
4680 fCoefChargeDB[2]->Draw("LEGO");
4681 //it doesn't matter anymore but....
4687 //_____________________________________________________________________________
4688 void AliTRDCalibraFit::PlotPHDB()
4691 // Plot the histos for fDebug = 3 and fDebug = 4 to visualise the detector
4694 TCanvas *cphdb = new TCanvas("cphdb","",50,50,600,800);
4696 if(fFitPol2On) nb++;
4698 if(fFitLagrPolOn) nb++;
4700 cphdb->Divide(nb,1);
4704 fCoefVdriftDB[0]->Draw("LEGO");
4709 fCoefVdriftDB[1]->Draw("LEGO");
4714 fCoefVdriftDB[2]->Draw("LEGO");
4720 //_____________________________________________________________________________
4721 void AliTRDCalibraFit::PlotT0DB()
4724 // Plot the histos for fDebug = 3 and fDebug = 4 to visualise the detector
4726 TCanvas *ct0db = new TCanvas("ct0db","",50,50,600,800);
4728 if(fFitPol2On) nb++;
4730 if(fFitLagrPolOn) nb++;
4732 ct0db->Divide(nb,1);
4736 fCoefT0DB[0]->Draw("LEGO");
4741 fCoefT0DB[1]->Draw("LEGO");
4746 fCoefT0DB[2]->Draw("LEGO");
4752 //_____________________________________________________________________________
4753 void AliTRDCalibraFit::PlotPRFDB()
4756 // Plot the histos for fDebug = 3 and fDebug = 4 to visualise the detector
4759 TCanvas *cprfdb = new TCanvas("cprfdb","",50,50,600,800);
4764 cprfdb->Divide(nb,1);
4768 fCoefPRFDB[0]->Draw("LEGO");
4773 fCoefPRFDB[1]->Draw("LEGO");
4780 //____________Write DB Histos__________________________________________________
4783 //_____________________________________________________________________________
4784 void AliTRDCalibraFit::WriteCHDB(TFile *fout)
4787 // If wanted, write the debug histos for fDebug = 3 and fDebug = 4
4790 fout->WriteTObject(fCoefChargeDB[0],fCoefChargeDB[0]->GetName(),(Option_t *) "OverWrite");
4793 fout->WriteTObject(fCoefChargeDB[3],fCoefChargeDB[0]->GetName(),(Option_t *) "OverWrite");
4795 if (fMeanChargeOn) {
4796 fout->WriteTObject(fCoefChargeDB[1],fCoefChargeDB[1]->GetName(),(Option_t *) "OverWrite");
4798 if (fFitChargeBisOn ) {
4799 fout->WriteTObject(fCoefChargeDB[2],fCoefChargeDB[2]->GetName(),(Option_t *) "OverWrite");
4804 //_____________________________________________________________________________
4805 void AliTRDCalibraFit::WritePHDB(TFile *fout)
4808 // If wanted, write the debug histos for fDebug = 3 and fDebug = 4
4812 fout->WriteTObject(fCoefVdriftDB[0],fCoefVdriftDB[0]->GetName(),(Option_t *) "OverWrite");
4815 fout->WriteTObject(fCoefVdriftDB[1],fCoefVdriftDB[1]->GetName(),(Option_t *) "OverWrite");
4818 fout->WriteTObject(fCoefVdriftDB[2],fCoefVdriftDB[2]->GetName(),(Option_t *) "OverWrite");
4823 //_____________________________________________________________________________
4824 void AliTRDCalibraFit::WriteT0DB(TFile *fout)
4827 // If wanted, write the debug histos for fDebug = 3 and fDebug = 4
4831 fout->WriteTObject(fCoefT0DB[0],fCoefT0DB[0]->GetName(),(Option_t *) "OverWrite");
4834 fout->WriteTObject(fCoefT0DB[1],fCoefT0DB[1]->GetName(),(Option_t *) "OverWrite");
4837 fout->WriteTObject(fCoefT0DB[2],fCoefT0DB[2]->GetName(),(Option_t *) "OverWrite");
4842 //_____________________________________________________________________________
4843 void AliTRDCalibraFit::WritePRFDB(TFile *fout)
4846 // If wanted, write the debug histos for fDebug = 3 and fDebug = 4
4849 fout->WriteTObject(fCoefPRFDB[0],fCoefPRFDB[0]->GetName(),(Option_t *) "OverWrite");
4852 fout->WriteTObject(fCoefPRFDB[1],fCoefPRFDB[1]->GetName(),(Option_t *) "OverWrite");
4858 //____________Calcul Coef Mean_________________________________________________
4861 //_____________________________________________________________________________
4862 Bool_t AliTRDCalibraFit::CalculT0CoefMean(Int_t dect, Int_t idect)
4865 // For the detector Dect calcul the mean time 0
4866 // for the calibration group idect from the choosen database
4869 AliTRDcalibDB *cal = AliTRDcalibDB::Instance();
4871 AliInfo("Could not get calibDB Manager");
4877 if ((fDebug != 2) && fAccCDB) {
4879 for (Int_t row = fCalibraMode->GetRowMin(1); row < fCalibraMode->GetRowMax(1); row++) {
4880 for (Int_t col = fCalibraMode->GetColMin(1); col < fCalibraMode->GetColMax(1); col++) {
4882 if ((fCalibraMode->GetNz(1) > 0) ||
4883 (fCalibraMode->GetNrphi(1) > 0)) {
4884 fT0Coef[2] += (Float_t) cal->GetT0(dect,col,row);
4888 fT0Coef[2] += (Float_t) cal->GetT0Average(dect);
4893 fT0Coef[2] = fT0Coef[2] / ((fCalibraMode->GetColMax(1)-fCalibraMode->GetColMin(1))
4894 * (fCalibraMode->GetRowMax(1)-fCalibraMode->GetRowMin(1)));
4895 if ((fDebug == 1) ||
4897 fCoefT0[2][idect] = fT0Coef[2];
4906 //_____________________________________________________________________________
4907 Bool_t AliTRDCalibraFit::CalculChargeCoefMean(Int_t dect, Int_t idect, Bool_t vrai)
4910 // For the detector Dect calcul the mean gain factor
4911 // for the calibration group idect from the choosen database
4914 AliTRDcalibDB *cal = AliTRDcalibDB::Instance();
4916 AliInfo("Could not get calibDB Manager");
4920 fChargeCoef[3] = 0.0;
4924 for (Int_t row = fCalibraMode->GetRowMin(0); row < fCalibraMode->GetRowMax(0); row++) {
4925 for (Int_t col = fCalibraMode->GetColMin(0); col < fCalibraMode->GetColMax(0); col++) {
4927 if ((fCalibraMode->GetNz(0) > 0) ||
4928 (fCalibraMode->GetNrphi(0) > 0)) {
4930 fChargeCoef[3] += (Float_t) cal->GetGainFactor(dect,col,row);
4932 if (vrai && fAccCDB) {
4933 fScaleFitFactor += (Float_t) cal->GetGainFactor(dect,col,row);
4936 fChargeCoef[3] += 1.0;
4938 if (vrai && (!fAccCDB)) {
4939 fScaleFitFactor += 1.0;
4945 fChargeCoef[3] += (Float_t) cal->GetGainFactorAverage(dect);
4947 if (vrai && fAccCDB) {
4948 fScaleFitFactor += ((Float_t) cal->GetGainFactorAverage(dect));
4951 fChargeCoef[3] += 1.0;
4953 if (vrai && (!fAccCDB)) {
4954 fScaleFitFactor += 1.0;
4960 fChargeCoef[3] = fChargeCoef[3] / ((fCalibraMode->GetColMax(0)-fCalibraMode->GetColMin(0))
4961 * (fCalibraMode->GetRowMax(0)-fCalibraMode->GetRowMin(0)));
4962 if ((fDebug == 1) ||
4964 fCoefCharge[3][idect]=fChargeCoef[3];
4973 //_____________________________________________________________________________
4974 Bool_t AliTRDCalibraFit::CalculPRFCoefMean(Int_t dect, Int_t idect)
4977 // For the detector Dect calcul the mean sigma of pad response
4978 // function for the calibration group idect from the choosen database
4981 AliTRDcalibDB *cal = AliTRDcalibDB::Instance();
4983 AliInfo("Could not get calibDB Manager");
4992 for (Int_t row = fCalibraMode->GetRowMin(2); row < fCalibraMode->GetRowMax(2); row++) {
4993 for (Int_t col = fCalibraMode->GetColMin(2); col < fCalibraMode->GetColMax(2); col++) {
4994 if ((fGeo->GetColMax(GetPlane(dect)) != (col+1)) && (col != 0)) {
4997 fPRFCoef[1] += (Float_t) cal->GetPRFWidth(dect,col,row);
5000 fPRFCoef[1] += GetPRFDefault(GetPlane(dect));
5007 fPRFCoef[1] = fPRFCoef[1]/cot;
5008 if ((fDebug == 1) ||
5010 fCoefPRF[1][idect] = fPRFCoef[1];
5014 if ((fDebug == 1) ||
5017 fCoefPRF[1][idect] = cal->GetPRFWidth(dect,fCalibraMode->GetColMin(2),fCalibraMode->GetRowMin(2));
5020 fCoefPRF[1][idect] = GetPRFDefault(GetPlane(dect));
5031 //_____________________________________________________________________________
5032 Bool_t AliTRDCalibraFit::CalculVdriftCoefMean(Int_t dect, Int_t idect)
5035 // For the detector dect calcul the mean drift velocity for the
5036 // calibration group idect from the choosen database
5039 AliTRDcalibDB *cal = AliTRDcalibDB::Instance();
5041 AliInfo("Could not get calibDB Manager");
5045 fVdriftCoef[2] = 0.0;
5048 for (Int_t row = fCalibraMode->GetRowMin(1); row < fCalibraMode->GetRowMax(1); row++) {
5049 for (Int_t col = fCalibraMode->GetColMin(1); col < fCalibraMode->GetColMax(1); col++) {
5051 if ((fCalibraMode->GetNz(1) > 0) ||
5052 (fCalibraMode->GetNrphi(1) > 0)) {
5054 fVdriftCoef[2] += (Float_t) cal->GetVdrift(dect,col,row);
5057 fVdriftCoef[2] += 1.5;
5063 fVdriftCoef[2] += (Float_t) cal->GetVdriftAverage(dect);
5066 fVdriftCoef[2] += 1.5;
5071 fVdriftCoef[2] = fVdriftCoef[2] / ((fCalibraMode->GetColMax(1)-fCalibraMode->GetColMin(1))
5072 * (fCalibraMode->GetRowMax(1)-fCalibraMode->GetRowMin(1)));
5073 if ((fDebug == 1) ||
5075 fCoefVdrift[2][idect] = fVdriftCoef[2];
5083 //_____________________________________________________________________________
5084 Float_t AliTRDCalibraFit::GetPRFDefault(Int_t plane) const
5087 // Default width of the PRF if there is no database as reference
5114 //____________Fit Methods______________________________________________________
5116 //_____________________________________________________________________________
5117 void AliTRDCalibraFit::FitPente(TH1* projPH, Int_t idect)
5120 // Slope methode for the drift velocity
5124 const Float_t kDrWidth = AliTRDgeometry::DrThick();
5131 Double_t vdriftCoefE = 0.0;
5132 Double_t t0CoefE = 0.0;
5133 fVdriftCoef[1] = 0.0;
5135 TLine *line = new TLine();
5138 TAxis *xpph = projPH->GetXaxis();
5139 Int_t nbins = xpph->GetNbins();
5140 Double_t lowedge = xpph->GetBinLowEdge(1);
5141 Double_t upedge = xpph->GetBinUpEdge(xpph->GetNbins());
5142 Double_t widbins = (upedge - lowedge) / nbins;
5143 Double_t limit = upedge + 0.5 * widbins;
5146 // Beginning of the signal
5147 TH1D *pentea = new TH1D("pentea","pentea",projPH->GetNbinsX(),0,(Float_t) limit);
5148 for (Int_t k = 1; k < projPH->GetNbinsX(); k++) {
5149 pentea->SetBinContent(k,(Double_t) (projPH->GetBinContent(k+1) - projPH->GetBinContent(k)));
5152 binmax = (Int_t) pentea->GetMaximumBin();
5153 if(fDebug == 2) AliInfo(Form("maximum positive bin for the positive slope %d",binmax));
5156 AliInfo("Put the binmax from 1 to 2 to enable the fit");
5158 if (binmax >= nbins) {
5161 AliInfo("Put the binmax from nbins-1 to nbins-2 to enable the fit");
5163 pentea->Fit("pol2","0MR","",TMath::Max(pentea->GetBinCenter(binmax-1),0.0),pentea->GetBinCenter(binmax+1));
5164 Float_t l3P1am = pentea->GetFunction("pol2")->GetParameter(1);
5165 Float_t l3P2am = pentea->GetFunction("pol2")->GetParameter(2);
5166 Float_t l3P1amE = pentea->GetFunction("pol2")->GetParError(1);
5167 Float_t l3P2amE = pentea->GetFunction("pol2")->GetParError(2);
5169 fPhd[0] = -(l3P1am / (2 * l3P2am));
5172 if((l3P1am != 0.0) && (l3P2am != 0.0)){
5173 t0CoefE = (l3P1amE/l3P1am + l3P2amE/l3P2am)*fPhd[0];
5176 if(fDebug == 2) AliInfo(Form("maximum extrapolated positive bin for the positive slope %f",fPhd[0]));
5178 // Amplification region
5181 for (Int_t kbin = 1; kbin < projPH->GetNbinsX(); kbin ++) {
5182 if (((projPH->GetBinContent(kbin+1) - projPH->GetBinContent(kbin)) <= 0.0) &&
5184 (kbin > (fPhd[0]/widbins))) {
5189 if(fDebug == 2) AliInfo(Form("For the amplification region binmax %d",binmax));
5192 AliInfo("Put the binmax from 1 to 2 to enable the fit");
5194 if (binmax >= nbins) {
5197 AliInfo("Put the binmax from nbins-1 to nbins-2 to enable the fit");
5199 projPH->Fit("pol2","0MR","",TMath::Max(projPH->GetBinCenter(binmax-1),0.0),projPH->GetBinCenter(binmax+1));
5200 Float_t l3P1amf = projPH->GetFunction("pol2")->GetParameter(1);
5201 Float_t l3P2amf = projPH->GetFunction("pol2")->GetParameter(2);
5202 Float_t l3P1amfE = projPH->GetFunction("pol2")->GetParError(1);
5203 Float_t l3P2amfE = projPH->GetFunction("pol2")->GetParError(2);
5206 fPhd[1] = -(l3P1amf / (2 * l3P2amf));
5208 if((l3P1amf != 0.0) && (l3P2amf != 0.0)){
5209 vdriftCoefE = (l3P1amfE/l3P1amf + l3P2amfE/l3P2amf)*fPhd[1];
5212 t0CoefE = vdriftCoefE;
5214 if(fDebug == 2) AliInfo(Form("For the amplification region extrapolated binmax %f",fPhd[1]));
5217 TH1D *pente = new TH1D("pente","pente",projPH->GetNbinsX(),0,(Float_t) limit);
5218 for (Int_t k = TMath::Min(binmax+4,projPH->GetNbinsX()); k < projPH->GetNbinsX(); k++) {
5219 pente->SetBinContent(k,(Double_t) (projPH->GetBinContent(k+1) - projPH->GetBinContent(k)));
5222 if(pente->GetEntries() > 0) binmin = (Int_t) pente->GetMinimumBin();
5225 AliInfo("Put the binmax from 1 to 2 to enable the fit");
5227 if (binmin >= nbins) {
5230 AliInfo("Put the binmax from nbins-1 to nbins-2 to enable the fit");
5232 if(fDebug == 2) AliInfo(Form("For the drift region binmin %d",binmin));
5236 ,TMath::Max(pente->GetBinCenter(binmin-1), 0.0)
5237 ,TMath::Min(pente->GetBinCenter(binmin+1),(Double_t) limit));
5238 Float_t l3P1dr = pente->GetFunction("pol2")->GetParameter(1);
5239 Float_t l3P2dr = pente->GetFunction("pol2")->GetParameter(2);
5240 Float_t l3P1drE = pente->GetFunction("pol2")->GetParError(1);
5241 Float_t l3P2drE = pente->GetFunction("pol2")->GetParError(2);
5243 fPhd[2] = -(l3P1dr / (2 * l3P2dr));
5245 if((l3P1dr != 0.0) && (l3P2dr != 0.0)){
5246 vdriftCoefE += (l3P1drE/l3P1dr + l3P2drE/l3P2dr)*fPhd[2];
5248 if(fDebug == 2) AliInfo(Form("For the drift region extrapolated binmax %f",fPhd[2]));
5250 Float_t fPhdt0 = 0.0;
5251 if(fTakeTheMaxPH) fPhdt0 = fPhd[1];
5252 else fPhdt0 = fPhd[0];
5254 if ((fPhd[2] > fPhd[0]) &&
5255 (fPhd[2] > fPhd[1]) &&
5256 (fPhd[1] > fPhd[0]) &&
5258 fVdriftCoef[1] = (kDrWidth) / (fPhd[2]-fPhd[1]);
5259 if(fFitPHNDB == 1) fNumberFitSuccess++;
5260 if (fPhdt0 >= 0.0) {
5261 fT0Coef[1] = (fPhdt0 - fT0Shift) / widbins;
5262 if (fT0Coef[1] < -1.0) {
5263 fT0Coef[1] = fT0Coef[2];
5267 fT0Coef[1] = fT0Coef[2];
5271 fVdriftCoef[1] = -TMath::Abs(fVdriftCoef[2]);
5272 fT0Coef[1] = fT0Coef[2];
5275 if ((fDebug == 1) ||
5277 fCoefVdrift[1][idect] = fVdriftCoef[1];
5278 fCoefVdriftE[1] [idect] = vdriftCoefE;
5279 fCoefT0[1][idect] = fT0Coef[1];
5280 fCoefT0E[1][idect] = t0CoefE;
5284 TCanvas *cpentei = new TCanvas("cpentei","cpentei",50,50,600,800);
5287 line->SetLineColor(2);
5288 line->DrawLine(fPhd[0],0,fPhd[0],projPH->GetMaximum());
5289 line->DrawLine(fPhd[1],0,fPhd[1],projPH->GetMaximum());
5290 line->DrawLine(fPhd[2],0,fPhd[2],projPH->GetMaximum());
5291 AliInfo(Form("fPhd[0] (beginning of the signal): %f" ,(Float_t) fPhd[0]));
5292 AliInfo(Form("fPhd[1] (end of the amplification region): %f" ,(Float_t) fPhd[1]));
5293 AliInfo(Form("fPhd[2] (end of the drift region): %f" ,(Float_t) fPhd[2]));
5294 AliInfo(Form("fVriftCoef[1] (with only the drift region(default)): %f",(Float_t) fVdriftCoef[1]));
5295 TCanvas *cpentei2 = new TCanvas("cpentei2","cpentei2",50,50,600,800);
5298 TCanvas *cpentei3 = new TCanvas("cpentei3","cpentei3",50,50,600,800);
5312 //_____________________________________________________________________________
5313 void AliTRDCalibraFit::FitLagrangePoly(TH1* projPH, Int_t idect)
5316 // Slope methode but with polynomes de Lagrange
5320 const Float_t kDrWidth = AliTRDgeometry::DrThick();
5323 Double_t *x = new Double_t[5];
5324 Double_t *y = new Double_t[5];
5339 Double_t vdriftCoefE = 0.0;
5340 Double_t t0CoefE = 1.0;
5341 fVdriftCoef[3] = 0.0;
5343 TLine *line = new TLine();
5344 TF1 * polynome = 0x0;
5345 TF1 * polynomea = 0x0;
5346 TF1 * polynomeb = 0x0;
5350 TAxis *xpph = projPH->GetXaxis();
5351 Int_t nbins = xpph->GetNbins();
5352 Double_t lowedge = xpph->GetBinLowEdge(1);
5353 Double_t upedge = xpph->GetBinUpEdge(xpph->GetNbins());
5354 Double_t widbins = (upedge - lowedge) / nbins;
5355 Double_t limit = upedge + 0.5 * widbins;
5360 // Beginning of the signal
5361 TH1D *pentea = new TH1D("pentea","pentea",projPH->GetNbinsX(),0,(Float_t) limit);
5362 for (Int_t k = 1; k < projPH->GetNbinsX(); k++) {
5363 pentea->SetBinContent(k,(Double_t) (projPH->GetBinContent(k+1) - projPH->GetBinContent(k)));
5366 binmax = (Int_t) pentea->GetMaximumBin();
5367 if(fDebug == 2) AliInfo(Form("maximum positive bin for the positive slope %d",binmax));
5370 Double_t minnn = 0.0;
5371 Double_t maxxx = 0.0;
5373 //Determination of minnn and maxxx
5374 //case binmax = nbins -1
5376 if(binmax == (nbins-1)){
5377 minnn = pentea->GetBinCenter(binmax-2);
5378 maxxx = pentea->GetBinCenter(binmax);
5379 x[0] = pentea->GetBinCenter(binmax-2);
5380 x[1] = pentea->GetBinCenter(binmax-1);
5381 x[2] = pentea->GetBinCenter(binmax);
5382 y[0] = pentea->GetBinContent(binmax-2);
5383 y[1] = pentea->GetBinContent(binmax-1);
5384 y[2] = pentea->GetBinContent(binmax);
5385 //Calcul the polynome de Lagrange
5386 c = CalculPolynomeLagrange2(x,y);
5387 AliInfo("At the limit for beginning!");
5389 //case binmax = nbins-2
5391 if(binmax == (nbins-2)){
5392 minnn = pentea->GetBinCenter(binmax-2);
5393 maxxx = pentea->GetBinCenter(binmax+1);
5394 x[0] = pentea->GetBinCenter(binmax-2);
5395 x[1] = pentea->GetBinCenter(binmax-1);
5396 x[2] = pentea->GetBinCenter(binmax);
5397 x[3] = pentea->GetBinCenter(binmax+1);
5398 y[0] = pentea->GetBinContent(binmax-2);
5399 y[1] = pentea->GetBinContent(binmax-1);
5400 y[2] = pentea->GetBinContent(binmax);
5401 y[3] = pentea->GetBinContent(binmax+1);
5402 //Calcul the polynome de Lagrange
5403 c = CalculPolynomeLagrange3(x,y);
5405 //case binmax <= nbins-3
5407 if(binmax <= (nbins-3)){
5408 if((binmax-2) >= 1){
5409 minnn = pentea->GetBinCenter(binmax-2);
5410 maxxx = pentea->GetBinCenter(binmax+2);
5411 x[0] = pentea->GetBinCenter(binmax-2);
5412 x[1] = pentea->GetBinCenter(binmax-1);
5413 x[2] = pentea->GetBinCenter(binmax);
5414 x[3] = pentea->GetBinCenter(binmax+1);
5415 x[4] = pentea->GetBinCenter(binmax+2);
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 y[4] = pentea->GetBinContent(binmax+2);
5421 //Calcul the polynome de Lagrange
5422 c = CalculPolynomeLagrange4(x,y);
5425 if((binmax-1) == 1){
5426 minnn = pentea->GetBinCenter(binmax-1);
5427 maxxx = pentea->GetBinCenter(binmax+2);
5428 x[0] = pentea->GetBinCenter(binmax-1);
5429 x[1] = pentea->GetBinCenter(binmax);
5430 x[2] = pentea->GetBinCenter(binmax+1);
5431 x[3] = pentea->GetBinCenter(binmax+2);
5432 y[0] = pentea->GetBinContent(binmax-1);
5433 y[1] = pentea->GetBinContent(binmax);
5434 y[2] = pentea->GetBinContent(binmax+1);
5435 y[3] = pentea->GetBinContent(binmax+2);
5436 //Calcul the polynome de Lagrange
5437 c = CalculPolynomeLagrange3(x,y);
5441 minnn = pentea->GetBinCenter(binmax);
5442 maxxx = pentea->GetBinCenter(binmax+2);
5443 x[0] = pentea->GetBinCenter(binmax);
5444 x[1] = pentea->GetBinCenter(binmax+1);
5445 x[2] = pentea->GetBinCenter(binmax+2);
5446 y[0] = pentea->GetBinContent(binmax);
5447 y[1] = pentea->GetBinContent(binmax+1);
5448 y[2] = pentea->GetBinContent(binmax+2);
5449 //Calcul the polynome de Lagrange
5450 c = CalculPolynomeLagrange2(x,y);
5453 //pass but should not happen
5454 if((binmax <= (nbins-3)) && (binmax < 1)){
5458 if(fDebug == 2) AliInfo(Form("For the beginning region binmax %d",binmax));
5461 polynomeb = new TF1("polb","[0]+[1]*x+[2]*x*x+[3]*x*x*x+[4]*x*x*x*x",minnn,maxxx);
5462 polynomeb->SetParameters(c[0],c[1],c[2],c[3],c[4]);
5464 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]));
5467 Double_t step = (maxxx-minnn)/10000;
5469 Double_t maxvalue = 0.0;
5470 Double_t placemaximum = minnn;
5471 for(Int_t o = 0; o < 10000; o++){
5472 if(o == 0) maxvalue = polynomeb->Eval(l);
5473 if(maxvalue < (polynomeb->Eval(l))){
5474 maxvalue = polynomeb->Eval(l);
5479 fPhd[0] = placemaximum;
5482 if(fDebug == 2) AliInfo(Form("maximum extrapolated positive bin for the positive slope %f",fPhd[0]));
5484 // Amplification region
5487 for (Int_t kbin = 1; kbin < projPH->GetNbinsX(); kbin ++) {
5488 if (((projPH->GetBinContent(kbin+1) - projPH->GetBinContent(kbin)) <= 0.0) &&
5490 (kbin > (fPhd[0]/widbins))) {
5495 if(fDebug == 2) AliInfo(Form("For the amplification region binmax %d",binmax));
5497 Double_t minn = 0.0;
5498 Double_t maxx = 0.0;
5500 //Determination of minn and maxx
5501 //case binmax = nbins
5503 if(binmax == nbins){
5504 minn = projPH->GetBinCenter(binmax-2);
5505 maxx = projPH->GetBinCenter(binmax);
5506 x[0] = projPH->GetBinCenter(binmax-2);
5507 x[1] = projPH->GetBinCenter(binmax-1);
5508 x[2] = projPH->GetBinCenter(binmax);
5509 y[0] = projPH->GetBinContent(binmax-2);
5510 y[1] = projPH->GetBinContent(binmax-1);
5511 y[2] = projPH->GetBinContent(binmax);
5512 //Calcul the polynome de Lagrange
5513 c = CalculPolynomeLagrange2(x,y);
5514 AliInfo("At the limit for the drift!");
5516 //case binmax = nbins-1
5518 if(binmax == (nbins-1)){
5519 minn = projPH->GetBinCenter(binmax-2);
5520 maxx = projPH->GetBinCenter(binmax+1);
5521 x[0] = projPH->GetBinCenter(binmax-2);
5522 x[1] = projPH->GetBinCenter(binmax-1);
5523 x[2] = projPH->GetBinCenter(binmax);
5524 x[3] = projPH->GetBinCenter(binmax+1);
5525 y[0] = projPH->GetBinContent(binmax-2);
5526 y[1] = projPH->GetBinContent(binmax-1);
5527 y[2] = projPH->GetBinContent(binmax);
5528 y[3] = projPH->GetBinContent(binmax+1);
5529 //Calcul the polynome de Lagrange
5530 c = CalculPolynomeLagrange3(x,y);
5532 //case binmax <= nbins-2
5534 if(binmax <= (nbins-2)){
5535 if((binmax-2) >= 1){
5536 minn = projPH->GetBinCenter(binmax-2);
5537 maxx = projPH->GetBinCenter(binmax+2);
5538 x[0] = projPH->GetBinCenter(binmax-2);
5539 x[1] = projPH->GetBinCenter(binmax-1);
5540 x[2] = projPH->GetBinCenter(binmax);
5541 x[3] = projPH->GetBinCenter(binmax+1);
5542 x[4] = projPH->GetBinCenter(binmax+2);
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 y[4] = projPH->GetBinContent(binmax+2);
5548 //Calcul the polynome de Lagrange
5549 c = CalculPolynomeLagrange4(x,y);
5552 if((binmax-1) == 1){
5553 minn = projPH->GetBinCenter(binmax-1);
5554 maxx = projPH->GetBinCenter(binmax+2);
5555 x[0] = projPH->GetBinCenter(binmax-1);
5556 x[1] = projPH->GetBinCenter(binmax);
5557 x[2] = projPH->GetBinCenter(binmax+1);
5558 x[3] = projPH->GetBinCenter(binmax+2);
5559 y[0] = projPH->GetBinContent(binmax-1);
5560 y[1] = projPH->GetBinContent(binmax);
5561 y[2] = projPH->GetBinContent(binmax+1);
5562 y[3] = projPH->GetBinContent(binmax+2);
5563 //Calcul the polynome de Lagrange
5564 c = CalculPolynomeLagrange3(x,y);
5568 minn = projPH->GetBinCenter(binmax);
5569 maxx = projPH->GetBinCenter(binmax+2);
5570 x[0] = projPH->GetBinCenter(binmax);
5571 x[1] = projPH->GetBinCenter(binmax+1);
5572 x[2] = projPH->GetBinCenter(binmax+2);
5573 y[0] = projPH->GetBinContent(binmax);
5574 y[1] = projPH->GetBinContent(binmax+1);
5575 y[2] = projPH->GetBinContent(binmax+2);
5576 //Calcul the polynome de Lagrange
5577 c = CalculPolynomeLagrange2(x,y);
5580 //pass but should not happen
5581 if((binmax <= (nbins-2)) && (binmax < 1)){
5585 if(fDebug == 2) AliInfo(Form("For the amplification region binmax %d",binmax));
5588 polynomea = new TF1("pola","[0]+[1]*x+[2]*x*x+[3]*x*x*x+[4]*x*x*x*x",minn,maxx);
5589 polynomea->SetParameters(c[0],c[1],c[2],c[3],c[4]);
5591 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]));
5594 Double_t step = (maxx-minn)/1000;
5596 Double_t maxvalue = 0.0;
5597 Double_t placemaximum = minn;
5598 for(Int_t o = 0; o < 1000; o++){
5599 if(o == 0) maxvalue = polynomea->Eval(l);
5600 if(maxvalue < (polynomea->Eval(l))){
5601 maxvalue = polynomea->Eval(l);
5606 fPhd[1] = placemaximum;
5609 if(fDebug == 2) AliInfo(Form("For the amplification region extrapolated binmax %f",fPhd[1]));
5612 TH1D *pente = new TH1D("pente","pente", projPH->GetNbinsX(),0,(Float_t) limit);
5613 for (Int_t k = TMath::Min(binmax+4, projPH->GetNbinsX()); k < projPH->GetNbinsX(); k++) {
5614 pente->SetBinContent(k,(Double_t) (projPH->GetBinContent(k+1) - projPH->GetBinContent(k)));
5617 if(pente->GetEntries() > 0) binmin = (Int_t) pente->GetMinimumBin();
5623 AliInfo("Put the binmax from 1 to 2 to enable the fit");
5627 if((projPH->GetBinContent(binmin)-projPH->GetBinError(binmin)) < (projPH->GetBinContent(binmin+1))) put = kFALSE;
5628 if((projPH->GetBinContent(binmin)+projPH->GetBinError(binmin)) > (projPH->GetBinContent(binmin-1))) put = kFALSE;
5631 AliInfo(Form("binmin %d BinContent %f BinError %f",binmin
5632 ,projPH->GetBinContent(binmin)
5633 ,projPH->GetBinError(binmin)));
5634 AliInfo(Form("binmin-1 %d BinContent %f BinError %f",binmin-1
5635 ,projPH->GetBinContent(binmin-1)
5636 ,projPH->GetBinError(binmin-1)));
5637 AliInfo(Form("binmin+1 %d BinContent %f BinError %f",binmin+1
5638 ,projPH->GetBinContent(binmin+1)
5639 ,projPH->GetBinError(binmin+1)));
5644 Bool_t case1 = kFALSE;
5645 Bool_t case2 = kFALSE;
5646 Bool_t case4 = kFALSE;
5648 //Determination of min and max
5649 //case binmin <= nbins-3
5651 if((binmin <= (nbins-3)) && ((binmin-2) >= TMath::Min(binmax+4, projPH->GetNbinsX()))){
5652 min = pente->GetBinCenter(binmin-2);
5653 max = pente->GetBinCenter(binmin+2);
5654 x[0] = pente->GetBinCenter(binmin-2);
5655 x[1] = pente->GetBinCenter(binmin-1);
5656 x[2] = pente->GetBinCenter(binmin);
5657 x[3] = pente->GetBinCenter(binmin+1);
5658 x[4] = pente->GetBinCenter(binmin+2);
5659 y[0] = pente->GetBinContent(binmin-2);
5660 y[1] = pente->GetBinContent(binmin-1);
5661 y[2] = pente->GetBinContent(binmin);
5662 y[3] = pente->GetBinContent(binmin+1);
5663 y[4] = pente->GetBinContent(binmin+2);
5664 //Calcul the polynome de Lagrange
5665 c = CalculPolynomeLagrange4(x,y);
5667 if((pente->GetBinContent(binmin+2) <= pente->GetBinContent(binmin+1)) &&
5668 (pente->GetBinContent(binmin-2) <= pente->GetBinContent(binmin-1))) put = kFALSE;
5669 if(((binmin+3) <= (nbins-1)) &&
5670 (pente->GetBinContent(binmin+3) <= pente->GetBinContent(binmin+2)) &&
5671 ((binmin-3) >= TMath::Min(binmax+4, projPH->GetNbinsX())) &&
5672 (pente->GetBinContent(binmin-3) <= pente->GetBinContent(binmin-2))) put = kFALSE;
5673 if((pente->GetBinContent(binmin+2) <= pente->GetBinContent(binmin+1)) &&
5674 (pente->GetBinContent(binmin-2) > pente->GetBinContent(binmin-1))) case1 = kTRUE;
5675 if((pente->GetBinContent(binmin+2) > pente->GetBinContent(binmin+1)) &&
5676 (pente->GetBinContent(binmin-2) <= pente->GetBinContent(binmin-1))) case4 = kTRUE;
5678 //case binmin = nbins-2
5680 if(((binmin == (nbins-2)) && ((binmin-2) >= TMath::Min(binmax+4, projPH->GetNbinsX()))) ||
5682 min = pente->GetBinCenter(binmin-2);
5683 max = pente->GetBinCenter(binmin+1);
5684 x[0] = pente->GetBinCenter(binmin-2);
5685 x[1] = pente->GetBinCenter(binmin-1);
5686 x[2] = pente->GetBinCenter(binmin);
5687 x[3] = pente->GetBinCenter(binmin+1);
5688 y[0] = pente->GetBinContent(binmin-2);
5689 y[1] = pente->GetBinContent(binmin-1);
5690 y[2] = pente->GetBinContent(binmin);
5691 y[3] = pente->GetBinContent(binmin+1);
5692 //Calcul the polynome de Lagrange
5693 c = CalculPolynomeLagrange3(x,y);
5694 //richtung +: nothing
5696 if((pente->GetBinContent(binmin-2) <= pente->GetBinContent(binmin-1))) case2 = kTRUE;
5699 if(((binmin <= (nbins-3)) && ((binmin-1) == TMath::Min(binmax+4, projPH->GetNbinsX()))) ||
5701 min = pente->GetBinCenter(binmin-1);
5702 max = pente->GetBinCenter(binmin+2);
5703 x[0] = pente->GetBinCenter(binmin-1);
5704 x[1] = pente->GetBinCenter(binmin);
5705 x[2] = pente->GetBinCenter(binmin+1);
5706 x[3] = pente->GetBinCenter(binmin+2);
5707 y[0] = pente->GetBinContent(binmin-1);
5708 y[1] = pente->GetBinContent(binmin);
5709 y[2] = pente->GetBinContent(binmin+1);
5710 y[3] = pente->GetBinContent(binmin+2);
5711 //Calcul the polynome de Lagrange
5712 c = CalculPolynomeLagrange3(x,y);
5714 if((pente->GetBinContent(binmin+2) <= pente->GetBinContent(binmin+1))) case2 = kTRUE;
5717 if((binmin <= (nbins-3)) && (binmin == TMath::Min(binmax+4, projPH->GetNbinsX()))){
5718 min = pente->GetBinCenter(binmin);
5719 max = pente->GetBinCenter(binmin+2);
5720 x[0] = pente->GetBinCenter(binmin);
5721 x[1] = pente->GetBinCenter(binmin+1);
5722 x[2] = pente->GetBinCenter(binmin+2);
5723 y[0] = pente->GetBinContent(binmin);
5724 y[1] = pente->GetBinContent(binmin+1);
5725 y[2] = pente->GetBinContent(binmin+2);
5726 //Calcul the polynome de Lagrange
5727 c = CalculPolynomeLagrange2(x,y);
5729 if((pente->GetBinContent(binmin+2) <= pente->GetBinContent(binmin+1))) put = kFALSE;
5732 if(((binmin == (nbins-2)) && ((binmin-1) == TMath::Min(binmax+4, projPH->GetNbinsX()))) ||
5734 min = pente->GetBinCenter(binmin-1);
5735 max = pente->GetBinCenter(binmin+1);
5736 x[0] = pente->GetBinCenter(binmin-1);
5737 x[1] = pente->GetBinCenter(binmin);
5738 x[2] = pente->GetBinCenter(binmin+1);
5739 y[0] = pente->GetBinContent(binmin-1);
5740 y[1] = pente->GetBinContent(binmin);
5741 y[2] = pente->GetBinContent(binmin+1);
5742 //Calcul the polynome de Lagrange
5743 c = CalculPolynomeLagrange2(x,y);
5744 //richtung +: nothing
5745 //richtung -: nothing
5747 //case binmin = nbins-1
5749 if((binmin == (nbins-1)) && ((binmin-2) >= TMath::Min(binmax+4, projPH->GetNbinsX()))){
5750 min = pente->GetBinCenter(binmin-2);
5751 max = pente->GetBinCenter(binmin);
5752 x[0] = pente->GetBinCenter(binmin-2);
5753 x[1] = pente->GetBinCenter(binmin-1);
5754 x[2] = pente->GetBinCenter(binmin);
5755 y[0] = pente->GetBinContent(binmin-2);
5756 y[1] = pente->GetBinContent(binmin-1);
5757 y[2] = pente->GetBinContent(binmin);
5758 //Calcul the polynome de Lagrange
5759 c = CalculPolynomeLagrange2(x,y);
5760 AliInfo("At the limit for the drift!");
5761 //fluctuation too big!
5762 //richtung +: nothing
5764 if((pente->GetBinContent(binmin-2) <= pente->GetBinContent(binmin-1))) put = kFALSE;
5766 if((binmin == (nbins-1)) && ((binmin-2) < TMath::Min(binmax+4, projPH->GetNbinsX()))) {
5768 AliInfo("At the limit for the drift and not usable!");
5772 if((binmin == (nbins-2)) && ((binmin-1) < TMath::Min(binmax+4, projPH->GetNbinsX()))){
5774 AliInfo("For the drift...problem!");
5777 //pass but should not happen
5778 if((binmin <= (nbins-3)) && (binmin < TMath::Min(binmax+4, projPH->GetNbinsX()))){
5780 AliInfo("For the drift...problem!");
5783 if(fDebug == 2) AliInfo(Form("For the drift region binmax %d",binmin));
5786 polynome = new TF1("pol","[0]+[1]*x+[2]*x*x+[3]*x*x*x+[4]*x*x*x*x",min,max);
5787 polynome->SetParameters(c[0],c[1],c[2],c[3],c[4]);
5789 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]));
5791 //AliInfo(Form("GetMinimum of the function %f",polynome->GetMinimumX()));
5792 Double_t step = (max-min)/1000;
5794 Double_t minvalue = 0.0;
5795 Double_t placeminimum = min;
5796 for(Int_t o = 0; o < 1000; o++){
5797 if(o == 0) minvalue = polynome->Eval(l);
5798 if(minvalue > (polynome->Eval(l))){
5799 minvalue = polynome->Eval(l);
5804 fPhd[2] = placeminimum;
5807 if(fDebug == 2) AliInfo(Form("For the drift region extrapolated binmax %f",fPhd[2]));
5809 Float_t fPhdt0 = 0.0;
5810 if(fTakeTheMaxPH) fPhdt0 = fPhd[1];
5811 else fPhdt0 = fPhd[0];
5813 if ((fPhd[2] > fPhd[0]) &&
5814 (fPhd[2] > fPhd[1]) &&
5815 (fPhd[1] > fPhd[0]) &&
5817 fVdriftCoef[3] = (kDrWidth) / (fPhd[2]-fPhd[1]);
5818 if(fFitPHNDB == 3) fNumberFitSuccess++;
5819 if (fPhdt0 >= 0.0) {
5820 fT0Coef[3] = (fPhdt0 - fT0Shift) / widbins;
5821 if (fT0Coef[3] < -1.0) {
5822 fT0Coef[3] = fT0Coef[2];
5826 fT0Coef[3] = fT0Coef[2];
5830 fVdriftCoef[3] = -TMath::Abs(fVdriftCoef[2]);
5831 fT0Coef[3] = fT0Coef[2];
5834 if ((fDebug == 1) ||
5836 fCoefVdrift[3][idect] = fVdriftCoef[3];
5837 fCoefVdriftE[2] [idect] = vdriftCoefE;
5838 fCoefT0[3][idect] = fT0Coef[3];
5839 fCoefT0E[2][idect] = t0CoefE;
5843 TCanvas *cpentei = new TCanvas("cpentei","cpentei",50,50,600,800);
5846 line->SetLineColor(2);
5847 line->DrawLine(fPhd[0],0,fPhd[0],projPH->GetMaximum());
5848 line->DrawLine(fPhd[1],0,fPhd[1],projPH->GetMaximum());
5849 line->DrawLine(fPhd[2],0,fPhd[2],projPH->GetMaximum());
5850 AliInfo(Form("fPhd[0] (beginning of the signal): %f" ,(Float_t) fPhd[0]));
5851 AliInfo(Form("fPhd[1] (end of the amplification region): %f" ,(Float_t) fPhd[1]));
5852 AliInfo(Form("fPhd[2] (end of the drift region): %f" ,(Float_t) fPhd[2]));
5853 AliInfo(Form("fVriftCoef[3] (with only the drift region(default)): %f",(Float_t) fVdriftCoef[3]));
5854 TCanvas *cpentei2 = new TCanvas("cpentei2","cpentei2",50,50,600,800);
5857 TCanvas *cpentei3 = new TCanvas("cpentei3","cpentei3",50,50,600,800);
5870 projPH->SetDirectory(0);
5874 //_____________________________________________________________________________
5875 void AliTRDCalibraFit::FitPH(TH1* projPH, Int_t idect)
5878 // Fit methode for the drift velocity
5882 const Float_t kDrWidth = AliTRDgeometry::DrThick();
5885 TAxis *xpph = projPH->GetXaxis();
5886 Double_t upedge = xpph->GetBinUpEdge(xpph->GetNbins());
5888 TF1 *fPH = new TF1("fPH",AliTRDCalibraFit::PH,-0.05,3.2,6);
5889 fPH->SetParameter(0,0.469); // Scaling
5890 fPH->SetParameter(1,0.18); // Start
5891 fPH->SetParameter(2,0.0857325); // AR
5892 fPH->SetParameter(3,1.89); // DR
5893 fPH->SetParameter(4,0.08); // QA/QD
5894 fPH->SetParameter(5,0.0); // Baseline
5896 TLine *line = new TLine();
5898 fVdriftCoef[0] = 0.0;
5900 Double_t vdriftCoefE = 0.0;
5901 Double_t t0CoefE = 0.0;
5903 if (idect%fFitPHPeriode == 0) {
5905 AliInfo(Form("<AliTRDCalibraFit::FitPH> The detector %d will be fitted",idect));
5906 fPH->SetParameter(0,(projPH->Integral()-(projPH->GetBinContent(1)*projPH->GetNbinsX())) * 0.00028); // Scaling
5907 fPH->SetParameter(1,fPhd[0] - 0.1); // Start
5908 fPH->SetParameter(2,fPhd[1] - fPhd[0]); // AR
5909 fPH->SetParameter(3,fPhd[2] - fPhd[1]); // DR
5910 fPH->SetParameter(4,0.225); // QA/QD
5911 fPH->SetParameter(5,(Float_t) projPH->GetBinContent(1));
5914 projPH->Fit(fPH,"0M","",0.0,upedge);
5918 TCanvas *cpente = new TCanvas("cpente","cpente",50,50,600,800);
5920 projPH->Fit(fPH,"M+","",0.0,upedge);
5922 line->SetLineColor(4);
5923 line->DrawLine(fPH->GetParameter(1)
5925 ,fPH->GetParameter(1)
5926 ,projPH->GetMaximum());
5927 line->DrawLine(fPH->GetParameter(1)+fPH->GetParameter(2)
5929 ,fPH->GetParameter(1)+fPH->GetParameter(2)
5930 ,projPH->GetMaximum());
5931 line->DrawLine(fPH->GetParameter(1)+fPH->GetParameter(2)+fPH->GetParameter(3)
5933 ,fPH->GetParameter(1)+fPH->GetParameter(2)+fPH->GetParameter(3)
5934 ,projPH->GetMaximum());
5937 if (fPH->GetParameter(3) != 0) {
5938 if(fFitPHNDB == 0) fNumberFitSuccess++;
5939 fVdriftCoef[0] = kDrWidth / (fPH->GetParameter(3));
5940 vdriftCoefE = (fPH->GetParError(3)/fPH->GetParameter(3))*fVdriftCoef[0];
5941 fT0Coef[0] = fPH->GetParameter(1);
5942 t0CoefE = fPH->GetParError(1);
5945 fVdriftCoef[0] = -TMath::Abs(fVdriftCoef[2]);
5946 fT0Coef[0] = fT0Coef[2];
5949 if ((fDebug == 1) ||
5951 fCoefVdrift[0][idect] = fVdriftCoef[0];
5952 fCoefVdriftE[0][idect] = vdriftCoefE;
5953 fCoefT0[0][idect] = fT0Coef[0];
5954 fCoefT0E[0][idect] = t0CoefE;
5957 AliInfo(Form("fVdriftCoef[0]: %f",(Float_t) fVdriftCoef[0]));
5964 // Put the default value
5965 if ((fDebug <= 1) ||
5967 fCoefVdrift[0][idect] = -TMath::Abs(fVdriftCoef[2]);
5968 fCoefT0[0][idect] = -TMath::Abs(fT0Coef[2]);
5979 //_____________________________________________________________________________
5980 void AliTRDCalibraFit::FitPRF(TH1 *projPRF, Int_t idect)
5983 // Fit methode for the sigma of the pad response function
5987 Double_t prfCoefE = 0.0;
5991 projPRF->Fit("gaus","0M","",-fRangeFitPRF,fRangeFitPRF);
5996 TCanvas *cfit = new TCanvas("cfit","cfit",50,50,600,800);
5998 projPRF->Fit("gaus","M+","",-fRangeFitPRF,fRangeFitPRF);
6003 fPRFCoef[0] = projPRF->GetFunction("gaus")->GetParameter(2);
6004 prfCoefE = projPRF->GetFunction("gaus")->GetParError(2);
6006 if(fPRFCoef[0] <= 0.0) fPRFCoef[0] = -fPRFCoef[1];
6008 if(fFitPRFNDB == 0) fNumberFitSuccess++;
6011 if ((fDebug == 1) ||
6013 fCoefPRF[0][idect] = fPRFCoef[0];
6014 fCoefPRFE[0][idect] = prfCoefE;
6017 AliInfo(Form("fPRFCoef[0]: %f",(Float_t) fPRFCoef[0]));
6022 //_____________________________________________________________________________
6023 void AliTRDCalibraFit::RmsPRF(TH1 *projPRF, Int_t idect)
6026 // Fit methode for the sigma of the pad response function
6030 Double_t prfCoefE = 0.0;
6034 TCanvas *cfit = new TCanvas("cfit","cfit",50,50,600,800);
6040 fPRFCoef[2] = projPRF->GetRMS();
6042 if(fPRFCoef[2] <= 0.0) fPRFCoef[2] = -fPRFCoef[1];
6044 if(fFitPRFNDB == 2) fNumberFitSuccess++;
6047 if ((fDebug == 1) ||
6049 fCoefPRF[2][idect] = fPRFCoef[2];
6050 fCoefPRFE[1][idect] = prfCoefE;
6053 AliInfo(Form("fPRFCoef[2]: %f",(Float_t) fPRFCoef[2]));
6058 //_____________________________________________________________________________
6059 void AliTRDCalibraFit::FitMean(TH1 *projch, Int_t idect, Double_t nentries)
6062 // Only mean methode for the gain factor
6065 Double_t chargeCoefE1 = 0.0;
6066 if(nentries > 0) chargeCoefE1 = projch->GetRMS()/TMath::Sqrt(nentries);
6069 TCanvas *cpmean = new TCanvas("cpmean","cpmean",50,50,600,800);
6074 if(fFitChargeNDB == 1){
6075 CalculChargeCoefMean(fCountDet[0],(Int_t) idect,kTRUE);
6076 fNumberFitSuccess++;
6078 if ((fDebug == 1) ||
6080 fCoefCharge[1][idect]= fChargeCoef[1];
6081 fCoefChargeE[1][idect]= chargeCoefE1;
6085 //_____________________________________________________________________________
6086 void AliTRDCalibraFit::FitMeanW(TH1 *projch, Int_t idect)
6089 // mean w methode for the gain factor
6093 Int_t nybins = projch->GetNbinsX();
6095 //The weight function
6096 Double_t a = 0.00228515;
6097 Double_t b = -0.00231487;
6098 Double_t c = 0.00044298;
6099 Double_t d = -0.00379239;
6100 Double_t e = 0.00338349;
6108 //A arbitrary error for the moment
6109 Double_t chargeCoefE4 = 0.0;
6110 fChargeCoef[4] = 0.0;
6113 Double_t sumw = 0.0;
6115 Int_t sumAll = (Int_t) projch->GetEntries();
6116 Int_t sumCurrent = 0;
6117 for(Int_t k = 0; k <nybins; k++){
6118 Double_t fraction = Float_t(sumCurrent)/Float_t(sumAll);
6119 if (fraction>0.95) break;
6120 Double_t weight = a + b*fraction + c*fraction*fraction + d *fraction*fraction*fraction+
6121 e*fraction*fraction*fraction*fraction;
6122 sumw += weight*projch->GetBinContent(k+1)*projch->GetBinCenter(k+1);
6123 sum += weight*projch->GetBinContent(k+1);
6124 sumCurrent += (Int_t) projch->GetBinContent(k+1);
6125 //printf("fraction %f, weight %f, bincontent %f\n",fraction,weight,projch->GetBinContent(k+1));
6127 if(sum > 0.0) fChargeCoef[4] = (sumw/sum);
6130 AliInfo(Form("fChargeCoef[4] is %f for the dect %d",fChargeCoef[4],idect));
6131 TCanvas *cpmeanw = new TCanvas("cpmeanw","cpmeanw",50,50,600,800);
6136 if(fFitChargeNDB == 4){
6137 fNumberFitSuccess++;
6138 CalculChargeCoefMean(fCountDet[0],(Int_t) idect,kTRUE);
6140 if ((fDebug == 1) ||
6142 fCoefCharge[4][idect]= fChargeCoef[4];
6143 fCoefChargeE[3][idect]= chargeCoefE4;
6148 //_____________________________________________________________________________
6149 void AliTRDCalibraFit::FitCH(TH1 *projch, Int_t idect)
6152 // Fit methode for the gain factor
6155 fChargeCoef[0] = 0.0;
6156 Double_t chargeCoefE0 = 0.0;
6157 Double_t chisqrl = 0.0;
6158 Double_t chisqrg = 0.0;
6159 Double_t chisqr = 0.0;
6160 TF1 *fLandauGaus = new TF1("fLandauGaus",FuncLandauGaus,0,300,5);
6162 projch->Fit("landau","0",""
6163 ,(Float_t) fChargeCoef[1]/fBeginFitCharge
6164 ,projch->GetBinCenter(projch->GetNbinsX()));
6165 Double_t l3P0 = projch->GetFunction("landau")->GetParameter(0);
6166 Double_t l3P1 = projch->GetFunction("landau")->GetParameter(1);
6167 Double_t l3P2 = projch->GetFunction("landau")->GetParameter(2);
6168 chisqrl = projch->GetFunction("landau")->GetChisquare();
6170 projch->Fit("gaus","0",""
6171 ,(Float_t) fChargeCoef[1]/fBeginFitCharge
6172 ,projch->GetBinCenter(projch->GetNbinsX()));
6173 Double_t g3P0 = projch->GetFunction("gaus")->GetParameter(0);
6174 Double_t g3P2 = projch->GetFunction("gaus")->GetParameter(2);
6175 chisqrg = projch->GetFunction("gaus")->GetChisquare();
6177 fLandauGaus->SetParameters(l3P0,l3P1,l3P2,g3P0,g3P2);
6178 if ((fDebug <= 1) ||
6180 projch->Fit("fLandauGaus","0",""
6181 ,(Float_t) fChargeCoef[1]/fBeginFitCharge
6182 ,projch->GetBinCenter(projch->GetNbinsX()));
6183 chisqr = projch->GetFunction("fLandauGaus")->GetChisquare();
6187 TCanvas *cp = new TCanvas("cp","cp",50,50,600,800);
6189 projch->Fit("fLandauGaus","+",""
6190 ,(Float_t) fChargeCoef[1]/fBeginFitCharge
6191 ,projch->GetBinCenter(projch->GetNbinsX()));
6192 chisqr = projch->GetFunction("fLandauGaus")->GetChisquare();
6194 fLandauGaus->Draw("same");
6197 if ((projch->GetFunction("fLandauGaus")->GetParameter(1) > 0) &&
6198 (projch->GetFunction("fLandauGaus")->GetParError(1) <
6199 (0.05*projch->GetFunction("fLandauGaus")->GetParameter(1))) &&
6200 (chisqr < chisqrl) &&
6201 (chisqr < chisqrg)) {
6202 // Calcul of "real" coef
6203 if(fFitChargeNDB == 0){
6204 fNumberFitSuccess++;
6205 CalculChargeCoefMean(fCountDet[0],(Int_t) idect,kTRUE);
6208 CalculChargeCoefMean(fCountDet[0],(Int_t) idect,kFALSE);
6210 fChargeCoef[0] = projch->GetFunction("fLandauGaus")->GetParameter(1);
6211 chargeCoefE0 = projch->GetFunction("fLandauGaus")->GetParError(1);
6214 // Calcul of "real" coef
6215 CalculChargeCoefMean(fCountDet[0],(Int_t) idect,kFALSE);
6216 fChargeCoef[0] = -TMath::Abs(fChargeCoef[3]);
6220 AliInfo(Form("fChargeCoef[0]: %f",(Float_t) fChargeCoef[0]));
6221 AliInfo(Form("fChargeCoef[1]: %f",(Float_t) fChargeCoef[1]));
6224 if ((fDebug == 1) ||
6226 if (fChargeCoef[0] > 0.0) {
6227 fCoefCharge[0][idect]= fChargeCoef[0];
6228 fCoefChargeE[0][idect]= chargeCoefE0;
6238 //_____________________________________________________________________________
6239 void AliTRDCalibraFit::FitBisCH(TH1* projch, Int_t idect)
6242 // Fit methode for the gain factor more time consuming
6245 //Some parameters to initialise
6246 Double_t widthLandau, widthGaus, MPV, Integral;
6247 Double_t chisquarel = 0.0;
6248 Double_t chisquareg = 0.0;
6250 projch->Fit("landau","0M+",""
6251 ,(Float_t) fChargeCoef[1]/6
6252 ,projch->GetBinCenter(projch->GetNbinsX()));
6253 widthLandau = projch->GetFunction("landau")->GetParameter(2);
6254 chisquarel = projch->GetFunction("landau")->GetChisquare();
6256 projch->Fit("gaus","0M+",""
6257 ,(Float_t) fChargeCoef[1]/6
6258 ,projch->GetBinCenter(projch->GetNbinsX()));
6259 widthGaus = projch->GetFunction("gaus")->GetParameter(2);
6260 chisquareg = projch->GetFunction("gaus")->GetChisquare();
6262 MPV = (projch->GetFunction("landau")->GetParameter(1))/2;
6263 Integral = (projch->GetFunction("gaus")->Integral(0.3*fChargeCoef[1],3*fChargeCoef[1])
6264 + projch->GetFunction("landau")->Integral(0.3*fChargeCoef[1],3*fChargeCoef[1]))/2;
6266 // Setting fit range and start values
6268 //Double_t sv[4] = { l3P2, fChargeCoef[1], projch->Integral("width"), fG3P2 };
6269 //Double_t sv[4] = { fL3P2, fChargeCoef[1], fL3P0, fG3P2 };
6270 Double_t sv[4] = { widthLandau, MPV, Integral, widthGaus};
6271 Double_t pllo[4] = { 0.001, 0.001, projch->Integral()/3, 0.001};
6272 Double_t plhi[4] = { 300.0, 300.0, 30*projch->Integral(), 300.0};
6273 Double_t fp[4] = { 1.0, 1.0, 1.0, 1.0 };
6274 Double_t fpe[4] = { 1.0, 1.0, 1.0, 1.0 };
6275 fr[0] = 0.3 * fChargeCoef[1];
6276 fr[1] = 3.0 * fChargeCoef[1];
6277 fChargeCoef[2] = 0.0;
6278 Double_t chargeCoefE2 = 0.0;
6282 TF1 *fitsnr = LanGauFit(projch,&fr[0],&sv[0]
6287 Double_t projchPeak;
6288 Double_t projchFWHM;
6289 LanGauPro(fp,projchPeak,projchFWHM);
6291 if ((fp[1] > 0) && ((fpe[1] < (0.05*fp[1])) && (chisqr < chisquarel) && (chisqr < chisquareg))) {
6292 //if ((fp[1] > 0) && ((chisqr < chisquarel) && (chisqr < chisquareg))) {
6293 if(fFitChargeNDB == 2){
6294 fNumberFitSuccess++;
6295 CalculChargeCoefMean(fCountDet[0],(Int_t) idect,kTRUE);
6298 CalculChargeCoefMean(fCountDet[0],(Int_t) idect,kFALSE);
6300 fChargeCoef[2] = fp[1];
6301 chargeCoefE2 = fpe[1];
6302 //chargeCoefE2 = chisqr;
6305 CalculChargeCoefMean(fCountDet[0],(Int_t) idect,kFALSE);
6306 fChargeCoef[2] = -TMath::Abs(fChargeCoef[3]);
6310 AliInfo(Form("fChargeCoef[2]: %f",(Float_t) fChargeCoef[2]));
6311 TCanvas *cpy = new TCanvas("cpy","cpy",50,50,600,800);
6314 fitsnr->Draw("same");
6317 if ((fDebug == 1) ||
6319 if (fChargeCoef[2] > 0.0) {
6320 fCoefCharge[2][idect]= fChargeCoef[2];
6321 fCoefChargeE[2][idect]= chargeCoefE2;
6331 //_____________________________________________________________________________
6332 Double_t *AliTRDCalibraFit::CalculPolynomeLagrange2(Double_t *x, Double_t *y)
6335 // Calcul the coefficients of the polynome passant par ces trois points de degre 2
6338 Double_t *c = new Double_t[5];
6339 Double_t x0 = y[0]/((x[0]-x[1])*(x[0]-x[2]));
6340 Double_t x1 = y[1]/((x[1]-x[0])*(x[1]-x[2]));
6341 Double_t x2 = y[2]/((x[2]-x[0])*(x[2]-x[1]));
6346 c[1] = -(x0*(x[1]+x[2])+x1*(x[0]+x[2])+x2*(x[0]+x[1]));
6347 c[0] = x0*x[1]*x[2]+x1*x[0]*x[2]+x2*x[0]*x[1];
6353 //_____________________________________________________________________________
6354 Double_t *AliTRDCalibraFit::CalculPolynomeLagrange3(Double_t *x, Double_t *y)
6357 // Calcul the coefficients of the polynome passant par ces quatre points de degre 3
6360 Double_t *c = new Double_t[5];
6361 Double_t x0 = y[0]/((x[0]-x[1])*(x[0]-x[2])*(x[0]-x[3]));
6362 Double_t x1 = y[1]/((x[1]-x[0])*(x[1]-x[2])*(x[1]-x[3]));
6363 Double_t x2 = y[2]/((x[2]-x[0])*(x[2]-x[1])*(x[2]-x[3]));
6364 Double_t x3 = y[3]/((x[3]-x[0])*(x[3]-x[1])*(x[3]-x[2]));
6368 c[2] = -(x0*(x[1]+x[2]+x[3])
6369 +x1*(x[0]+x[2]+x[3])
6370 +x2*(x[0]+x[1]+x[3])
6371 +x3*(x[0]+x[1]+x[2]));
6372 c[1] = (x0*(x[1]*x[2]+x[1]*x[3]+x[2]*x[3])
6373 +x1*(x[0]*x[2]+x[0]*x[3]+x[2]*x[3])
6374 +x2*(x[0]*x[1]+x[0]*x[3]+x[1]*x[3])
6375 +x3*(x[0]*x[1]+x[0]*x[2]+x[1]*x[2]));
6377 c[0] = -(x0*x[1]*x[2]*x[3]
6380 +x3*x[0]*x[1]*x[2]);
6386 //_____________________________________________________________________________
6387 Double_t *AliTRDCalibraFit::CalculPolynomeLagrange4(Double_t *x, Double_t *y)
6390 // Calcul the coefficients of the polynome passant par ces cinqs points de degre 4
6393 Double_t *c = new Double_t[5];
6394 Double_t x0 = y[0]/((x[0]-x[1])*(x[0]-x[2])*(x[0]-x[3])*(x[0]-x[4]));
6395 Double_t x1 = y[1]/((x[1]-x[0])*(x[1]-x[2])*(x[1]-x[3])*(x[1]-x[4]));
6396 Double_t x2 = y[2]/((x[2]-x[0])*(x[2]-x[1])*(x[2]-x[3])*(x[2]-x[4]));
6397 Double_t x3 = y[3]/((x[3]-x[0])*(x[3]-x[1])*(x[3]-x[2])*(x[3]-x[4]));
6398 Double_t x4 = y[4]/((x[4]-x[0])*(x[4]-x[1])*(x[4]-x[2])*(x[4]-x[3]));
6400 c[4] = x0+x1+x2+x3+x4;
6401 c[3] = -(x0*(x[1]+x[2]+x[3]+x[4])
6402 +x1*(x[0]+x[2]+x[3]+x[4])
6403 +x2*(x[0]+x[1]+x[3]+x[4])
6404 +x3*(x[0]+x[1]+x[2]+x[4])
6405 +x4*(x[0]+x[1]+x[2]+x[3]));
6406 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])
6407 +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])
6408 +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])
6409 +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])
6410 +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]));
6412 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])
6413 +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])
6414 +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])
6415 +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])
6416 +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]));
6418 c[0] = (x0*x[1]*x[2]*x[3]*x[4]
6419 +x1*x[0]*x[2]*x[3]*x[4]
6420 +x2*x[0]*x[1]*x[3]*x[4]
6421 +x3*x[0]*x[1]*x[2]*x[4]
6422 +x4*x[0]*x[1]*x[2]*x[3]);
6428 //_____________________________________________________________________________
6429 void AliTRDCalibraFit::NormierungCharge()
6432 // Normalisation of the gain factor resulting for the fits
6435 // Calcul of the mean of choosen method by fFitChargeNDB
6437 //printf("total number of entries %d\n",fVectorFitCH->GetEntriesFast());
6438 for (Int_t k = 0; k < (Int_t) fVectorFitCH->GetEntriesFast(); k++) {
6440 Int_t detector = ((AliTRDFitCHInfo *) fVectorFitCH->At(k))->GetDetector();
6441 Float_t *coef = ((AliTRDFitCHInfo *) fVectorFitCH->At(k))->GetCoef();
6442 //printf("detector %d coef[0] %f\n",detector,coef[0]);
6443 if (GetChamber(detector) == 2) {
6446 if (GetChamber(detector) != 2) {
6449 for (Int_t j = 0; j < total; j++) {
6457 fScaleFitFactor = fScaleFitFactor / sum;
6460 fScaleFitFactor = 1.0;
6463 if ((fDebug == 3) ||
6465 if ((fFitChargeOn) &&
6466 (fCoefChargeDB[0]->GetEntries() > 0.0) &&
6467 (fCoefChargeDB[0]->GetSumOfWeights() > 0.0)) {
6468 fCoefChargeDB[0]->Scale(fCoefChargeDB[0]->GetEntries() / fCoefChargeDB[0]->GetSumOfWeights());
6470 if ((fMeanChargeOn) &&
6471 (fCoefChargeDB[1]->GetEntries() > 0.0) &&
6472 (fCoefChargeDB[1]->GetSumOfWeights() > 0.0)) {
6473 fCoefChargeDB[1]->Scale(fCoefChargeDB[1]->GetEntries() / fCoefChargeDB[1]->GetSumOfWeights());
6475 if ((fFitChargeBisOn) &&
6476 (fCoefChargeDB[2]->GetEntries() > 0.0) &&
6477 (fCoefChargeDB[2]->GetSumOfWeights() > 0.0)) {
6478 fCoefChargeDB[2]->Scale(fCoefChargeDB[2]->GetEntries() / fCoefChargeDB[2]->GetSumOfWeights());
6480 if ((fFitMeanWOn) &&
6481 (fCoefChargeDB[3]->GetEntries() > 0.0) &&
6482 (fCoefChargeDB[3]->GetSumOfWeights() > 0.0)) {
6483 fCoefChargeDB[3]->Scale(fCoefChargeDB[3]->GetEntries() / fCoefChargeDB[3]->GetSumOfWeights());
6489 //_____________________________________________________________________________
6490 TH1I *AliTRDCalibraFit::ReBin(TH1I *hist) const
6493 // Rebin of the 1D histo for the gain calibration if needed.
6494 // you have to choose fRebin, divider of fNumberBinCharge
6497 TAxis *xhist = hist->GetXaxis();
6498 TH1I *rehist = new TH1I("projrebin","",(Int_t) xhist->GetNbins()/fRebin
6499 ,xhist->GetBinLowEdge(1)
6500 ,xhist->GetBinUpEdge(xhist->GetNbins()));
6502 AliInfo(Form("fRebin: %d",fRebin));
6504 for (Int_t k = 1; k <= (Int_t) xhist->GetNbins()/fRebin; k++) {
6506 for (Int_t ji = i; ji < i+fRebin; ji++) {
6507 sum += hist->GetBinContent(ji);
6510 rehist->SetBinContent(k,sum);
6515 TCanvas *crebin = new TCanvas("crebin","",50,50,600,800);
6524 //_____________________________________________________________________________
6525 TH1F *AliTRDCalibraFit::ReBin(TH1F *hist) const
6528 // Rebin of the 1D histo for the gain calibration if needed
6529 // you have to choose fRebin divider of fNumberBinCharge
6532 TAxis *xhist = hist->GetXaxis();
6533 TH1F *rehist = new TH1F("projrebin","",(Int_t) xhist->GetNbins()/fRebin
6534 ,xhist->GetBinLowEdge(1)
6535 ,xhist->GetBinUpEdge(xhist->GetNbins()));
6537 AliInfo(Form("fRebin: %d",fRebin));
6539 for (Int_t k = 1; k <= (Int_t) xhist->GetNbins()/fRebin; k++) {
6541 for (Int_t ji = i; ji < i+fRebin; ji++) {
6542 sum += hist->GetBinContent(ji);
6545 rehist->SetBinContent(k,sum);
6550 TCanvas *crebin = new TCanvas("crebin","",50,50,600,800);
6559 //_____________________________________________________________________________
6560 TH1F *AliTRDCalibraFit::CorrectTheError(TGraphErrors *hist)
6563 // In the case of the vectors method the trees contains TGraphErrors for PH and PRF
6564 // to be able to add them after
6565 // We convert it to a TH1F to be able to applied the same fit function method
6566 // After having called this function you can not add the statistics anymore
6571 Int_t nbins = hist->GetN();
6572 Double_t *x = hist->GetX();
6573 Double_t *entries = hist->GetEX();
6574 Double_t *mean = hist->GetY();
6575 Double_t *square = hist->GetEY();
6576 fEntriesCurrent = 0;
6582 Double_t step = x[1] - x[0];
6583 Double_t minvalue = x[0] - step/2;
6584 Double_t maxvalue = x[(nbins-1)] + step/2;
6586 rehist = new TH1F("projcorrecterror","",nbins,minvalue,maxvalue);
6588 for (Int_t k = 0; k < nbins; k++) {
6589 rehist->SetBinContent(k+1,mean[k]);
6590 if (entries[k] > 0.0) {
6591 fEntriesCurrent += (Int_t) entries[k];
6592 Double_t d = TMath::Abs(square[k] - (mean[k]*mean[k]));
6593 rehist->SetBinError(k+1,TMath::Sqrt(d/entries[k]));
6596 rehist->SetBinError(k+1,0.0);
6605 //____________Some basic geometry function_____________________________________
6608 //_____________________________________________________________________________
6609 Int_t AliTRDCalibraFit::GetPlane(Int_t d) const
6612 // Reconstruct the plane number from the detector number
6615 return ((Int_t) (d % 6));
6619 //_____________________________________________________________________________
6620 Int_t AliTRDCalibraFit::GetChamber(Int_t d) const
6623 // Reconstruct the chamber number from the detector number
6627 return ((Int_t) (d % 30) / fgkNplan);
6631 //_____________________________________________________________________________
6632 Int_t AliTRDCalibraFit::GetSector(Int_t d) const
6635 // Reconstruct the sector number from the detector number
6639 return ((Int_t) (d / fg));
6644 //____________Fill and Init tree Gain, PRF, Vdrift and T0______________________
6647 //_____________________________________________________________________________
6648 void AliTRDCalibraFit::InitTreePRF()
6651 // Init the tree where the coefficients from the fit methods can be stored
6655 fPRFPad = new Float_t[2304];
6656 fPRF = new TTree("PRF","PRF");
6657 fPRF->Branch("detector",&fPRFDetector,"detector/I");
6658 fPRF->Branch("width" ,fPRFPad ,"width[2304]/F");
6660 // Set to default value for the plane 0 supposed to be the first one
6661 for (Int_t k = 0; k < 2304; k++) {
6668 //_____________________________________________________________________________
6669 void AliTRDCalibraFit::FillTreePRF(Int_t countdet)
6672 // Fill the tree with the sigma of the pad response function for the detector countdet
6675 Int_t numberofgroup = 0;
6676 fPRFDetector = countdet;
6679 if (GetChamber((Int_t)(countdet+1)) == 2) {
6680 numberofgroup = 1728;
6683 numberofgroup = 2304;
6686 // Reset to default value for the next
6687 for (Int_t k = 0; k < numberofgroup; k++) {
6688 if (GetPlane((Int_t) (countdet+1)) == 0) {
6691 if (GetPlane((Int_t) (countdet+1)) == 1) {
6694 if (GetPlane((Int_t) (countdet+1)) == 2) {
6697 if (GetPlane((Int_t) (countdet+1)) == 3) {
6700 if (GetPlane((Int_t) (countdet+1)) == 4) {
6703 if (GetPlane((Int_t) (countdet+1)) == 5) {
6712 //_____________________________________________________________________________
6713 void AliTRDCalibraFit::ConvertVectorFitCHTree()
6716 // Convert the vector stuff to a tree of 1D histos if the user
6717 // want to write it after the fill functions
6720 Int_t detector = -1;
6721 Int_t numberofgroup = 1;
6722 Float_t gainPad[2304];
6724 fGain = new TTree("Gain","Gain");
6725 fGain->Branch("detector",&detector,"detector/I");
6726 fGain->Branch("gainPad" ,gainPad ,"gainPad[2304]/F");
6728 Int_t loop = (Int_t) fVectorFitCH->GetEntriesFast();
6729 for (Int_t k = 0; k < loop; k++) {
6730 detector = ((AliTRDFitCHInfo *) fVectorFitCH->At(k))->GetDetector();
6731 if (GetChamber((Int_t) ((AliTRDFitCHInfo *) fVectorFitCH->At(k))->GetDetector()) == 2) {
6732 numberofgroup = 1728;
6735 numberofgroup = 2304;
6737 for (Int_t i = 0; i < numberofgroup; i++) {
6738 if (((AliTRDFitCHInfo *) fVectorFitCH->At(k))->GetCoef()[i] >= 0) {
6739 gainPad[i] = ((AliTRDFitCHInfo *) fVectorFitCH->At(k))->GetCoef()[i] * fScaleFitFactor;
6742 gainPad[i] = (Float_t) ((AliTRDFitCHInfo *) fVectorFitCH->At(k))->GetCoef()[i];
6746 if(numberofgroup < 2304){
6747 for(Int_t i = numberofgroup; i < 2304; i++){
6748 gainPad[i] = -100.0;
6756 //_____________________________________________________________________________
6757 void AliTRDCalibraFit::FillTreeVdrift(Int_t countdet)
6760 // Fill the tree with the drift velocities for the detector countdet
6763 Int_t numberofgroup = 0;
6764 fVdriftDetector = countdet;
6767 if (GetChamber((Int_t)(countdet+1)) == 2) {
6768 numberofgroup = 1728;
6771 numberofgroup = 2304;
6773 // Reset to default value the gain coef
6774 for (Int_t k = 0; k < numberofgroup; k++) {
6775 fVdriftPad[k] = -1.5;
6777 fVdriftDetector = -1;
6781 //_____________________________________________________________________________
6782 void AliTRDCalibraFit::InitTreePH()
6785 // Init the tree where the coefficients from the fit methods can be stored
6789 fVdriftPad = new Float_t[2304];
6790 fVdrift = new TTree("Vdrift","Vdrift");
6791 fVdrift->Branch("detector",&fVdriftDetector,"detector/I");
6792 fVdrift->Branch("vdrift" ,fVdriftPad ,"vdrift[2304]/F");
6793 // Set to default value for the plane 0 supposed to be the first one
6794 for (Int_t k = 0; k < 2304; k++) {
6795 fVdriftPad[k] = -1.5;
6797 fVdriftDetector = -1;
6801 //_____________________________________________________________________________
6802 void AliTRDCalibraFit::FillTreeT0(Int_t countdet)
6805 // Fill the tree with the t0 value for the detector countdet
6808 Int_t numberofgroup = 0;
6810 fT0Detector = countdet;
6813 if (GetChamber((Int_t) (countdet+1)) == 2) {
6814 numberofgroup = 1728;
6817 numberofgroup = 2304;
6819 // Reset to default value
6820 for (Int_t k = 0; k < numberofgroup; k++) {
6827 //_____________________________________________________________________________
6828 void AliTRDCalibraFit::InitTreeT0()
6831 // Init the tree where the coefficients from the fit methods can be stored
6835 fT0Pad = new Float_t[2304];
6836 fT0 = new TTree("T0","T0");
6837 fT0->Branch("detector",&fT0Detector,"detector/I");
6838 fT0->Branch("t0",fT0Pad,"t0[2304]/F");
6839 //Set to default value for the plane 0 supposed to be the first one
6840 for(Int_t k = 0; k < 2304; k++){
6848 //____________Private Functions________________________________________________
6851 //_____________________________________________________________________________
6852 Double_t AliTRDCalibraFit::PH(Double_t *x, Double_t *par)
6855 // Function for the fit
6858 //TF1 *fAsymmGauss = new TF1("fAsymmGauss",AsymmGauss,0,4,6);
6860 //PARAMETERS FOR FIT PH
6862 //fAsymmGauss->SetParameter(0,0.113755);
6863 //fAsymmGauss->SetParameter(1,0.350706);
6864 //fAsymmGauss->SetParameter(2,0.0604244);
6865 //fAsymmGauss->SetParameter(3,7.65596);
6866 //fAsymmGauss->SetParameter(4,1.00124);
6867 //fAsymmGauss->SetParameter(5,0.870597); // No tail cancelation
6875 Double_t dx = 0.005;
6876 Double_t xs = par[1];
6878 Double_t paras[2] = { 0.0, 0.0 };
6881 if ((xs >= par[1]) &&
6882 (xs < (par[1]+par[2]))) {
6883 //fAsymmGauss->SetParameter(0,par[0]);
6884 //fAsymmGauss->SetParameter(1,xs);
6885 //ss += fAsymmGauss->Eval(xx);
6888 ss += AsymmGauss(&xx,paras);
6890 if ((xs >= (par[1]+par[2])) &&
6891 (xs < (par[1]+par[2]+par[3]))) {
6892 //fAsymmGauss->SetParameter(0,par[0]*par[4]);
6893 //fAsymmGauss->SetParameter(1,xs);
6894 //ss += fAsymmGauss->Eval(xx);
6895 paras[0] = par[0]*par[4];
6897 ss += AsymmGauss(&xx,paras);
6906 //_____________________________________________________________________________
6907 Double_t AliTRDCalibraFit::AsymmGauss(Double_t *x, Double_t *par)
6910 // Function for the fit
6913 //par[0] = normalization
6921 Double_t par1save = par[1];
6922 //Double_t par2save = par[2];
6923 Double_t par2save = 0.0604244;
6924 //Double_t par3save = par[3];
6925 Double_t par3save = 7.65596;
6926 //Double_t par5save = par[5];
6927 Double_t par5save = 0.870597;
6928 Double_t dx = x[0] - par1save;
6930 Double_t sigma2 = par2save*par2save;
6931 Double_t sqrt2 = TMath::Sqrt(2.0);
6932 Double_t exp1 = par3save * TMath::Exp(-par3save * (dx - 0.5 * par3save * sigma2))
6933 * (1.0 - TMath::Erf((par3save * sigma2 - dx) / (sqrt2 * par2save)));
6934 Double_t exp2 = par5save * TMath::Exp(-par5save * (dx - 0.5 * par5save * sigma2))
6935 * (1.0 - TMath::Erf((par5save * sigma2 - dx) / (sqrt2 * par2save)));
6937 //return par[0]*(exp1+par[4]*exp2);
6938 return par[0] * (exp1 + 1.00124 * exp2);
6942 //_____________________________________________________________________________
6943 Double_t AliTRDCalibraFit::FuncLandauGaus(Double_t *x, Double_t *par)
6946 // Sum Landau + Gaus with identical mean
6949 Double_t valLandau = par[0] * TMath::Landau(x[0],par[1],par[2]);
6950 //Double_t valGaus = par[3] * TMath::Gaus(x[0],par[4],par[5]);
6951 Double_t valGaus = par[3] * TMath::Gaus(x[0],par[1],par[4]);
6952 Double_t val = valLandau + valGaus;
6958 //_____________________________________________________________________________
6959 Double_t AliTRDCalibraFit::LanGauFun(Double_t *x, Double_t *par)
6962 // Function for the fit
6965 // par[0]=Width (scale) parameter of Landau density
6966 // par[1]=Most Probable (MP, location) parameter of Landau density
6967 // par[2]=Total area (integral -inf to inf, normalization constant)
6968 // par[3]=Width (sigma) of convoluted Gaussian function
6970 // In the Landau distribution (represented by the CERNLIB approximation),
6971 // the maximum is located at x=-0.22278298 with the location parameter=0.
6972 // This shift is corrected within this function, so that the actual
6973 // maximum is identical to the MP parameter.
6976 // Numeric constants
6977 Double_t invsq2pi = 0.3989422804014; // (2 pi)^(-1/2)
6978 Double_t mpshift = -0.22278298; // Landau maximum location
6980 // Control constants
6981 Double_t np = 100.0; // Number of convolution steps
6982 Double_t sc = 5.0; // Convolution extends to +-sc Gaussian sigmas
6994 // MP shift correction
6995 mpc = par[1] - mpshift * par[0];
6997 // Range of convolution integral
6998 xlow = x[0] - sc * par[3];
6999 xupp = x[0] + sc * par[3];
7001 step = (xupp - xlow) / np;
7003 // Convolution integral of Landau and Gaussian by sum
7004 for (i = 1.0; i <= np/2; i++) {
7006 xx = xlow + (i-.5) * step;
7007 fland = TMath::Landau(xx,mpc,par[0]) / par[0];
7008 sum += fland * TMath::Gaus(x[0],xx,par[3]);
7010 xx = xupp - (i-.5) * step;
7011 fland = TMath::Landau(xx,mpc,par[0]) / par[0];
7012 sum += fland * TMath::Gaus(x[0],xx,par[3]);
7016 return (par[2] * step * sum * invsq2pi / par[3]);
7020 //_____________________________________________________________________________
7021 TF1 *AliTRDCalibraFit::LanGauFit(TH1 *his, Double_t *fitrange, Double_t *startvalues
7022 , Double_t *parlimitslo, Double_t *parlimitshi
7023 , Double_t *fitparams, Double_t *fiterrors
7024 , Double_t *chiSqr, Int_t *ndf)
7027 // Function for the fit
7031 Char_t funname[100];
7033 TF1 *ffitold = (TF1 *) gROOT->GetListOfFunctions()->FindObject(funname);
7038 TF1 *ffit = new TF1(funname,LanGauFun,fitrange[0],fitrange[1],4);
7039 ffit->SetParameters(startvalues);
7040 ffit->SetParNames("Width","MP","Area","GSigma");
7042 for (i = 0; i < 4; i++) {
7043 ffit->SetParLimits(i,parlimitslo[i],parlimitshi[i]);
7046 his->Fit(funname,"RB0"); // Fit within specified range, use ParLimits, do not plot
7048 ffit->GetParameters(fitparams); // Obtain fit parameters
7049 for (i = 0; i < 4; i++) {
7050 fiterrors[i] = ffit->GetParError(i); // Obtain fit parameter errors
7052 chiSqr[0] = ffit->GetChisquare(); // Obtain chi^2
7053 ndf[0] = ffit->GetNDF(); // Obtain ndf
7055 return (ffit); // Return fit function
7059 //_____________________________________________________________________________
7060 Int_t AliTRDCalibraFit::LanGauPro(Double_t *params, Double_t &maxx, Double_t &fwhm)
7063 // Function for the fit
7076 Int_t maxcalls = 10000;
7078 // Search for maximum
7079 p = params[1] - 0.1 * params[0];
7080 step = 0.05 * params[0];
7084 while ((l != lold) && (i < maxcalls)) {
7088 l = LanGauFun(&x,params);
7090 step = -step / 10.0;
7095 if (i == maxcalls) {
7101 // Search for right x location of fy
7102 p = maxx + params[0];
7108 while ( (l != lold) && (i < maxcalls) ) {
7113 l = TMath::Abs(LanGauFun(&x,params) - fy);
7126 // Search for left x location of fy
7128 p = maxx - 0.5 * params[0];
7134 while ((l != lold) && (i < maxcalls)) {
7138 l = TMath::Abs(LanGauFun(&x,params) - fy);
7140 step = -step / 10.0;
7145 if (i == maxcalls) {
7156 //_____________________________________________________________________________
7157 Double_t AliTRDCalibraFit::GausConstant(Double_t *x, Double_t *par)
7160 // Gaus with identical mean
7163 Double_t gauss = par[0] * TMath::Gaus(x[0],0.0,par[1])+par[2];