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.
24 // It can be used for the calibration per chamber but also per group of pads and eventually per pad.
25 // The user has to choose with the functions SetNz and SetNrphi the precision of the calibration.
31 <TR><TD><center>Nz</center></TD><TD><center> 0 </center></TD><TD><center> 1 </center></TD><TD><center> 2 </center></TD><TD><center> 3 </center></TD><TD><center> 4 </center></TD></TR>
32 <TR><TD><CENTER>group of row pads per detector</CENTER></TD><TD><CENTER>1</CENTER></TD><TD><CENTER>2</CENTER></TD><TD><CENTER>4</CENTER></TD><TD><CENTER>6(chamb2)<br> 8(others chambers)</CENTER></TD><TD><CENTER>12 (chamb2)<br> 16 (chamb0)</CENTER></TD></TR>
33 <TR><TD><CENTER>row pads per group</CENTER></TD><TD><CENTER>12 (chamb2)<br> 16 (chamb0)</CENTER></TD><TD><CENTER>6 (chamb2)<br> 8 (chamb0)</CENTER></TD><TD><CENTER>3 (chamb2)<br> 4 (chamb0)</CENTER></TD><TD><CENTER>2</CENTER></TD><TD><CENTER>1</CENTER></TD></TR>
34 <TR><TD><CENTER>~distance [cm]</CENTER></TD><TD><CENTER>106 (chamb2)<br> 130 (chamb0)</CENTER></TD><TD><CENTER>53 (chamb2)<br> 65 (chamb0)</CENTER></TD><TD><CENTER>26.5 (chamb2)<br> 32.5 (chamb0)</CENTER></TD><TD><CENTER>17 (chamb2)<br> 17 (chamb0)</CENTER></TD><TD><CENTER>9 (chamb2)<br> 9 (chamb0)</CENTER></TD></TR>
35 <CAPTION>In the z direction</CAPTION>
41 <TR><TD><center>Nrphi</center></TD><TD><center> 0 </center></TD><TD><center> 1 </center></TD><TD><center> 2 </center></TD><TD><center> 3 </center></TD><TD><center> 4 </center></TD><TD><center> 5 </center></TD><TD><center> 6 </center></TD></TR>
42 <TR><TD><CENTER>group of col pads per detector</CENTER></TD><TD><CENTER>1</CENTER></TD><TD><CENTER>2</CENTER></TD><TD><CENTER>4</CENTER></TD><TD><CENTER>8</CENTER></TD><TD><CENTER>16</CENTER></TD><TD><center>36</center></TD><TD><center>144</center></TD></TR>
43 <TR><TD><CENTER>col pads per group</CENTER></TD><TD><CENTER>144</CENTER></TD><TD><CENTER>72</CENTER></TD><TD><CENTER>36</CENTER></TD><TD><CENTER>18</CENTER></TD><TD><CENTER>9</CENTER></TD><TD><center>4</center></TD><TD><center>1</center></TD></TR>
44 <TR><TD><CENTER>~distance [cm]</CENTER></TD><TD><CENTER>113.4</CENTER></TD><TD><CENTER>56.7</CENTER></TD><TD><CENTER>25.3</CENTER></TD><TD><CENTER>14.3</CENTER></TD><TD><CENTER>7.25</CENTER></TD><TD><center>3.2</center></TD><TD><center>0.8</center></TD></TR>
45 <CAPTION>In the rphi direction</CAPTION>
52 // Fill histograms or vectors
53 //----------------------------
55 // 2D Histograms (Histo2d) or vectors (Vector2d), then converted in Trees, will be filled
56 // from RAW DATA in a run or from reconstructed TRD tracks during the offline tracking
57 // in the function "FollowBackProlongation" (AliTRDtracker)
58 // Per default the functions to fill are off.
62 Example of 2D histos for the relative gain (left) and the drift velocity (right) calibration of the sector 13 <br>
64 <img src="./gif/2dhisto.gif" width="600" height="350"><br>
69 // Fit the histograms to find the coefficients
70 //---------------------------------------------
72 // These 2D histograms or vectors (first converted in 1D histos) will be fitted
73 // if they have enough entries, otherwise the (default) value of the choosen database
74 // will be put. For the relative gain calibration the resulted factors will be globally
75 // normalized to the gain factors of the choosen database. It unables to precise
76 // previous calibration procedure.
77 // The function SetDebug enables the user to see:
78 // _fDebug = 0: nothing, only the values are written in the tree if wanted
79 // _fDebug = 1: a comparaison of the coefficients found and the default values
80 // in the choosen database.
81 // fCoef , histogram of the coefs as function of the calibration group number
82 // fDelta , histogram of the relative difference of the coef with the default
83 // value in the database as function of the calibration group number
84 // fError , dirstribution of this relative difference
85 // _fDebug = 2: only the fit of the choosen calibration group fFitVoir (SetFitVoir)
86 // _fDebug = 3: The coefficients in the choosen detector fDet (SetDet) as function of the
87 // pad row and col number
88 // _fDebug = 4; The coeffcicients in the choosen detector fDet (SetDet) like in the 3 but with
89 // also the comparaison histograms of the 1 for this detector
93 Example of fCoef for the relative gain calibration of the sector 13 <br>
95 <img src="./gif/coef.gif" width="400" height="460">
97 Example of fDelta (right) and fError (left) for the relative gain calibration of the sector 13 <br>
99 <img src="./gif/detlaerror.gif" width="550" height="380"><br>
105 // R. Bailhache (R.Bailhache@gsi.de)
107 //////////////////////////////////////////////////////////////////////////////////////
113 #include <TProfile2D.h>
116 #include <TGraphErrors.h>
118 #include <TObjArray.h>
126 #include <TStopwatch.h>
129 #include <TDirectory.h>
133 #include "AliCDBManager.h"
135 #include "AliRunLoader.h"
136 #include "AliLoader.h"
137 #include "AliRawReaderFile.h"
138 #include "AliRawReader.h"
140 #include "AliTRDCalibra.h"
141 #include "AliTRDcalibDB.h"
142 #include "AliTRDCommonParam.h"
143 #include "AliTRDmcmTracklet.h"
144 #include "AliTRDpadPlane.h"
145 #include "AliTRDcluster.h"
146 #include "AliTRDtrack.h"
147 #include "AliTRDdigit.h"
148 #include "AliTRDdigitsManager.h"
150 #include "AliTRDgeometry.h"
151 #include "./Cal/AliTRDCalROC.h"
152 #include "./Cal/AliTRDCalPad.h"
153 #include "./Cal/AliTRDCalDet.h"
154 #include "AliTRDrawData.h"
156 ClassImp(AliTRDCalibra)
158 AliTRDCalibra* AliTRDCalibra::fgInstance = 0;
159 Bool_t AliTRDCalibra::fgTerminated = kFALSE;
161 //_____________singleton implementation_________________________________________________
162 AliTRDCalibra *AliTRDCalibra::Instance()
165 // Singleton implementation
168 if (fgTerminated != kFALSE) {
172 if (fgInstance == 0) {
173 fgInstance = new AliTRDCalibra();
180 //______________________________________________________________________________________
181 void AliTRDCalibra::Terminate()
184 // Singleton implementation
185 // Deletes the instance of this class
188 fgTerminated = kTRUE;
190 if (fgInstance != 0) {
197 //______________________________________________________________________________________
198 AliTRDCalibra::AliTRDCalibra()
201 ,fMcmTracking(kFALSE)
202 ,fMcmCorrectAngle(kFALSE)
209 ,fCountRelativeScale(0)
210 ,fRelativeScaleAuto(kFALSE)
212 ,fThresholdClusterPRF1(0.0)
213 ,fThresholdClusterPRF2(0.0)
214 ,fCenterOfflineCluster(kFALSE)
220 ,fBeginFitCharge(0.0)
222 ,fMeanChargeOn(kFALSE)
223 ,fFitChargeBisOn(kFALSE)
242 ,fDetectorAliTRDtrack(kFALSE)
243 ,fChamberAliTRDtrack(-1)
244 ,fDetectorPreviousTrack(-1)
257 ,fScaleFitFactor(0.0)
278 // Default constructor
281 for (Int_t i = 0; i < 3; i++) {
286 for (Int_t k = 0; k < 3; k++) {
293 for (Int_t i = 0; i < 3; i++) {
294 fWriteCoef[i] = kFALSE;
299 for (Int_t k = 0; k < 3; k++) {
303 for (Int_t i = 0; i < 3; i++) {
312 //______________________________________________________________________________________
313 AliTRDCalibra::AliTRDCalibra(const AliTRDCalibra &c)
316 ,fMcmTracking(kFALSE)
317 ,fMcmCorrectAngle(kFALSE)
324 ,fCountRelativeScale(0)
325 ,fRelativeScaleAuto(kFALSE)
327 ,fThresholdClusterPRF1(0.0)
328 ,fThresholdClusterPRF2(0.0)
329 ,fCenterOfflineCluster(kFALSE)
335 ,fBeginFitCharge(0.0)
337 ,fMeanChargeOn(kFALSE)
338 ,fFitChargeBisOn(kFALSE)
357 ,fDetectorAliTRDtrack(kFALSE)
358 ,fChamberAliTRDtrack(-1)
359 ,fDetectorPreviousTrack(-1)
372 ,fScaleFitFactor(0.0)
398 //____________________________________________________________________________________
399 AliTRDCalibra::~AliTRDCalibra()
402 // AliTRDCalibra destructor
410 //_____________________________________________________________________________
411 void AliTRDCalibra::Destroy()
424 //_____________________________________________________________________________
425 void AliTRDCalibra::ClearHistos()
446 //_____________________________________________________________________________
447 void AliTRDCalibra::ClearTree()
472 //_____________________________________________________________________________
473 void AliTRDCalibra::Init()
476 // Init some default values
479 // How to fill the 2D
481 fThresholdClusterPRF1 = 2.0;
482 fThresholdClusterPRF2 = 3.0;
485 fNumberBinCharge = 100;
489 fWriteName = "TRD.calibration.root";
490 fWriteNameCoef = "TRD.coefficient.root";
494 fBeginFitCharge = 3.5;
499 // Internal variables
501 // Fill the 2D histos in the offline tracking
502 fDetectorPreviousTrack = -1;
503 fChamberAliTRDtrack = -1;
508 fNumberClusters = 18;
510 fNumberUsedCh[0] = 0;
511 fNumberUsedCh[1] = 0;
512 fNumberUsedPh[0] = 0;
513 fNumberUsedPh[1] = 0;
515 // Variables in the loop
516 for (Int_t k = 0; k < 4; k++) {
517 fChargeCoef[k] = 1.0;
518 fVdriftCoef[k] = 1.5;
521 for (Int_t i = 0; i < 2; i++) {
526 for (Int_t i = 0; i < 3; i++) {
538 // Local database to be changed
543 //____________Functions fit Online CH2d________________________________________
544 Bool_t AliTRDCalibra::FitCHOnline(TH2I *ch)
547 // Fit the 1D histos, projections of the 2D ch on the Xaxis, for each
548 // calibration group normalized the resulted coefficients (to 1 normally)
549 // and write the results in a tree
552 // Number of Xbins (detectors or groups of pads)
553 TAxis *xch = ch->GetXaxis();
554 Int_t nbins = xch->GetNbins();
555 TAxis *yph = ch->GetYaxis();
556 Int_t nybins = yph->GetNbins();
557 if (!InitFit(nbins,0)) {
560 fStatisticMean = 0.0;
572 // Init fCountDet and fCount
573 InitfCountDetAndfCount(0);
575 // Beginning of the loop betwwen dect1 and dect2
576 for (Int_t idect = fDect1[0]; idect < fDect2[0]; idect++) {
578 TH1I *projch = (TH1I *) ch->ProjectionY("projch",idect+1,idect+1,(Option_t *)"e");
579 projch->SetDirectory(0);
581 // Determination of fNnZ, fNnRphi, fNfragZ and fNfragRphi
582 UpdatefCountDetAndfCount(idect,0);
584 // Reconstruction of the row and pad group: rowmin, row max ...
585 ReconstructFitRowMinRowMax(idect, 0);
587 // Number of entries for this calibration group
588 Double_t nentries = 0.0;
589 for (Int_t k = 0; k < nybins; k++) {
590 nentries += ch->GetBinContent(ch->GetBin(idect+1,k+1));
596 // Rebin and statistic stuff
599 projch = ReBin((TH1I *) projch);
601 // This detector has not enough statistics or was off
602 if (nentries < fMinEntries) {
603 // Fill with the default infos
604 NotEnoughStatistic(idect,0);
612 // Statistics of the group fitted
613 AliInfo(Form("For the group number %d there are %f stats",idect,nentries));
614 fStatisticMean += nentries;
617 // Method Mean and fit
618 // idect is egal for fDebug = 0 and 2, only to fill the hist
619 FitCH((TH1 *) projch,(Int_t) (idect-fDect1[0]));
621 // idect is egal for fDebug = 0 and 2, only to fill the hist
622 if (fFitChargeBisOn) {
623 FitBisCH((TH1 *) projch,(Int_t) (idect-fDect1[0]));
626 // Visualise the detector for fDebug 3 or 4
627 // Here is the reconstruction of the pad and row group is used!
632 FillInfosFit(idect,0);
647 // 0 no plot, 1 and 4 error plot, 3 and 4 DB plot
658 if (fNumberFit > 0) {
659 AliInfo(Form("There is a mean statistic of: %d",(Int_t) fStatisticMean/fNumberFit));
660 fStatisticMean = fStatisticMean / fNumberFit;
663 AliInfo("There is no fit!");
667 ConvertVectorFitCHTree();
676 //____________Functions fit Online CH2d________________________________________
677 Bool_t AliTRDCalibra::FitCHOnline()
680 // Reconstruct a 1D histo from the vectorCH for each calibration group,
681 // fit the histo, normalized the resulted coefficients (to 1 normally)
682 // and write the results in a tree
685 // Number of Xbins (detectors or groups of pads)
689 fStatisticMean = 0.0;
693 // Init fCountDet and fCount
694 InitfCountDetAndfCount(0);
696 // Beginning of the loop between dect1 and dect2
697 for (Int_t idect = fDect1[0]; idect < fDect2[0]; idect++) {
699 // Search if the group is in the VectorCH
700 Int_t place = SearchInVector(idect,0);
708 AliTRDCTInfo *fCHInfo = new AliTRDCTInfo();
710 fCHInfo = ((AliTRDCTInfo *) fVectorCH->At(place));
711 projch = ConvertVectorCTHisto(fCHInfo,(const char *) name);
712 projch->SetDirectory(0);
716 // Determination of fNnZ, fNnRphi, fNfragZ and fNfragRphi
717 UpdatefCountDetAndfCount(idect,0);
719 // Reconstruction of the row and pad group: rowmin, row max ...
720 ReconstructFitRowMinRowMax(idect,0);
723 Double_t nentries = 0.0;
725 for (Int_t k = 0; k < fNumberBinCharge; k++) {
726 nentries += projch->GetBinContent(k+1);
733 // Rebin and statistic stuff
737 projch = ReBin((TH1F *) projch);
740 // This detector has not enough statistics or was not found in VectorCH
743 (nentries < fMinEntries))) {
745 // Fill with the default infos
746 NotEnoughStatistic(idect,0);
757 // Statistic of the histos fitted
758 AliInfo(Form("For the group number %d there are %f stats",idect,nentries));
760 fStatisticMean += nentries;
762 // Method Mean and fit
763 // idect is egal for fDebug = 0 and 2, only to fill the hist
764 FitCH((TH1 *) projch,(Int_t) (idect-fDect1[0]));
767 // idect is egal for fDebug = 0 and 2, only to fill the hist
768 if (fFitChargeBisOn) {
769 FitBisCH((TH1 *) projch,(Int_t) (idect-fDect1[0]));
772 // Visualise the detector for fDebug 3 or 4
773 // Here is the reconstruction of the pad and row group is used!
779 FillInfosFit(idect,0);
795 // 0 no plot, 1 and 4 error plot, 3 and 4 DB plot
806 if (fNumberFit > 0) {
807 AliInfo(Form("There is a mean statistic of: %d",(Int_t) fStatisticMean / fNumberFit));
808 fStatisticMean = fStatisticMean / fNumberFit;
811 AliInfo("There is no fit!");
815 ConvertVectorFitCHTree();
824 //____________Functions fit Online CH2d________________________________________
825 Bool_t AliTRDCalibra::FitCHOnline(TTree *tree)
828 // Look if the calibration group can be found in the tree, if yes take the
829 // histo, fit it, normalized the resulted coefficients (to 1 normally) and
830 // write the results in a tree
833 // Number of Xbins (detectors or groups of pads)
837 fStatisticMean = 0.0;
849 // Init fCountDet and fCount
850 InitfCountDetAndfCount(0);
852 tree->SetBranchAddress("histo",&projch);
853 TObjArray *vectorplace = ConvertTreeVector(tree);
855 // Beginning of the loop between dect1 and dect2
856 for (Int_t idect = fDect1[0]; idect < fDect2[0]; idect++) {
858 //Search if the group is in the VectorCH
859 Int_t place = SearchInTreeVector(vectorplace,idect);
864 tree->GetEntry(place);
867 // Determination of fNnZ, fNnRphi, fNfragZ and fNfragRphi
868 UpdatefCountDetAndfCount(idect,0);
870 // Reconstruction of the row and pad group: rowmin, row max ...
871 ReconstructFitRowMinRowMax(idect,0);
874 Double_t nentries = 0.0;
876 for (Int_t k = 0; k < projch->GetXaxis()->GetNbins(); k++) {
877 nentries += projch->GetBinContent(k+1);
884 // Rebin and statistic stuff
888 projch = ReBin((TH1F *) projch);
891 // This detector has not enough statistics or was not found in VectorCH
894 (nentries < fMinEntries))) {
896 // Fill with the default infos
897 NotEnoughStatistic(idect,0);
903 // Statistics of the group fitted
904 AliInfo(Form("For the group number %d there are %f stats",idect,nentries));
906 fStatisticMean += nentries;
908 // Method Mean and fit
909 // idect is egal for fDebug = 0 and 2, only to fill the hist
910 FitCH((TH1 *) projch,(Int_t) (idect-fDect1[0]));
913 // idect is egal for fDebug = 0 and 2, only to fill the hist
914 if (fFitChargeBisOn) {
915 FitBisCH((TH1 *) projch,(Int_t) (idect-fDect1[0]));
918 // Visualise the detector for fDebug 3 or 4
919 // Here is the reconstruction of the pad and row group is used!
925 FillInfosFit(idect,0);
936 // 0 no plot, 1 and 4 error plot, 3 and 4 DB plot
947 if (fNumberFit > 0) {
948 AliInfo(Form("There is a mean statistic of: %d",(Int_t) fStatisticMean / fNumberFit));
949 fStatisticMean = fStatisticMean / fNumberFit;
952 AliInfo("There is no fit!");
956 ConvertVectorFitCHTree();
965 //________________functions fit Online PH2d____________________________________
966 Bool_t AliTRDCalibra::FitPHOnline(TProfile2D *ph)
969 // Take the 1D profiles (average pulse height), projections of the 2D PH
970 // on the Xaxis, for each calibration group
971 // Fit or use the slope of the average pulse height to reconstruct the
972 // drift velocity write the results in a tree
973 // A first calibration of T0 is also made using the same method (slope method)
976 // Number of Xbins (detectors or groups of pads)
977 TAxis *xph = ph->GetXaxis();
978 TAxis *yph = ph->GetYaxis();
979 Int_t nbins = xph->GetNbins();
980 Int_t nybins = yph->GetNbins();
981 if (!InitFit(nbins,1)) {
984 fStatisticMean = 0.0;
996 // Init fCountDet and fCount
997 InitfCountDetAndfCount(1);
999 // Beginning of the loop
1000 for (Int_t idect = fDect1[1]; idect < fDect2[1]; idect++) {
1002 TH1D *projph = (TH1D *) ph->ProjectionY("projph",idect+1,idect+1,(Option_t *) "e");
1003 projph->SetDirectory(0);
1005 // Number of entries for this calibration group
1006 Double_t nentries = 0;
1007 for (Int_t k = 0; k < nybins; k++) {
1008 nentries += ph->GetBinEntries(ph->GetBin(idect+1,k+1));
1014 // Determination of fNnZ, fNnRphi, fNfragZ and fNfragRphi
1015 UpdatefCountDetAndfCount(idect,1);
1017 // Reconstruction of the row and pad group: rowmin, row max ...
1018 ReconstructFitRowMinRowMax(idect,1);
1020 // Rebin and statistic stuff
1021 // This detector has not enough statistics or was off
1022 if (nentries < fMinEntries) {
1024 // Fill with the default values
1025 NotEnoughStatistic(idect,1);
1036 // Statistics of the histos fitted
1037 AliInfo(Form("For the group number %d there are %f stats",idect,nentries));
1039 fStatisticMean += nentries;
1041 // Calcul of "real" coef
1042 CalculVdriftCoefMean(fCountDet[1],(Int_t) (idect - fDect1[1]));
1043 CalculT0CoefMean(fCountDet[1],(Int_t) (idect - fDect1[1]));
1045 // Method Mean and fit
1046 // idect is egal for fDebug = 0 and 2, only to fill the hist
1047 FitPente((TH1 *) projph,(Int_t) (idect - fDect1[1]));
1050 // idect is egal for fDebug = 0 and 2, only to fill the hist
1052 FitPH((TH1 *) projph,(Int_t) (idect - fDect1[1]));
1055 // Visualise the detector for fDebug 3 or 4
1056 // Here is the reconstruction of the pad and row group is used!
1062 // Fill the tree if end of a detector or only the pointer to the branch!!!
1063 FillInfosFit(idect,1);
1074 // 0 no plot, 1 and 4 error plot, 3 and 4 DB plot
1075 if ((fDebug == 1) ||
1080 if ((fDebug == 4) ||
1087 if (fNumberFit > 0) {
1088 AliInfo(Form("There is a mean statistic of: %d",(Int_t) fStatisticMean / fNumberFit));
1089 fStatisticMean = fStatisticMean / fNumberFit;
1092 AliInfo("There is no fit!");
1095 // Write the things!
1104 //____________Functions fit Online PH2d________________________________________
1105 Bool_t AliTRDCalibra::FitPHOnline()
1108 // Reconstruct the average pulse height from the vectorPH for each
1109 // calibration group
1110 // Fit or use the slope of the average pulse height to reconstruct the
1111 // drift velocity write the results in a tree
1112 // A first calibration of T0 is also made using the same method (slope method)
1115 // Number of Xbins (detectors or groups of pads)
1116 if (!InitFit(0,1)) {
1119 fStatisticMean = 0.0;
1123 // Init fCountDet and fCount
1124 InitfCountDetAndfCount(1);
1126 // Beginning of the loop
1127 for (Int_t idect = fDect1[1]; idect < fDect2[1]; idect++) {
1129 // Search if the group is in the VectorCH
1130 Int_t place = SearchInVector(idect,1);
1140 AliTRDPInfo *fPHInfo = new AliTRDPInfo();
1142 fPHInfo = (AliTRDPInfo *) fVectorPH->At(place);
1143 projph = CorrectTheError((TGraphErrors *) ConvertVectorPHisto(fPHInfo,(const char *) name));
1144 projph->SetDirectory(0);
1148 // Determination of fNnZ, fNnRphi, fNfragZ and fNfragRphi
1149 UpdatefCountDetAndfCount(idect,1);
1151 // Reconstruction of the row and pad group: rowmin, row max ...
1152 ReconstructFitRowMinRowMax(idect,1);
1154 // Rebin and statistic stuff
1155 // This detector has not enough statistics or was off
1156 if ((place == -1) ||
1158 (fEntriesCurrent < fMinEntries))) {
1160 // Fill with the default values
1161 NotEnoughStatistic(idect,1);
1172 // Statistic of the histos fitted
1173 AliInfo(Form("For the group number %d there are %d stats",idect,fEntriesCurrent));
1175 fStatisticMean += fEntriesCurrent;
1177 // Calcul of "real" coef
1178 CalculVdriftCoefMean(fCountDet[1],(Int_t) (idect - fDect1[1]));
1179 CalculT0CoefMean(fCountDet[1],(Int_t) (idect - fDect1[1]));
1181 // Method Mean and fit
1182 // idect is egal for fDebug = 0 and 2, only to fill the hist
1183 FitPente((TH1 *) projph,(Int_t) (idect - fDect1[1]));
1186 // idect is egal for fDebug = 0 and 2, only to fill the hist
1188 FitPH((TH1 *) projph,(Int_t) (idect - fDect1[1]));
1191 // Visualise the detector for fDebug 3 or 4
1192 // Here is the reconstruction of the pad and row group is used!
1198 // Fill the tree if end of a detector or only the pointer to the branch!!!
1199 FillInfosFit(idect,1);
1210 // 0 no plot, 1 and 4 error plot, 3 and 4 DB plot
1211 if ((fDebug == 1) ||
1216 if ((fDebug == 4) ||
1223 if (fNumberFit > 0) {
1224 AliInfo(Form("There is a mean statistic of: %d",(Int_t) fStatisticMean / fNumberFit));
1225 fStatisticMean = fStatisticMean / fNumberFit;
1228 AliInfo("There is no fit!");
1231 // Write the things!
1232 if (fWriteCoef[1]) {
1240 //____________Functions fit Online PH2d________________________________________
1241 Bool_t AliTRDCalibra::FitPHOnline(TTree *tree)
1244 // Look if the calibration group can be found in the tree, if yes take the
1245 // histo, fit it, and write the results in a tree
1246 // A first calibration of T0 is also made using the same method (slope method)
1249 // Number of Xbins (detectors or groups of pads)
1250 if (!InitFit(0,1)) {
1253 fStatisticMean = 0.0;
1265 // Init fCountDet and fCount
1266 InitfCountDetAndfCount(1);
1267 TGraphErrors *projphtree = 0x0;
1268 tree->SetBranchAddress("histo",&projphtree);
1269 TObjArray *vectorplace = ConvertTreeVector(tree);
1271 // Beginning of the loop
1272 for (Int_t idect = fDect1[1]; idect < fDect2[1]; idect++) {
1274 // Search if the group is in the VectorCH
1275 Int_t place = SearchInTreeVector(vectorplace,idect);
1283 tree->GetEntry(place);
1284 projph = CorrectTheError(projphtree);
1287 // Determination of fNnZ, fNnRphi, fNfragZ and fNfragRphi
1288 UpdatefCountDetAndfCount(idect,1);
1290 // Reconstruction of the row and pad group: rowmin, row max ...
1291 ReconstructFitRowMinRowMax(idect,1);
1293 // Rebin and statistic stuff
1294 // This detector has not enough statistics or was off
1297 (fEntriesCurrent < fMinEntries))) {
1299 // Fill with the default values
1300 NotEnoughStatistic(idect,1);
1311 // Statistics of the histos fitted
1312 AliInfo(Form("For the group number %d there are %d stats",idect,fEntriesCurrent));
1314 fStatisticMean += fEntriesCurrent;
1316 // Calcul of "real" coef
1317 CalculVdriftCoefMean(fCountDet[1],(Int_t) (idect - fDect1[1]));
1318 CalculT0CoefMean(fCountDet[1],(Int_t) (idect - fDect1[1]));
1320 // Method Mean and fit
1321 // idect is egal for fDebug = 0 and 2, only to fill the hist
1322 FitPente((TH1 *) projph,(Int_t) (idect - fDect1[1]));
1324 // idect is egal for fDebug = 0 and 2, only to fill the hist
1326 FitPH((TH1 *) projph,(Int_t) (idect - fDect1[1]));
1329 // Visualise the detector for fDebug 3 or 4
1330 // Here is the reconstruction of the pad and row group is used!
1336 // Fill the tree if end of a detector or only the pointer to the branch!!!
1337 FillInfosFit(idect,1);
1347 // 0 no plot, 1 and 4 error plot, 3 and 4 DB plot
1348 if ((fDebug == 1) ||
1353 if ((fDebug == 4) ||
1360 if (fNumberFit > 0) {
1361 AliInfo(Form("There is a mean statistic of: %d",(Int_t) fStatisticMean / fNumberFit));
1362 fStatisticMean = fStatisticMean / fNumberFit;
1365 AliInfo("There is no fit!");
1368 // Write the things!
1369 if (fWriteCoef[1]) {
1377 //____________Functions fit Online PRF2d_______________________________________
1378 Bool_t AliTRDCalibra::FitPRFOnline(TProfile2D *prf)
1381 // Take the 1D profiles (pad response function), projections of the 2D PRF
1382 // on the Xaxis, for each calibration group
1383 // Fit with a gaussian to reconstruct the sigma of the pad response function
1384 // write the results in a tree
1387 // Number of Xbins (detectors or groups of pads)
1388 TAxis *xprf = prf->GetXaxis();
1389 TAxis *yprf = prf->GetYaxis();
1390 Int_t nybins = yprf->GetNbins();
1391 Int_t nbins = xprf->GetNbins();
1392 if (!InitFit(nbins,2)) {
1395 fStatisticMean = 0.0;
1401 fVectorPRF->Clear();
1407 // Init fCountDet and fCount
1408 InitfCountDetAndfCount(2);
1410 // Beginning of the loop
1411 for (Int_t idect = fDect1[2]; idect < fDect2[2]; idect++) {
1413 TH1D *projprf = (TH1D *) prf->ProjectionY("projprf",idect+1,idect+1,(Option_t *) "e");
1414 projprf->SetDirectory(0);
1416 // Number of entries for this calibration group
1417 Double_t nentries = 0;
1418 for (Int_t k = 0; k < nybins; k++) {
1419 nentries += prf->GetBinEntries(prf->GetBin(idect+1,k+1));
1421 if(nentries > 0) fNumberEnt++;
1423 // Determination of fNnZ, fNnRphi, fNfragZ and fNfragRphi
1424 UpdatefCountDetAndfCount(idect,2);
1426 // Reconstruction of the row and pad group: rowmin, row max ...
1427 ReconstructFitRowMinRowMax(idect,2);
1429 // Rebin and statistic stuff
1430 // This detector has not enough statistics or was off
1431 if (nentries < fMinEntries) {
1433 // Fill with the default values
1434 NotEnoughStatistic(idect,2);
1445 // Statistics of the histos fitted
1446 AliInfo(Form("For the group number %d there are %f stats",idect,nentries));
1448 fStatisticMean += nentries;
1450 // Calcul of "real" coef
1451 if ((fDebug == 1) ||
1453 CalculPRFCoefMean(fCountDet[2],(Int_t) (idect - fDect1[2]));
1456 // Method Mean and fit
1457 // idect is egal for fDebug = 0 and 2, only to fill the hist
1458 FitPRF((TH1 *) projprf,(Int_t) (idect - fDect1[2]));
1460 // Visualise the detector for fDebug 3 or 4
1461 // Here is the reconstruction of the pad and row group is used!
1466 // Fill the tree if end of a detector or only the pointer to the branch!!!
1467 FillInfosFit(idect,2);
1477 // No plot, 1 and 4 error plot, 3 and 4 DB plot
1478 if ((fDebug == 1) ||
1482 if ((fDebug == 4) ||
1488 if (fNumberFit > 0) {
1489 AliInfo(Form("There is a mean statistic of: %d",(Int_t) fStatisticMean / fNumberFit));
1490 fStatisticMean = fStatisticMean / fNumberFit;
1493 AliInfo("There is no fit!");
1496 // Write the things!
1497 if (fWriteCoef[2]) {
1505 //____________Functions fit Online PRF2d_______________________________________
1506 Bool_t AliTRDCalibra::FitPRFOnline(TTree *tree)
1509 // Look if the calibration group can be found in the tree, if yes take
1510 // the histo, fit it, and write the results in a tree
1513 // Number of Xbins (detectors or groups of pads)
1514 if (!InitFit(0,2)) {
1517 fStatisticMean = 0.0;
1523 fVectorPRF->Clear();
1529 // Init fCountDet and fCount
1530 InitfCountDetAndfCount(2);
1531 TGraphErrors *projprftree = 0x0;
1532 tree->SetBranchAddress("histo",&projprftree);
1533 TObjArray *vectorplace = ConvertTreeVector(tree);
1535 // Beginning of the loop
1536 for (Int_t idect = fDect1[2]; idect < fDect2[2]; idect++) {
1538 // Search if the group is in the VectorCH
1539 Int_t place = SearchInTreeVector(vectorplace,idect);
1542 TH1F *projprf = 0x0;
1547 tree->GetEntry(place);
1548 projprf = CorrectTheError(projprftree);
1551 // Determination of fNnZ, fNnRphi, fNfragZ and fNfragRphi
1552 UpdatefCountDetAndfCount(idect,2);
1554 // Reconstruction of the row and pad group: rowmin, row max ...
1555 ReconstructFitRowMinRowMax(idect,2);
1557 // Rebin and statistic stuff
1558 // This detector has not enough statistics or was off
1559 if ((place == -1) ||
1561 (fEntriesCurrent < fMinEntries))) {
1563 // Fill with the default values
1564 NotEnoughStatistic(idect,2);
1575 // Statistics of the histos fitted
1576 AliInfo(Form("For the group number %d there are %d stats",idect,fEntriesCurrent));
1578 fStatisticMean += fEntriesCurrent;
1580 // Calcul of "real" coef
1581 if ((fDebug == 1) ||
1583 CalculPRFCoefMean(fCountDet[2],(Int_t) (idect - fDect1[2]));
1586 // Method Mean and fit
1587 // idect is egal for fDebug = 0 and 2, only to fill the hist
1588 FitPRF((TH1 *) projprf,(Int_t) (idect - fDect1[2]));
1589 // Visualise the detector for fDebug 3 or 4
1590 // Here is the reconstruction of the pad and row group is used!
1594 // Fill the tree if end of a detector or only the pointer to the branch!!!
1595 FillInfosFit(idect,2);
1605 // No plot, 1 and 4 error plot, 3 and 4 DB plot
1606 if ((fDebug == 1) ||
1610 if ((fDebug == 4) ||
1616 if (fNumberFit > 0) {
1617 AliInfo(Form("There is a mean statistic of: %d",(Int_t) fStatisticMean / fNumberFit));
1618 fStatisticMean = fStatisticMean / fNumberFit;
1621 AliInfo("There is no fit!");
1624 // Write the things!
1625 if (fWriteCoef[2]) {
1633 //____________Functions fit Online PRF2d_______________________________________
1634 Bool_t AliTRDCalibra::FitPRFOnline()
1637 // Reconstruct the 1D histo (pad response function) from the vectorPRD for
1638 // each calibration group
1639 // Fit with a gaussian to reconstruct the sigma of the pad response function
1640 // write the results in a tree
1643 // Number of Xbins (detectors or groups of pads)
1644 if (!InitFit(0,2)) {
1647 fStatisticMean = 0.0;
1651 // Init fCountDet and fCount
1652 InitfCountDetAndfCount(2);
1654 // Beginning of the loop
1655 for (Int_t idect = fDect1[2]; idect < fDect2[2]; idect++) {
1657 // Search if the group is in the VectorCH
1658 Int_t place = SearchInVector(idect,2);
1661 TH1F *projprf = 0x0;
1662 TString name("PRF");
1668 AliTRDPInfo *fPRFInfo = new AliTRDPInfo();
1670 fPRFInfo = (AliTRDPInfo *) fVectorPRF->At(place);
1671 projprf = CorrectTheError((TGraphErrors *) ConvertVectorPHisto(fPRFInfo,(const char *)name));
1672 projprf->SetDirectory(0);
1676 // Determination of fNnZ, fNnRphi, fNfragZ and fNfragRphi
1677 UpdatefCountDetAndfCount(idect,2);
1679 // Reconstruction of the row and pad group: rowmin, row max ...
1680 ReconstructFitRowMinRowMax(idect,2);
1682 // Rebin and statistic stuff
1683 // This detector has not enough statistics or was off
1684 if ((place == -1) ||
1686 (fEntriesCurrent < fMinEntries))) {
1688 // Fill with the default values
1689 NotEnoughStatistic(idect,2);
1700 // Statistic of the histos fitted
1701 AliInfo(Form("For the group number %d there are %d stats", idect,fEntriesCurrent));
1703 fStatisticMean += fEntriesCurrent;
1705 // Calcul of "real" coef
1706 if ((fDebug == 1) ||
1708 CalculPRFCoefMean(fCountDet[2],(Int_t) (idect-fDect1[2]));
1711 // Method Mean and fit
1712 // idect is egal for fDebug = 0 and 2, only to fill the hist
1713 FitPRF((TH1 *) projprf,(Int_t) (idect-fDect1[2]));
1714 // Visualise the detector for fDebug 3 or 4
1715 // Here is the reconstruction of the pad and row group is used!
1719 // Fill the tree if end of a detector or only the pointer to the branch!!!
1720 FillInfosFit(idect,2);
1730 // No plot, 1 and 4 error plot, 3 and 4 DB plot
1731 if ((fDebug == 1) ||
1735 if ((fDebug == 4) ||
1741 if (fNumberFit > 0) {
1742 AliInfo(Form("There is a mean statistic of: %d",(Int_t) fStatisticMean / fNumberFit));
1745 AliInfo("There is no fit!");
1748 // Write the things!
1749 if (fWriteCoef[2]) {
1757 //____________Functions for initialising the AliTRDCalibra in the code_________
1758 Bool_t AliTRDCalibra::Init2Dhistos()
1761 // For the offline tracking
1762 // This function will be called in the function AliReconstruction::Run()
1763 // Init the calibration mode (Nz, Nrphi), the 2D histograms if fHisto2d = kTRUE,
1768 AliTRDcalibDB *cal = AliTRDcalibDB::Instance();
1770 AliInfo("Could not get calibDB");
1773 AliTRDCommonParam *parCom = AliTRDCommonParam::Instance();
1775 AliInfo("Could not get CommonParam");
1780 fTimeMax = cal->GetNumberOfTimeBins();
1781 fSf = parCom->GetSamplingFrequency();
1782 if (fRelativeScaleAuto) {
1786 fRelativeScale = 20;
1789 // Create the 2D histos corresponding to the pad groupCalibration mode
1792 AliInfo(Form("The pad calibration mode for the relative gain calibration: Nz %d, and Nrphi %d"
1796 // Calcul the number of Xbins
1798 ModePadCalibration(2,0);
1799 ModePadFragmentation(0,2,0,0);
1800 fDetChamb2[0] = fNfragZ[0] * fNfragRphi[0];
1802 AliInfo(Form("For the chamber 2: %d",fDetChamb2[0]));
1804 fNtotal[0] += 6 * 18 * fDetChamb2[0];
1805 ModePadCalibration(0,0);
1806 ModePadFragmentation(0,0,0,0);
1807 fDetChamb0[0] = fNfragZ[0] * fNfragRphi[0];
1809 AliInfo(Form("For the other chamber 0: %d",fDetChamb0[0]));
1811 fNtotal[0] += 6 * 4 * 18 * fDetChamb0[0];
1812 AliInfo(Form("Total number of Xbins: %d",fNtotal[0]));
1814 // Create the 2D histo
1816 CreateCH2d(fNtotal[0]);
1819 fVectorCH = new TObjArray();
1820 fPlaCH = new TObjArray();
1824 fAmpTotal = new Float_t[TMath::Max(fDetChamb2[0],fDetChamb0[0])];
1825 for (Int_t k = 0; k < TMath::Max(fDetChamb2[0],fDetChamb0[0]); k++) {
1833 AliInfo(Form("The pad calibration mode for the drift velocity calibration: Nz %d, and Nrphi %d"
1837 // Calcul the number of Xbins
1839 ModePadCalibration(2,1);
1840 ModePadFragmentation(0,2,0,1);
1841 fDetChamb2[1] = fNfragZ[1]*fNfragRphi[1];
1843 AliInfo(Form("For the chamber 2: %d",fDetChamb2[1]));
1845 fNtotal[1] += 6 * 18 * fDetChamb2[1];
1846 ModePadCalibration(0,1);
1847 ModePadFragmentation(0,0,0,1);
1848 fDetChamb0[1] = fNfragZ[1] * fNfragRphi[1];
1850 AliInfo(Form("For the chamber 0: %d",fDetChamb0[1]));
1852 fNtotal[1] += 6 * 4 * 18 * fDetChamb0[1];
1853 AliInfo(Form("Total number of Xbins: %d",fNtotal[1]));
1855 // Create the 2D histo
1857 CreatePH2d(fNtotal[1]);
1860 fVectorPH = new TObjArray();
1861 fPlaPH = new TObjArray();
1865 fPHPlace = new Short_t[fTimeMax];
1866 for (Int_t k = 0; k < fTimeMax; k++) {
1869 fPHValue = new Float_t[fTimeMax];
1870 for (Int_t k = 0; k < fTimeMax; k++) {
1878 AliInfo(Form("The pad calibration mode for the PRF calibration: Nz %d, and Nrphi %d"
1882 // Calcul the number of Xbins
1884 ModePadCalibration(2,2);
1885 ModePadFragmentation(0,2,0,2);
1886 fDetChamb2[2] = fNfragZ[2] * fNfragRphi[2];
1888 AliInfo(Form("For the chamber 2: %d",fDetChamb2[2]));
1890 fNtotal[2] += 6 * 18 * fDetChamb2[2];
1891 ModePadCalibration(0,2);
1892 ModePadFragmentation(0,0,0,2);
1893 fDetChamb0[2] = fNfragZ[2] * fNfragRphi[2];
1895 AliInfo(Form("For the chamber 0: %d",fDetChamb0[2]));
1897 fNtotal[2] += 6 * 4 * 18 * fDetChamb0[2];
1898 AliInfo(Form("Total number of Xbins: %d",fNtotal[2]));
1900 // Create the 2D histo
1902 CreatePRF2d(fNtotal[2]);
1905 fVectorPRF = new TObjArray();
1906 fPlaPRF = new TObjArray();
1915 //____________Functions for filling the histos in the code_____________________
1917 //____________Offine tracking in the AliTRDtracker_____________________________
1918 Bool_t AliTRDCalibra::ResetTrack()
1921 // For the offline tracking
1922 // This function will be called in the function
1923 // AliTRDtracker::FollowBackPropagation() at the beginning
1924 // Reset the parameter to know we have a new TRD track
1927 fDetectorAliTRDtrack = kFALSE;
1932 //____________Offline tracking in the AliTRDtracker____________________________
1933 Bool_t AliTRDCalibra::UpdateHistograms(AliTRDcluster *cl, AliTRDtrack *t)
1936 // For the offline tracking
1937 // This function will be called in the function
1938 // AliTRDtracker::FollowBackPropagation() in the loop over the clusters
1940 // Fill the 2D histos or the vectors with the info of the clusters at
1941 // the end of a detectors if the track is "good"
1944 // Get the parameter object
1945 AliTRDCommonParam *parCom = AliTRDCommonParam::Instance();
1947 AliInfo("Could not get CommonParam");
1951 // Get the parameter object
1952 AliTRDcalibDB *cal = AliTRDcalibDB::Instance();
1954 AliInfo("Could not get calibDB");
1958 // Localisation of the detector
1959 Int_t detector = cl->GetDetector();
1960 Int_t chamber = GetChamber(detector);
1961 Int_t plane = GetPlane(detector);
1963 // Fill the infos for the previous clusters if not the same
1964 // detector anymore or if not the same track
1965 if (((detector != fDetectorPreviousTrack) || (!fDetectorAliTRDtrack)) &&
1966 (fDetectorPreviousTrack != -1)) {
1970 // If the same track, then look if the previous detector is in
1971 // the same plane, if yes: not a good track
1972 if (fDetectorAliTRDtrack &&
1973 (GetPlane(detector) <= GetPlane(fDetectorPreviousTrack))) {
1974 fGoodTrack = kFALSE;
1977 // Fill only if the track doesn't touch a masked pad or doesn't
1978 // appear in the middle (fGoodTrack)
1983 FillTheInfoOfTheTrackCH();
1988 FillTheInfoOfTheTrackPH();
1991 } // if a good track
1995 } // Fill at the end the charge
1997 // Calcul the position of the detector
1998 if (detector != fDetectorPreviousTrack) {
1999 LocalisationDetectorXbins(detector);
2002 // Reset the good track for the PRF
2003 Bool_t good = kTRUE;
2005 // Localisation of the cluster
2006 Double_t pos[3] = { 0.0, 0.0, 0.0 };
2007 pos[0] = cl->GetX();
2008 pos[1] = cl->GetY();
2009 pos[2] = cl->GetZ();
2010 Int_t time = cl->GetLocalTimeBin();
2012 // Reset the detector
2013 fDetectorPreviousTrack = detector;
2014 fDetectorAliTRDtrack = kTRUE;
2016 // Position of the cluster
2017 AliTRDpadPlane *padplane = parCom->GetPadPlane(plane,chamber);
2018 Int_t row = padplane->GetPadRowNumber(pos[2]);
2019 Double_t offsetz = padplane->GetPadRowOffset(row,pos[2]);
2020 Double_t offsettilt = padplane->GetTiltOffset(offsetz);
2021 Int_t col = padplane->GetPadColNumber(pos[1] + offsettilt,offsetz);
2023 // See if we are not near a masked pad
2024 if (!IsPadOn(detector,col,row)) {
2026 fGoodTrack = kFALSE;
2030 if (!IsPadOn(detector,col-1,row)) {
2031 fGoodTrack = kFALSE;
2037 if (!IsPadOn(detector,col+1,row)) {
2038 fGoodTrack = kFALSE;
2043 // Row of the cluster and position in the pad groups
2044 Int_t posr[3] = { 0, 0, 0 };
2045 if ((fCH2dOn) && (fNnZ[0] != 0)) {
2046 posr[0] = (Int_t) row / fNnZ[0];
2048 if ((fPH2dOn) && (fNnZ[1] != 0)) {
2049 posr[1] = (Int_t) row / fNnZ[1];
2051 if ((fPRF2dOn) && (fNnZ[2] != 0)) {
2052 posr[2] = (Int_t) row / fNnZ[2];
2055 // Col of the cluster and position in the pad groups
2056 Int_t posc[3] = { 0, 0, 0 };
2057 if ((fCH2dOn) && (fNnRphi[0] != 0)) {
2058 posc[0] = (Int_t) col / fNnRphi[0];
2060 if ((fPH2dOn) && (fNnRphi[1] != 0)) {
2061 posc[1] = (Int_t) col / fNnRphi[1];
2063 if ((fPRF2dOn) && (fNnRphi[2] != 0)) {
2064 posc[2] = (Int_t) col / fNnRphi[2];
2067 // Charge in the cluster
2068 // For the moment take the abs
2069 Float_t q = TMath::Abs(cl->GetQ());
2070 Short_t *signals = cl->GetSignals();
2072 // Correction due to the track angle
2073 Float_t correction = 1.0;
2074 Float_t normalisation = 6.67;
2075 if ((q >0) && (t->GetNdedx() > 0)) {
2076 correction = t->GetClusterdQdl((t->GetNdedx() - 1)) / (q * normalisation);
2079 // Fill the fAmpTotal with the charge
2082 fAmpTotal[(Int_t) (posc[0]*fNfragZ[0]+posr[0])] += q * correction;
2085 fAmpTotal[(Int_t) (posc[0]*fNfragZ[0]+posr[0])] += ((Float_t) signals[3]) * correction;
2089 // Fill the fPHPlace and value
2091 fPHPlace[time] = posc[1]*fNfragZ[1]+posr[1];
2093 fPHValue[time] = q * correction;
2096 fPHValue[time] = ((Float_t) signals[3]) * correction;
2100 // Fill direct the PRF
2101 if ((fPRF2dOn) && (good)) {
2103 Float_t yminus = 0.0;
2104 Float_t xcenter = 0.0;
2105 Float_t ycenter = 0.0;
2107 Bool_t echec = kFALSE;
2109 if ((cl->From3pad()) && (!cl->IsUsed())) {
2111 // Center 3 balanced
2112 if ((((Float_t) signals[3]) > fThresholdClusterPRF2) &&
2113 (((Float_t) signals[2]) > fThresholdClusterPRF2) &&
2114 (((Float_t) signals[4]) > fThresholdClusterPRF2) &&
2115 (((Float_t) signals[1]) < fThresholdClusterPRF1) &&
2116 (((Float_t) signals[5]) < fThresholdClusterPRF1) &&
2117 ((((Float_t) signals[2])*((Float_t) signals[4])/(((Float_t) signals[3])*((Float_t) signals[3]))) < 0.06)) {
2118 // Col correspond to signals[3]
2119 if (fCenterOfflineCluster) {
2120 xcenter = cl->GetCenter();
2123 // Security of the denomiateur is 0
2124 if ((((Float_t) (((Float_t) signals[3]) * ((Float_t) signals[3]))) /
2125 ((Float_t) (((Float_t) signals[2]) * ((Float_t) signals[4])))) != 1.0) {
2126 xcenter = 0.5 * (TMath::Log((Float_t) (((Float_t) signals[4]) / ((Float_t) signals[2]))))
2127 / (TMath::Log(((Float_t) (((Float_t) signals[3]) * ((Float_t) signals[3])))
2128 / ((Float_t) (((Float_t) signals[2]) * ((Float_t) signals[4])))));
2134 if ((xcenter > -0.5) && (xcenter < 0.5)) {
2135 ycenter = (Float_t) (((Float_t) signals[3])
2136 / (((Float_t) signals[2]) + ((Float_t) signals[3]) + (((Float_t) signals[4]))));
2137 yminus = (Float_t) (((Float_t) signals[2])
2138 / (((Float_t) signals[2]) + ((Float_t) signals[3]) + (((Float_t) signals[4]))));
2139 ymax = (Float_t) (((Float_t) signals[4])
2140 / (((Float_t) signals[2]) + ((Float_t) signals[3]) + (((Float_t) signals[4]))));
2141 if ((TMath::Abs(((Float_t) signals[2]) + ((Float_t) signals[3]) + (((Float_t) signals[4])) - q) < 10.0)) {
2147 // Fill only if it is in the drift region!
2148 if ((((Float_t) (((Float_t) time) / fSf)) > 0.3) && (echec)) {
2150 fPRF2d->Fill((fXbins[2]+posc[2]*fNfragZ[2]+posr[2]+0.5),xcenter,ycenter);
2151 if (xcenter < 0.0) {
2152 fPRF2d->Fill((fXbins[2]+posc[2]*fNfragZ[2]+posr[2]+0.5),-(xcenter+1.0),yminus);
2154 if (xcenter > 0.0) {
2155 fPRF2d->Fill((fXbins[2]+posc[2]*fNfragZ[2]+posr[2]+0.5),1.0-xcenter,ymax);
2159 UpdateVectorPRF(fXbins[2]+posc[2]*fNfragZ[2]+posr[2],xcenter,ycenter);
2160 if (xcenter < 0.0) {
2161 UpdateVectorPRF(fXbins[2]+posc[2]*fNfragZ[2]+posr[2],-(xcenter+1.0),yminus);
2163 if (xcenter > 0.0) {
2164 UpdateVectorPRF(fXbins[2]+posc[2]*fNfragZ[2]+posr[2],1.0-xcenter,ymax);
2167 } // If in the drift region
2177 //____________Online trackling in AliTRDtrigger________________________________
2178 Bool_t AliTRDCalibra::UpdateHistogramcm(AliTRDmcmTracklet *trk)
2182 // This function will be called in the function AliTRDtrigger::TestTracklet
2183 // before applying the pt cut on the tracklets
2184 // Fill the infos for the tracklets fTrkTest if the tracklets is "good"
2187 // Localisation of the Xbins involved
2188 Int_t idect = trk->GetDetector();
2189 LocalisationDetectorXbins(idect);
2191 // Get the parameter object
2192 AliTRDcalibDB *cal = AliTRDcalibDB::Instance();
2194 AliInfo("Could not get calibDB");
2201 // Row of the tracklet and position in the pad groups
2202 Int_t row = trk->GetRow();
2203 Int_t posr[3] = { 0, 0, 0 };
2204 if ((fCH2dOn) && (fNnZ[0] != 0)) {
2205 posr[0] = (Int_t) row / fNnZ[0];
2207 if ((fPH2dOn) && (fNnZ[1] != 0)) {
2208 posr[1] = (Int_t) row / fNnZ[1];
2210 if ((fPRF2dOn) && (fNnZ[2] != 0)) {
2211 posr[2] = (Int_t) row / fNnZ[2];
2214 // Eventuelle correction due to track angle in z direction
2215 Float_t correction = 1.0;
2216 if (fMcmCorrectAngle) {
2217 Float_t z = trk->GetRowz();
2218 Float_t r = trk->GetTime0();
2219 correction = r / TMath::Sqrt((r*r+z*z));
2222 // Boucle sur les clusters
2223 // Condition on number of cluster: don't come from the middle of the detector
2224 if (trk->GetNclusters() >= fNumberClusters) {
2226 for (Int_t icl = 0; icl < trk->GetNclusters(); icl++) {
2228 Float_t amp[3] = { 0.0, 0.0, 0.0 };
2229 Int_t time = trk->GetClusterTime(icl);
2230 Int_t col = trk->GetClusterCol(icl);
2232 amp[0] = trk->GetClusterADC(icl)[0] * correction;
2233 amp[1] = trk->GetClusterADC(icl)[1] * correction;
2234 amp[2] = trk->GetClusterADC(icl)[2] * correction;
2236 if ((amp[0] < 0.0) ||
2242 // Col of cluster and position in the pad groups
2243 Int_t posc[3] = { 0, 0, 0 };
2244 if ((fCH2dOn) && (fNnRphi[0] != 0)) {
2245 posc[0] = (Int_t) col / fNnRphi[0];
2247 if ((fPH2dOn) && (fNnRphi[1] != 0)) {
2248 posc[1] = (Int_t) col / fNnRphi[1];
2250 if ((fPRF2dOn) && (fNnRphi[2] != 0)) {
2251 posc[2] = (Int_t) col / fNnRphi[2];
2254 // See if we are not near a masked pad
2255 Bool_t good = kTRUE;
2256 if (!IsPadOn(idect,col,row)) {
2257 fGoodTrack = kFALSE;
2262 if (!IsPadOn(idect,col-1,row)) {
2263 fGoodTrack = kFALSE;
2269 if (!IsPadOn(idect,col+1,row)) {
2270 fGoodTrack = kFALSE;
2277 fPHPlace[time] = posc[1] * fNfragZ[1] + posr[1];
2282 fAmpTotal[(Int_t) (posc[0]*fNfragZ[0]+posr[0])] += (Float_t) (amp[0]+amp[1]+amp[2]);
2285 fPHValue[time] = (Float_t) (amp[0]+amp[1]+amp[2]);
2290 fAmpTotal[(Int_t) (posc[0]*fNfragZ[0]+posr[0])] += (Float_t) amp[1];
2293 fPHValue[time] = amp[1];
2298 if (fPRF2dOn && good) {
2299 if ((amp[0] > fThresholdClusterPRF2) &&
2300 (amp[1] > fThresholdClusterPRF2) &&
2301 (amp[2] > fThresholdClusterPRF2) &&
2302 ((amp[0]*amp[2]/(amp[1]*amp[1])) < 0.06)) {
2303 // Security of the denomiateur is 0
2304 if ((((Float_t) (((Float_t) amp[1]) * ((Float_t) amp[1])))
2305 / ((Float_t) (((Float_t) amp[0]) * ((Float_t) amp[2])))) != 1.0) {
2306 Float_t xcenter = 0.5 * (TMath::Log(amp[2] / amp[0]))
2307 / (TMath::Log((amp[1]*amp[1]) / (amp[0]*amp[2])));
2308 Float_t ycenter = amp[1] / (amp[0] + amp[1] + amp[2]);
2309 if ((xcenter > -0.5) &&
2311 Float_t yminus = amp[0] / (amp[0]+amp[1]+amp[2]);
2312 Float_t ymax = amp[2] / (amp[0]+amp[1]+amp[2]);
2313 // Fill only if it is in the drift region!
2314 if (((Float_t) time / fSf) > 0.3) {
2316 fPRF2d->Fill((fXbins[2]+posc[2]*fNfragZ[2]+posr[2]+0.5),xcenter,ycenter);
2317 if (xcenter < 0.0) {
2318 fPRF2d->Fill((fXbins[2]+posc[2]*fNfragZ[2]+posr[2]+0.5),-(xcenter+1.0),yminus);
2320 if (xcenter > 0.0) {
2321 fPRF2d->Fill((fXbins[2]+posc[2]*fNfragZ[2]+posr[2]+0.5),(1.0-xcenter),ymax);
2325 UpdateVectorPRF((fXbins[2]+posc[2]*fNfragZ[2]+posr[2]),xcenter,ycenter);
2326 if (xcenter < 0.0) {
2327 UpdateVectorPRF(fXbins[2]+posc[2]*fNfragZ[2]+posr[2],-(xcenter+1.0),yminus);
2329 if (xcenter > 0.0) {
2330 UpdateVectorPRF(fXbins[2]+posc[2]*fNfragZ[2]+posr[2],(1.0-xcenter),ymax);
2339 } // Boucle clusters
2342 if (fCH2dOn && fGoodTrack) {
2343 FillTheInfoOfTheTrackCH();
2347 if (fPH2dOn && fGoodTrack) {
2348 FillTheInfoOfTheTrackPH();
2351 } // Condition on number of clusters
2357 //____________Functions for seeing if the pad is really okey___________________
2359 //_____________________________________________________________________________
2360 Bool_t AliTRDCalibra::SetModeCalibrationFromTObject(TObject *object, Int_t i)
2363 // Set fNz[i] and fNrphi[i] of the AliTRDCalibra::Instance()
2364 // corresponding to the given TObject
2367 const char *nametitle = object->GetTitle();
2370 const Char_t *patternz0 = "Nz0";
2371 const Char_t *patternz1 = "Nz1";
2372 const Char_t *patternz2 = "Nz2";
2373 const Char_t *patternz3 = "Nz3";
2374 const Char_t *patternz4 = "Nz4";
2375 const Char_t *patternrphi0 = "Nrphi0";
2376 const Char_t *patternrphi1 = "Nrphi1";
2377 const Char_t *patternrphi2 = "Nrphi2";
2378 const Char_t *patternrphi3 = "Nrphi3";
2379 const Char_t *patternrphi4 = "Nrphi4";
2380 const Char_t *patternrphi5 = "Nrphi5";
2381 const Char_t *patternrphi6 = "Nrphi6";
2384 UShort_t testrphi = 0;
2387 if (strstr(nametitle,patternz0)) {
2391 if (strstr(nametitle,patternz1)) {
2395 if (strstr(nametitle,patternz2)) {
2399 if (strstr(nametitle,patternz3)) {
2403 if (strstr(nametitle,patternz4)) {
2409 if (strstr(nametitle,patternrphi0)) {
2413 if (strstr(nametitle,patternrphi1)) {
2417 if (strstr(nametitle,patternrphi2)) {
2421 if (strstr(nametitle,patternrphi3)) {
2425 if (strstr(nametitle,patternrphi4)) {
2429 if (strstr(nametitle,patternrphi5)) {
2433 if (strstr(nametitle,patternrphi6)) {
2438 // Look if all is okey
2451 //_____________________________________________________________________________
2452 Bool_t AliTRDCalibra::IsPadOn(Int_t detector, Int_t col, Int_t row) const
2455 // Look in the choosen database if the pad is On.
2456 // If no the track will be "not good"
2459 // Get the parameter object
2460 AliTRDcalibDB *cal = AliTRDcalibDB::Instance();
2462 AliInfo("Could not get calibDB");
2466 if (!cal->IsChamberInstalled(detector) ||
2467 cal->IsChamberMasked(detector) ||
2468 cal->IsPadMasked(detector,col,row)) {
2477 //____________Functions for plotting the 2D____________________________________
2479 //_____________________________________________________________________________
2480 void AliTRDCalibra::Plot2d()
2483 // Plot the 2D histos
2498 //____________Writing the 2D___________________________________________________
2500 //_____________________________________________________________________________
2501 Bool_t AliTRDCalibra::Write2d()
2504 // Write the 2D histograms or the vectors converted in trees in the file
2505 // "TRD.calibration.root"
2508 TFile *fout = TFile::Open(fWriteName,"RECREATE");
2509 // Check if the file could be opened
2510 if (!fout || !fout->IsOpen()) {
2511 AliInfo("No File found!");
2514 AliInfo(Form("Numbertrack: %d Numberusedch[0]: %d, Numberusedch[1]: %d Numberusedph[0]: %d, Numberusedph[1]: %d"
2519 ,fNumberUsedPh[1]));
2521 TStopwatch stopwatch;
2525 if ((fCH2dOn ) && (fWrite[0])) {
2527 fout->WriteTObject(fCH2d);
2534 TTree *treeCH2d = ConvertVectorCTTreeHisto(fVectorCH,fPlaCH,"treeCH2d",(const char *) name);
2535 fout->WriteTObject(treeCH2d);
2538 if ((fPH2dOn ) && (fWrite[1])) {
2540 fout->WriteTObject(fPH2d);
2547 TTree *treePH2d = ConvertVectorPTreeHisto(fVectorPH,fPlaPH,"treePH2d",(const char *) name);
2548 fout->WriteTObject(treePH2d);
2551 if ((fPRF2dOn ) && (fWrite[2])) {
2553 fout->WriteTObject(fPRF2d);
2560 TTree *treePRF2d = ConvertVectorPTreeHisto(fVectorPRF,fPlaPRF,"treePRF2d",(const char *) name);
2561 fout->WriteTObject(treePRF2d);
2567 AliInfo(Form("Execution time Write2d: R:%.2fs C:%.2fs"
2568 ,stopwatch.RealTime(),stopwatch.CpuTime()));
2574 //_____________________________________________________________________________
2575 AliTRDCalDet *AliTRDCalibra::CreateDetObjectTree(TTree *tree, Int_t i)
2578 // It creates the AliTRDCalDet object from the tree of the coefficient
2579 // for the calibration i (i != 2)
2580 // It takes the mean value of the coefficients per detector
2581 // This object has to be written in the database
2584 // Create the DetObject
2585 AliTRDCalDet *object = 0x0;
2587 object = new AliTRDCalDet("ChamberGainFactor","GainFactor (detector value)");
2590 object = new AliTRDCalDet("ChamberVdrift","TRD drift velocities (detector value)");
2593 object = new AliTRDCalDet("ChamberT0","T0 (detector value)");
2597 Int_t detector = -1;
2598 Float_t values[2304];
2599 tree->SetBranchAddress("detector",&detector);
2601 tree->SetBranchAddress("gainPad",values);
2604 tree->SetBranchAddress("vdrift" ,values);
2607 tree->SetBranchAddress("t0" ,values);
2610 // For calculating the mean
2613 Int_t numberofentries = tree->GetEntries();
2615 if (numberofentries != 540) {
2616 AliInfo("The tree is not complete");
2619 for (Int_t det = 0; det < numberofentries; ++det) {
2620 tree->GetEntry(det);
2621 if (GetChamber(detector) == 2) {
2629 for (Int_t k = 0; k < nto; k++) {
2630 mean += TMath::Abs(values[k]) / nto;
2634 for (Int_t k = 0; k < nto; k++) {
2635 if(k == 0) mean = values[k];
2636 if(mean > values[k]) mean = values[k];
2639 object->SetValue(detector,mean);
2646 //_____________________________________________________________________________
2647 TObject *AliTRDCalibra::CreatePadObjectTree(TTree *tree, Int_t i
2648 , AliTRDCalDet *detobject)
2651 // It Creates the AliTRDCalPad object from the tree of the
2652 // coefficient for the calibration i (i != 2)
2653 // You need first to create the object for the detectors,
2654 // where the mean value is put.
2655 // This object has to be written in the database
2658 // Create the DetObject
2659 AliTRDCalPad *object = 0x0;
2661 object = new AliTRDCalPad("GainFactor","GainFactor (local variations)");
2664 object = new AliTRDCalPad("LocalVdrift","TRD drift velocities (local variations)");
2667 object = new AliTRDCalPad("LocalT0","T0 (local variations)");
2671 Int_t detector = -1;
2672 Float_t values[2304];
2673 tree->SetBranchAddress("detector",&detector);
2675 tree->SetBranchAddress("gainPad",values);
2678 tree->SetBranchAddress("vdrift" ,values);
2681 tree->SetBranchAddress("t0" ,values);
2686 Int_t numberofentries = tree->GetEntries();
2688 if (numberofentries != 540) {
2689 AliInfo("The tree is not complete");
2692 for (Int_t det = 0; det < numberofentries; ++det) {
2693 tree->GetEntry(det);
2694 AliTRDCalROC *calROC = object->GetCalROC(detector);
2695 mean = detobject->GetValue(detector);
2696 if ((mean == 0) && (i != 3)) {
2699 Int_t rowMax = calROC->GetNrows();
2700 Int_t colMax = calROC->GetNcols();
2701 for (Int_t row = 0; row < rowMax; ++row) {
2702 for (Int_t col = 0; col < colMax; ++col) {
2703 if(i != 3) calROC->SetValue(col,row,TMath::Abs(values[(Int_t) (col*rowMax+row)])/mean);
2704 else calROC->SetValue(col,row,values[(Int_t) (col*rowMax+row)]-mean);
2714 //_____________________________________________________________________________
2715 TObject *AliTRDCalibra::CreatePadObjectTree(TTree *tree)
2718 // It Creates the AliTRDCalPad object from the tree of the
2719 // coefficient for the calibration PRF (i = 2)
2720 // This object has to be written in the database
2723 // Create the DetObject
2724 AliTRDCalPad *object = new AliTRDCalPad("PRFWidth","PRFWidth");
2727 Int_t detector = -1;
2728 Float_t values[2304];
2729 tree->SetBranchAddress("detector",&detector);
2730 tree->SetBranchAddress("width" ,values);
2733 Int_t numberofentries = tree->GetEntries();
2735 if (numberofentries != 540) {
2736 AliInfo("The tree is not complete");
2739 for (Int_t det = 0; det < numberofentries; ++det) {
2740 tree->GetEntry(det);
2741 AliTRDCalROC *calROC = object->GetCalROC(detector);
2742 Int_t rowMax = calROC->GetNrows();
2743 Int_t colMax = calROC->GetNcols();
2744 for (Int_t row = 0; row < rowMax; ++row) {
2745 for (Int_t col = 0; col < colMax; ++col) {
2746 calROC->SetValue(col,row,TMath::Abs(values[(Int_t) (col*rowMax+row)]));
2755 //_____________________________________________________________________________
2756 void AliTRDCalibra::SetRelativeScale(Float_t RelativeScale)
2759 // Set the factor that will divide the deposited charge
2760 // to fit in the histo range [0,300]
2763 if (RelativeScale > 0.0) {
2764 fRelativeScale = RelativeScale;
2767 AliInfo("RelativeScale must be strict positif!");
2772 //_____________________________________________________________________________
2773 void AliTRDCalibra::SetNz(Int_t i, Short_t Nz)
2776 // Set the mode of calibration group in the z direction for the parameter i
2784 AliInfo("You have to choose between 0 and 4");
2789 //_____________________________________________________________________________
2790 void AliTRDCalibra::SetNrphi(Int_t i, Short_t Nrphi)
2793 // Set the mode of calibration group in the rphi direction for the parameter i
2801 AliInfo("You have to choose between 0 and 6");
2806 //_____________________________________________________________________________
2807 void AliTRDCalibra::SetPeriodeFitPH(Int_t periodeFitPH)
2810 // Set FitPH if 1 then each detector will be fitted
2813 if (periodeFitPH > 0) {
2814 fFitPHPeriode = periodeFitPH;
2817 AliInfo("periodeFitPH must be higher than 0!");
2822 //_____________________________________________________________________________
2823 void AliTRDCalibra::SetBeginFitCharge(Float_t beginFitCharge)
2826 // The fit of the deposited charge distribution begins at
2827 // histo->Mean()/beginFitCharge
2828 // You can here set beginFitCharge
2831 if (beginFitCharge > 0) {
2832 fBeginFitCharge = beginFitCharge;
2835 AliInfo("beginFitCharge must be strict positif!");
2840 //_____________________________________________________________________________
2841 void AliTRDCalibra::SetT0Shift(Float_t t0Shift)
2844 // The t0 calculated with the maximum positif slope is shift from t0Shift
2845 // You can here set t0Shift
2852 AliInfo("t0Shift must be strict positif!");
2857 //_____________________________________________________________________________
2858 void AliTRDCalibra::SetRangeFitPRF(Float_t rangeFitPRF)
2861 // The fit of the PRF is from -rangeFitPRF to rangeFitPRF
2862 // You can here set rangeFitPRF
2865 if ((rangeFitPRF > 0) &&
2866 (rangeFitPRF <= 1.0)) {
2867 fRangeFitPRF = rangeFitPRF;
2870 AliInfo("rangeFitPRF must be between 0 and 1.0");
2875 //_____________________________________________________________________________
2876 void AliTRDCalibra::SetRebin(Short_t rebin)
2879 // Rebin with rebin time less bins the Ch histo
2880 // You can set here rebin that should divide the number of bins of CH histo
2885 AliInfo("You have to be sure that fRebin divides fNumberBinCharge used!");
2888 AliInfo("You have to choose a positiv value!");
2893 //_____________________________________________________________________________
2894 TTree *AliTRDCalibra::Sum2Trees(const Char_t *filename1
2895 , const Char_t *filename2
2896 , const Char_t *variablecali)
2899 // It returns the sum of two trees with the name variablecali
2900 // in the files filenam1 and filename2 equivalent of merging two 2D histos
2901 // The name of the resulting tree is the same as the two input trees
2902 // variablecali can be treeCH2d, treePH2d or treePRF2d
2906 TChain *treeChain = new TChain(variablecali);
2907 TObjArray *vectorplace = new TObjArray();
2908 TObjArray *where = new TObjArray();
2912 TFile *file1 = new TFile(filename1,"READ");
2913 TTree *tree1 = (TTree *) file1->Get(variablecali);
2918 vectorplace = ConvertTreeVector(tree1);
2920 // Say where it is in tree 1
2921 for (Int_t jui = 0; jui < (Int_t) vectorplace->GetEntriesFast(); jui++) {
2922 AliTRDPlace *placejui = new AliTRDPlace();
2923 placejui->SetPlace(jui);
2924 TObjArray *chainplace = new TObjArray();
2925 chainplace->Add((TObject *) placejui);
2926 where->Add((TObject *) chainplace);
2930 treeChain->Add(filename1);
2935 TFile *file2 = new TFile(filename2,"READ");
2936 TTree *tree2 = (TTree *) file2->Get(variablecali);
2941 TObjArray *vector2 = ConvertTreeVector(tree2);
2942 Int_t j = treeChain->GetEntries();
2944 for (Int_t jui = 0; jui < (Int_t) vector2->GetEntriesFast(); jui++) {
2945 // Search if already found
2946 Int_t place = SearchInTreeVector(vectorplace,((AliTRDPlace *) vector2->At(jui))->GetPlace());
2947 // Create a new element in the two std vectors
2949 AliTRDPlace *placejjui = new AliTRDPlace();
2950 placejjui->SetPlace((j+jui));
2951 TObjArray *chainplace = new TObjArray();
2952 chainplace->Add((TObject *) placejjui);
2953 vectorplace->Add((TObject *) (vector2->At(jui)));
2954 where->Add((TObject *) chainplace);
2956 // Update the element at the place "place" in the std vector whereinthechain
2958 AliTRDPlace *placejjui = new AliTRDPlace();
2959 placejjui->SetPlace((j+jui));
2960 TObjArray *chainplace = ((TObjArray *) where->At(place));
2961 chainplace->Add((TObject *) placejjui);
2962 where->AddAt((TObject *) chainplace,place);
2967 treeChain->Add(filename2);
2970 // Take care of the profile
2971 const Char_t *pattern = "P";
2974 if (!strstr(variablecali,pattern)) {
2976 // Ready to read the chain
2978 treeChain->SetBranchAddress("histo",&his);
2980 // Initialise the final tree
2982 TH1F *histsum = 0x0;
2984 tree = new TTree(variablecali,variablecali);
2985 tree->Branch("groupnumber",&group,"groupnumber/I");
2986 tree->Branch("histo","TH1F",&histsum,32000,0);
2989 if (treeChain->GetEntries() < 1) {
2993 for (Int_t h = 0; h < (Int_t) vectorplace->GetEntriesFast(); h++) {
2994 group = ((AliTRDPlace *) vectorplace->At(h))->GetPlace();
2995 TObjArray *chainplace = ((TObjArray *) where->At(h));
2996 treeChain->GetEntry(((AliTRDPlace *) chainplace->At(0))->GetPlace());
2997 //Init for the first time
2999 histsum = new TH1F("","",his->GetXaxis()->GetNbins()
3000 ,his->GetXaxis()->GetBinLowEdge(1)
3001 ,his->GetXaxis()->GetBinUpEdge(his->GetXaxis()->GetNbins()));
3004 // Reset for each new group
3005 histsum->SetEntries(0.0);
3006 for (Int_t l = 0; l <= histsum->GetXaxis()->GetNbins(); l++) {
3007 histsum->SetBinContent(l,0.0);
3008 histsum->SetBinError(l,0.0);
3010 histsum->Add(his,1);
3011 if ((Int_t) chainplace->GetEntriesFast() > 1) {
3012 for (Int_t s = 1; s < (Int_t) chainplace->GetEntriesFast(); s++) {
3013 treeChain->GetEntry(((AliTRDPlace *) chainplace->At(s))->GetPlace());
3014 histsum->Add(his,1);
3023 // Ready to read the chain
3024 TGraphErrors *his = 0x0;
3025 treeChain->SetBranchAddress("histo",&his);
3027 // Initialise the final tree
3029 TGraphErrors *histsum = 0x0;
3030 Double_t *xref = 0x0;
3032 tree = new TTree(variablecali,variablecali);
3033 tree->Branch("groupnumber",&group,"groupnumber/I");
3034 tree->Branch("histo","TGraphErrors",&histsum,32000,0);
3037 if (treeChain->GetEntries() < 1) {
3041 for (Int_t h = 0; h < (Int_t) vectorplace->GetEntriesFast(); h++) {
3043 group = ((AliTRDPlace *) vectorplace->At(h))->GetPlace();
3044 TObjArray *chainplace = ((TObjArray *) where->At(h));
3045 treeChain->GetEntry(((AliTRDPlace *) chainplace->At(0))->GetPlace());
3046 //Init or reset for a new group
3047 Int_t nbins = his->GetN();
3049 x = new Double_t[nbins];
3052 ex = new Double_t[nbins];
3054 y = new Double_t[nbins];
3056 ey = new Double_t[nbins];
3058 for (Int_t lo = 0; lo < nbins; lo++) {
3065 histsum = new TGraphErrors(nbins,x,y,ex,ey);
3068 histsum = AddProfiles(his,histsum);
3069 if ((Int_t) chainplace->GetEntriesFast() > 1) {
3070 for (Int_t s = 1; s < (Int_t) chainplace->GetEntriesFast(); s++) {
3071 treeChain->GetEntry(((AliTRDPlace *) chainplace->At(s))->GetPlace());
3072 histsum = AddProfiles(his,histsum);
3086 //____________Function fill 2D for the moment out of the code__________________
3088 //____________Function fill 2D all objects from digits_________________________
3089 Bool_t AliTRDCalibra::Create2DDiSimOnline(Int_t iev1, Int_t iev2)
3092 // Only for simulations, after the simulation, create the 2D histos
3093 // from the digits stored in the file "TRD.Digits.root"
3094 // Only for CH and PH
3097 const Int_t kNplan = 6;
3098 const Int_t kNcham = 5;
3100 // RunLoader and so on
3102 delete gAlice->GetRunLoader();
3107 AliRunLoader *rl = AliRunLoader::Open("galice.root");
3113 gAlice = rl->GetAliRun();
3118 // Import the Trees for the event nEvent in the file
3119 rl->LoadKinematics();
3123 AliLoader *loader = rl->GetLoader("TRDLoader");
3125 AliInfo("No TRDLLoader found!");
3129 // Get the pointer to the TRD detector
3130 AliTRD *trd = (AliTRD *) gAlice->GetDetector("TRD");
3132 AliInfo("No TRD detector found");
3136 // Get the pointer to the geometry object
3137 AliTRDgeometry *geo;
3139 geo = trd->GetGeometry();
3142 AliInfo("No TRD geometry found");
3147 AliCDBManager *man = AliCDBManager::Instance();
3149 AliInfo("Could not get CDB Manager");
3153 // Get the parameter object
3154 AliTRDCommonParam *parCom = AliTRDCommonParam::Instance();
3156 AliInfo("Could not get CommonParam");
3159 AliTRDcalibDB *cal = AliTRDcalibDB::Instance();
3161 AliInfo("Could not get calibDB");
3166 fTimeMax = cal->GetNumberOfTimeBins();
3167 fSf = (Float_t) parCom->GetSamplingFrequency();
3168 if (fRelativeScaleAuto) {
3172 if (fRelativeScale <= 0.0) {
3173 AliInfo("You have to set the relativescale factor per hand!");
3178 // Create the 2D histos corresponding to the pad group calibration mode
3181 AliInfo(Form("We will fill the CH2d histo with the pad calibration mode: Nz %d, and Nrphi %d"
3185 // Calcul the number of Xbins
3187 ModePadCalibration(2,0);
3188 ModePadFragmentation(0,2,0,0);
3189 fDetChamb2[0] = fNfragZ[0] * fNfragRphi[0];
3190 fNtotal[0] += 6 * 18 * fDetChamb2[0];
3191 ModePadCalibration(0,0);
3192 ModePadFragmentation(0,0,0,0);
3193 fDetChamb0[0] = fNfragZ[0] * fNfragRphi[0];
3194 fNtotal[0] += 6 * 4 * 18 * fDetChamb0[0];
3195 AliInfo(Form("Total number of Xbins: %d",fNtotal[0]));
3197 // Create the 2D histo
3199 CreateCH2d(fNtotal[0]);
3202 fVectorCH = new TObjArray();
3203 fPlaCH = new TObjArray();
3210 AliInfo(Form("We will fill the PH2d histo with the pad calibration mode: Nz %d, and Nrphi %d"
3214 // Calcul the number of Xbins
3216 ModePadCalibration(2,1);
3217 ModePadFragmentation(0,2,0,1);
3218 fDetChamb2[1] = fNfragZ[1] * fNfragRphi[1];
3219 fNtotal[1] += 6 * 18 * fDetChamb2[1];
3220 ModePadCalibration(0,1);
3221 ModePadFragmentation(0,0,0,1);
3222 fDetChamb0[1] = fNfragZ[1] * fNfragRphi[1];
3223 fNtotal[1] += 6 * 4 * 18 * fDetChamb0[1];
3224 AliInfo(Form("Total number of Xbins: %d",fNtotal[1]));
3226 // Create the 2D histo
3228 CreatePH2d(fNtotal[1]);
3231 fVectorPH = new TObjArray();
3232 fPlaPH = new TObjArray();
3237 loader->LoadDigits();
3238 AliInfo("LoadDigits ");
3239 AliTRDdigitsManager *digitsManager = new AliTRDdigitsManager();
3241 //iev2 egal to the max if 0
3243 iev2 = rl->GetNumberOfEvents();
3244 AliInfo(Form("Total number of events: %d",iev2));
3248 for (Int_t ievent = iev1; ievent < iev2; ievent++) {
3249 AliInfo(Form("Process event %d",ievent));
3250 rl->GetEvent(ievent);
3251 if (!loader->TreeD()) {
3252 AliInfo("loader Loading Digits ... ");
3253 loader->LoadDigits();
3255 digitsManager->ReadDigits(loader->TreeD());
3256 AliInfo("digitsManager Read Digits Done");
3257 // Read the digits from the file
3258 if (!(digitsManager->ReadDigits(loader->TreeD()))) {
3263 for (Int_t iSect = 0; iSect < 18; iSect++) {
3264 for (Int_t iChamb = 0; iChamb < kNcham; iChamb++) {
3265 for (Int_t iPlane = 0; iPlane < kNplan; iPlane++) {
3267 // A little geometry:
3268 Int_t iDet = geo->GetDetector(iPlane,iChamb,iSect);
3269 Int_t rowMax = parCom->GetRowMax(iPlane,iChamb,iSect);
3270 Int_t colMax = parCom->GetColMax(iPlane);
3272 // Variables for the group
3273 LocalisationDetectorXbins(iDet);
3275 // In the cas of charge
3277 amptotal = new Float_t[fNfragRphi[0]*fNfragZ[0]];
3279 for (Int_t k = 0; k < fNfragRphi[0]*fNfragZ[0]; k++) {
3284 // Loop through the detector pixel
3285 for (Int_t time = 0; time < fTimeMax; time++) {
3286 for (Int_t col = 0; col < colMax; col++) {
3287 for (Int_t row = 0; row < rowMax; row++) {
3289 // Amplitude and position in pad group
3290 AliTRDdigit *digit = digitsManager->GetDigit(row,col,time,iDet);
3291 Int_t amp = digit->GetAmp();
3292 Int_t posr[2] = {0,0};
3293 Int_t posc[2] = {0,0};
3296 posr[0] = (Int_t) row / fNnZ[0];
3299 (fNnRphi[0] != 0)) {
3300 posc[0] = (Int_t) col / fNnRphi[0];
3304 posr[1] = (Int_t) row / fNnZ[1];
3307 (fNnRphi[1] != 0)) {
3308 posc[1] = (Int_t) col / fNnRphi[1];
3313 if (amp < fThresholdDigit) {
3316 amptotal[(Int_t) (posc[0]*fNfragZ[0]+posr[0])] += amp;
3320 fPH2d->Fill((fXbins[1]+posc[1]*fNfragZ[1]+posr[1])+0.5,(Float_t) time/fSf,(Double_t) amp);
3323 UpdateVectorPH((fXbins[1]+posc[1]*fNfragZ[1]+posr[1]),time,(Double_t) amp);
3336 // If automatic scale
3337 if ((fCountRelativeScale < 100) && (fRelativeScaleAuto)) {
3338 // Take only the one zone track
3339 for (Int_t k = 0; k < fNfragZ[0]*fNfragRphi[0]; k++) {
3340 if ((fCountRelativeScale < 100) && (amptotal[k] > 2.0)) {
3341 fRelativeScale += amptotal[k]*0.014*0.01;
3342 fCountRelativeScale++;
3347 // We fill the CH2d after having scale with the first 100
3348 if ((fCountRelativeScale >= 100) && (fRelativeScaleAuto)) {
3350 for (Int_t k = 0; k < fNfragZ[0]*fNfragRphi[0]; k++) {
3352 (amptotal[k] > 0.0)) {
3353 fCH2d->Fill(fXbins[0]+k+0.5,amptotal[k]/fRelativeScale);
3356 (amptotal[k] > 0.0)) {
3357 UpdateVectorCH(fXbins[0]+k ,amptotal[k]/fRelativeScale);
3362 // No relative salce
3363 if (!fRelativeScaleAuto) {
3364 for (Int_t k = 0; k < fNfragZ[0]*fNfragRphi[0]; k++) {
3366 (amptotal[k] > 0.0)) {
3367 fCH2d->Fill(fXbins[0]+k+0.5, amptotal[k]/fRelativeScale);
3370 (amptotal[k] > 0.0)) {
3371 UpdateVectorCH(fXbins[0]+k, amptotal[k]/fRelativeScale);
3384 loader->UnloadDigits();
3389 if (fPH2dOn && fHisto2d) {
3392 if (fCH2dOn && fHisto2d) {
3397 if (fWrite[0] || fWrite[1]) {
3399 TFile *fout = TFile::Open(fWriteName,"RECREATE");
3400 // Check if the file could be opened
3401 if (!fout || !fout->IsOpen()) {
3402 AliInfo("<No File found!");
3406 if (fCH2dOn && fHisto2d && fWrite[0]) {
3407 fout->WriteTObject(fCH2d);
3409 if (fPH2dOn && fHisto2d && fWrite[1]) {
3410 fout->WriteTObject(fPH2d);
3413 if (fVector2d && fCH2dOn && fWrite[0]) {
3418 TTree *treeCH2d = ConvertVectorCTTreeHisto(fVectorCH,fPlaCH,"treeCH2d",(const char *) name);
3419 fout->WriteTObject(treeCH2d);
3422 if (fVector2d && fPH2dOn && fWrite[1]) {
3427 TTree *treePH2d = ConvertVectorPTreeHisto(fVectorPH,fPlaPH,"treePH2d",(const char *) name);
3428 fout->WriteTObject(treePH2d);
3439 //____________Function fill 2D all objects from Raw Data_______________________
3440 Bool_t AliTRDCalibra::Create2DRaDaOnline(Int_t iev1, Int_t iev2)
3443 // After having written the RAW DATA in the current directory, create the
3444 // 2D histos from these RAW DATA
3445 // Only for CH and PH
3448 const Int_t kNplan = 6;
3449 const Int_t kNcham = 5;
3450 TString dirname(".");
3453 AliCDBManager *man = AliCDBManager::Instance();
3455 AliInfo("Could not get CDB Manager");
3459 // Get the parameter object
3460 AliTRDCommonParam *parCom = AliTRDCommonParam::Instance();
3462 AliInfo("Could not get CommonParam");
3466 AliTRDcalibDB *cal = AliTRDcalibDB::Instance();
3468 AliInfo("Could not get calibDB");
3473 fTimeMax = cal->GetNumberOfTimeBins();
3474 fSf = (Float_t) parCom->GetSamplingFrequency();
3475 if (fRelativeScaleAuto) {
3479 if (fRelativeScale <= 0.0) {
3480 AliInfo("You have to set the relativescale factor per hand!");
3485 // Create the 2D histo corresponding to the pad group calibration mode
3488 AliInfo(Form("We will fill the CH2d histo with the pad calibration mode: Nz %d, and Nrphi %d"
3492 // Calcul the number of Xbins
3494 ModePadCalibration(2,0);
3495 ModePadFragmentation(0,2,0,0);
3496 fDetChamb2[0] = fNfragZ[0] * fNfragRphi[0];
3497 fNtotal[0] += 6 * 18 * fDetChamb2[0];
3498 ModePadCalibration(0,0);
3499 ModePadFragmentation(0,0,0,0);
3500 fDetChamb0[0] = fNfragZ[0] * fNfragRphi[0];
3501 fNtotal[0] += 6 * 4 * 18 * fDetChamb0[0];
3502 AliInfo(Form("Total number of Xbins: %d",fNtotal[0]));
3504 // Create the 2D histo
3506 CreateCH2d(fNtotal[0]);
3509 fVectorCH = new TObjArray();
3510 fPlaCH = new TObjArray();
3517 AliInfo(Form("We will fill the PH2d histo with the pad calibration mode: Nz %d, and Nrphi %d"
3521 // Calcul the number of Xbins
3523 ModePadCalibration(2,1);
3524 ModePadFragmentation(0,2,0,1);
3525 fDetChamb2[1] = fNfragZ[1] * fNfragRphi[1];
3526 fNtotal[1] += 6 * 18 * fDetChamb2[1];
3527 ModePadCalibration(0,1);
3528 ModePadFragmentation(0,0,0,1);
3529 fDetChamb0[1] = fNfragZ[1] * fNfragRphi[1];
3530 fNtotal[1] += 6 * 4 * 18 * fDetChamb0[1];
3531 AliInfo(Form("Total number of Xbins: %d",fNtotal[1]));
3533 // Create the 2D histo
3535 CreatePH2d(fNtotal[1]);
3538 fVectorPH = new TObjArray();
3539 fPlaPH = new TObjArray();
3544 AliTRDrawData *rawdata = new AliTRDrawData();
3545 AliInfo("AliTRDrawData object created ");
3548 for (Int_t ievent = iev1; ievent < iev2; ievent++) {
3551 AliRawReaderFile *readerfile = new AliRawReaderFile(dirname,ievent);
3553 AliInfo("No readerfile found!");
3557 AliTRDdigitsManager *digitsManager = rawdata->Raw2Digits((AliRawReader *) readerfile);
3558 if (!digitsManager) {
3559 AliInfo("No DigitsManager done!");
3563 // Loop on detectors
3564 for (Int_t iSect = 0; iSect < 18; iSect++) {
3565 for (Int_t iPlane = 0; iPlane < kNplan; iPlane++) {
3566 for (Int_t iChamb = 0; iChamb < kNcham; iChamb++) {
3568 // A little geometry:
3569 Int_t iDet = AliTRDgeometry::GetDetector(iPlane,iChamb,iSect);
3570 Int_t rowMax = parCom->GetRowMax(iPlane,iChamb,iSect);
3571 Int_t colMax = parCom->GetColMax(iPlane);
3573 // Variables for the group
3574 LocalisationDetectorXbins(iDet);
3576 // In the cas of charge
3578 amptotal = new Float_t[fNfragRphi[0]*fNfragZ[0]];
3580 for (Int_t k = 0; k < fNfragRphi[0]*fNfragZ[0]; k++) {
3585 // Loop through the detector pixel
3586 for (Int_t time = 0; time < fTimeMax; time++) {
3587 for (Int_t col = 0; col < colMax; col++) {
3588 for (Int_t row = 0; row < rowMax; row++) {
3590 // Amplitude and position of the digit
3591 AliTRDdigit *digit = digitsManager->GetDigit(row,col,time,iDet);
3592 Int_t amp = digit->GetAmp();
3593 Int_t posr[2] = { 0, 0 };
3594 Int_t posc[2] = { 0, 0 };
3597 posr[0] = (Int_t) row / fNnZ[0];
3600 (fNnRphi[0] != 0)) {
3601 posc[0] = (Int_t) col / fNnRphi[0];
3605 posr[1] = (Int_t) row / fNnZ[1];
3608 (fNnRphi[1] != 0)) {
3609 posc[1] = (Int_t) col / fNnRphi[1];
3614 if (amp < fThresholdDigit) {
3617 amptotal[(Int_t) (posc[0]*fNfragZ[0]+posr[0])] += amp;
3622 fPH2d->Fill((fXbins[1]+posc[1]*fNfragZ[1]+posr[1])+0.5,(Float_t)time/fSf,amp);
3625 UpdateVectorPH(fXbins[1]+posc[1]*fNfragZ[1]+posr[1],time,amp);
3637 // If automatic scale
3638 if ((fCountRelativeScale < 100) && (fRelativeScaleAuto)) {
3639 // Take only the one zone track
3640 for (Int_t k = 0; k < fNfragZ[0]*fNfragRphi[0]; k++) {
3641 if ((fCountRelativeScale < 100) && (amptotal[k] > 2.0)) {
3642 fRelativeScale += amptotal[k] * 0.014 * 0.01;
3643 fCountRelativeScale++;
3648 // We fill the CH2d after having scale with the first 100
3649 if ((fCountRelativeScale >= 100) && (fRelativeScaleAuto)) {
3651 for (Int_t k = 0; k < fNfragZ[0]*fNfragRphi[0]; k++) {
3652 if (fHisto2d && (amptotal[k] > 0.0)) {
3653 fCH2d->Fill(fXbins[0]+k+0.5,amptotal[k]/fRelativeScale);
3655 if (fVector2d && (amptotal[k] > 0.0)) {
3656 UpdateVectorCH(fXbins[0]+k, amptotal[k]/fRelativeScale);
3661 // No relative salce
3662 if (!fRelativeScaleAuto) {
3663 for (Int_t k = 0; k < fNfragZ[0]*fNfragRphi[0]; k++) {
3665 (amptotal[k] > 0.0)) {
3666 fCH2d->Fill(fXbins[0]+k+0.5,amptotal[k]/fRelativeScale);
3669 (amptotal[k] > 0.0)) {
3670 UpdateVectorCH(fXbins[0]+k, amptotal[k]/fRelativeScale);
3683 delete digitsManager;
3689 if (fPH2dOn && fHisto2d) {
3692 if (fCH2dOn && fHisto2d) {
3697 if (fWrite[0] || fWrite[1]) {
3699 TFile *fout = TFile::Open(fWriteName,"UPDATE");
3700 // Check if the file could be opened
3701 if (!fout || !fout->IsOpen()) {
3702 AliInfo("<No File found!");
3706 if (fCH2dOn && fHisto2d && fWrite[0]) {
3707 fout->WriteTObject(fCH2d);
3709 if (fPH2dOn && fHisto2d && fWrite[1]) {
3710 fout->WriteTObject(fPH2d);
3713 if (fVector2d && fCH2dOn && fWrite[0]) {
3718 TTree *treeCH2d = ConvertVectorCTTreeHisto(fVectorCH,fPlaCH,"treeCH2d",(const Char_t *) name);
3719 fout->WriteTObject(treeCH2d);
3722 if (fVector2d && fPH2dOn && fWrite[1]) {
3727 TTree *treePH2d = ConvertVectorPTreeHisto(fVectorPH,fPlaPH,"treePH2d",(const Char_t *) name);
3728 fout->WriteTObject(treePH2d);
3737 //____________Pad Calibration Public___________________________________________
3739 //____________Define the number of pads per group for one detector and one calibration
3740 void AliTRDCalibra::ModePadCalibration(Int_t iChamb, Int_t i)
3743 // Definition of the calibration mode
3744 // from Nz and Nrphi, the number of row and col pads per calibration groups are setted
3751 if ((fNz[i] == 0) && (iChamb == 2)) {
3754 if ((fNz[i] == 0) && (iChamb != 2)) {
3757 if ((fNz[i] == 1) && (iChamb == 2)) {
3760 if ((fNz[i] == 1) && (iChamb != 2)) {
3763 if ((fNz[i] == 2) && (iChamb == 2)) {
3766 if ((fNz[i] == 2) && (iChamb != 2)) {
3776 if (fNrphi[i] == 0) {
3779 if (fNrphi[i] == 1) {
3782 if (fNrphi[i] == 2) {
3785 if (fNrphi[i] == 3) {
3788 if (fNrphi[i] == 4) {
3791 if (fNrphi[i] == 5) {
3794 if (fNrphi[i] == 6) {
3800 //____________Define the number of pad groups in one detector for one calibration
3801 Bool_t AliTRDCalibra::ModePadFragmentation(Int_t iPlane,Int_t iChamb, Int_t iSect, Int_t i)
3804 // Definition of the calibration mode
3805 // From the number of row and col pads per calibration groups the
3806 // number of calibration groups are setted
3812 AliTRDCommonParam *parCom = AliTRDCommonParam::Instance();
3814 AliInfo("Could not get CommonParam Manager");
3818 // A little geometry:
3819 Int_t rowMax = parCom->GetRowMax(iPlane,iChamb,iSect);
3820 Int_t colMax = parCom->GetColMax(iPlane);
3822 // The fragmentation
3824 fNfragZ[i] = (Int_t) rowMax / fNnZ[i];
3827 if (fNnRphi[i] != 0) {
3828 fNfragRphi[i] = (Int_t) colMax / fNnRphi[i];
3835 //____________Protected Functions______________________________________________
3836 //____________Create the 2D histo to be filled online__________________________
3839 //_____________________________________________________________________________
3840 void AliTRDCalibra::CreatePRF2d(Int_t nn)
3843 // Create the 2D histos
3851 fPRF2d = new TProfile2D("PRF2d",(const Char_t *) name
3852 ,nn,0,nn,fNumberBinPRF,-1.0,1.0);
3853 fPRF2d->SetXTitle("Det/pad groups");
3854 fPRF2d->SetYTitle("Position x/W [pad width units]");
3855 fPRF2d->SetZTitle("Q_{i}/Q_{total}");
3856 fPRF2d->SetStats(0);
3860 //_____________________________________________________________________________
3861 void AliTRDCalibra::CreatePH2d(Int_t nn)
3864 // Create the 2D histos
3872 fPH2d = new TProfile2D("PH2d",(const Char_t *) name
3874 ,-0.5/fSf,(Float_t) (fTimeMax-0.5)/fSf);
3875 fPH2d->SetXTitle("Det/pad groups");
3876 fPH2d->SetYTitle("time [#mus]");
3877 fPH2d->SetZTitle("<PH> [a.u.]");
3882 //_____________________________________________________________________________
3883 void AliTRDCalibra::CreateCH2d(Int_t nn)
3886 // Create the 2D histos
3894 fCH2d = new TH2I("CH2d",(const Char_t *) name
3895 ,nn,0,nn,fNumberBinCharge,0,300);
3896 fCH2d->SetXTitle("Det/pad groups");
3897 fCH2d->SetYTitle("charge deposit [a.u]");
3898 fCH2d->SetZTitle("counts");
3904 //____________Offine tracking in the AliTRDtracker_____________________________
3905 void AliTRDCalibra::FillTheInfoOfTheTrackCH()
3908 // For the offline tracking or mcm tracklets
3909 // This function will be called in the functions UpdateHistogram...
3910 // to fill the info of a track for the relativ gain calibration
3913 Int_t nb = 0; // Nombre de zones traversees
3914 Int_t fd = -1; // Premiere zone non nulle
3917 // See if the track goes through different zones
3918 for (Int_t k = 0; k < fNfragZ[0]*fNfragRphi[0]; k++) {
3919 if (fAmpTotal[k] > 0.0) {
3927 // If automatic scale
3928 if ((fCountRelativeScale < 100) && (fRelativeScaleAuto)) {
3929 // Take only the one zone track
3931 fRelativeScale += fAmpTotal[fd] * 0.014 * 0.01;
3932 fCountRelativeScale++;
3936 // We fill the CH2d after having scale with the first 100
3937 if ((fCountRelativeScale >= 100) && (fRelativeScaleAuto)) {
3938 // Case of track with only one zone
3941 fCH2d->Fill(fXbins[0]+fd+0.5,fAmpTotal[fd]/fRelativeScale);
3944 UpdateVectorCH(fXbins[0]+fd,fAmpTotal[fd]/fRelativeScale);
3947 // Case of track with two zones
3949 // Two zones voisines sinon rien!
3950 if ((fAmpTotal[fd] > 0.0) &&
3951 (fAmpTotal[fd+1] > 0.0)) {
3952 // One of the two very big
3953 if (fAmpTotal[fd] > fProcent*fAmpTotal[fd+1]) {
3955 fCH2d->Fill(fXbins[0]+fd+0.5,fAmpTotal[fd]/fRelativeScale);
3958 UpdateVectorCH(fXbins[0]+fd,fAmpTotal[fd]/fRelativeScale);
3961 if (fAmpTotal[fd+1] > fProcent*fAmpTotal[fd]) {
3963 fCH2d->Fill(fXbins[0]+fd+1.5,fAmpTotal[fd+1]/fRelativeScale);
3966 UpdateVectorCH(fXbins[0]+fd,fAmpTotal[fd+1]/fRelativeScale);
3973 // Fill with no automatic scale
3974 if (!fRelativeScaleAuto) {
3975 // Case of track with only one zone
3979 fCH2d->Fill(fXbins[0]+fd+0.5,fAmpTotal[fd]/fRelativeScale);
3982 UpdateVectorCH(fXbins[0]+fd,fAmpTotal[fd]/fRelativeScale);
3985 // Case of track with two zones
3987 // Two zones voisines sinon rien!
3989 if ((fAmpTotal[fd] > 0.0) &&
3990 (fAmpTotal[fd+1] > 0.0)) {
3991 // One of the two very big
3992 if (fAmpTotal[fd] > fProcent*fAmpTotal[fd+1]) {
3994 fCH2d->Fill(fXbins[0]+fd+0.5,fAmpTotal[fd]/fRelativeScale);
3997 UpdateVectorCH(fXbins[0]+fd,fAmpTotal[fd]/fRelativeScale);
4001 if (fAmpTotal[fd+1] > fProcent*fAmpTotal[fd]) {
4003 fCH2d->Fill(fXbins[0]+fd+1.5,fAmpTotal[fd+1]/fRelativeScale);
4006 UpdateVectorCH(fXbins[0]+fd+1,fAmpTotal[fd+1]/fRelativeScale);
4012 if (fNfragZ[0] > 1) {
4013 if (fAmpTotal[fd] > 0.0) {
4014 if ((fd+fNfragZ[0]) < (fNfragZ[0]*fNfragRphi[0])) {
4015 if (fAmpTotal[fd+fNfragZ[0]] > 0.0) {
4016 // One of the two very big
4017 if (fAmpTotal[fd] > fProcent*fAmpTotal[fd+fNfragZ[0]]) {
4019 fCH2d->Fill(fXbins[0]+fd+0.5,fAmpTotal[fd]/fRelativeScale);
4022 UpdateVectorCH(fXbins[0]+fd,fAmpTotal[fd]/fRelativeScale);
4026 if (fAmpTotal[fd+fNfragZ[0]] > fProcent*fAmpTotal[fd]) {
4028 fCH2d->Fill(fXbins[0]+fd+fNfragZ[0]+0.5,fAmpTotal[fd+fNfragZ[0]]/fRelativeScale);
4032 UpdateVectorCH(fXbins[0]+fd+fNfragZ[0],fAmpTotal[fd+fNfragZ[0]]/fRelativeScale);
4045 //____________Offine tracking in the AliTRDtracker_____________________________
4046 void AliTRDCalibra::ResetfVariables()
4049 // Reset values of fAmpTotal, fPHValue and fPHPlace for
4050 // the updateHistogram... functions
4053 // Reset the good track
4056 // Reset the fAmpTotal where we put value
4058 for (Int_t k = 0; k < fNfragZ[0]*fNfragRphi[0]; k++) {
4063 // Reset the fPHValue
4065 for (Int_t k = 0; k < fTimeMax; k++) {
4073 //____________Offine tracking in the AliTRDtracker_____________________________
4074 void AliTRDCalibra::FillTheInfoOfTheTrackPH()
4077 // For the offline tracking or mcm tracklets
4078 // This function will be called in the functions UpdateHistogram...
4079 // to fill the info of a track for the drift velocity calibration
4082 Int_t nb = 1; // Nombre de zones traversees 1, 2 ou plus de 3
4083 Int_t fd1 = -1; // Premiere zone non nulle
4084 Int_t fd2 = -1; // Deuxieme zone non nulle
4085 Int_t k1 = -1; // Debut de la premiere zone
4086 Int_t k2 = -1; // Debut de la seconde zone
4088 // See if the track goes through different zones
4089 for (Int_t k = 0; k < fTimeMax; k++) {
4090 if (fPHValue[k] > 0.0) {
4095 if (fPHPlace[k] != fd1) {
4101 if (fPHPlace[k] != fd2) {
4109 // Case of track with only one zone
4112 for (Int_t i = 0; i < fTimeMax; i++) {
4113 if (fPHValue[i] > 0.0) {
4115 fPH2d->Fill((fXbins[1]+fPHPlace[i])+0.5,(Float_t) i/fSf,(Float_t) fPHValue[i]);
4118 AliInfo(Form("WRITE nb %d ,place final: %d, fPHPlace[i]: %d, i: %d, fPHValue[i]: %f"
4119 ,nb,fXbins[1]+fPHPlace[i],fPHPlace[i],i,fPHValue[i]));
4122 UpdateVectorPH(fXbins[1]+fPHPlace[i],i,fPHValue[i]);
4127 // Case of track with two zones
4129 // Two zones voisines sinon rien!
4131 if ((fd1 == fd2+1) ||
4133 // One of the two fast all the think
4134 if (k2 > (k1+fDifference)) {
4136 for (Int_t i = k1; i < k2; i++) {
4137 if (fPHValue[i] > 0.0) {
4139 fPH2d->Fill((fXbins[1]+fPHPlace[i])+0.5,(Float_t) i/fSf,(Float_t) fPHValue[i]);
4142 UpdateVectorPH(fXbins[1]+fPHPlace[i],i,fPHValue[i]);
4147 if ((k2+fDifference) < fTimeMax) {
4149 for (Int_t i = k2; i < fTimeMax; i++) {
4150 if (fPHValue[i] > 0.0) {
4152 fPH2d->Fill((fXbins[1]+fPHPlace[i])+0.5,(Float_t) i/fSf,(Float_t) fPHValue[i]);
4155 UpdateVectorPH(fXbins[1]+fPHPlace[i],i,fPHValue[i]);
4161 // Two zones voisines sinon rien!
4162 if (fNfragZ[1] > 1) {
4164 if ((fd1+fNfragZ[1]) < (fNfragZ[1]*fNfragRphi[1])) {
4165 if (fd2 == (fd1+fNfragZ[1])) {
4166 // One of the two fast all the think
4167 if (k2 > (k1+fDifference)) {
4169 for (Int_t i = k1; i < k2; i++) {
4170 if (fPHValue[i] > 0.0) {
4172 fPH2d->Fill((fXbins[1]+fPHPlace[i])+0.5,(Float_t) i/fSf,(Float_t) fPHValue[i]);
4175 UpdateVectorPH(fXbins[1]+fPHPlace[i],i,fPHValue[i]);
4180 if ((k2+fDifference) < fTimeMax) {
4182 for (Int_t i = k2; i < fTimeMax; i++) {
4183 if (fPHValue[i] > 0.0) {
4185 fPH2d->Fill((fXbins[1]+fPHPlace[i])+0.5,(Float_t) i/fSf,(Float_t) fPHValue[i]);
4188 UpdateVectorPH(fXbins[1]+fPHPlace[i],i,fPHValue[i]);
4195 // Two zones voisines sinon rien!
4197 if ((fd1 - fNfragZ[1]) >= 0) {
4198 if (fd2 == (fd1 - fNfragZ[1])) {
4199 // One of the two fast all the think
4200 if (k2 > (k1 + fDifference)) {
4202 for (Int_t i = k1; i < k2; i++) {
4203 if (fPHValue[i] > 0.0) {
4205 fPH2d->Fill((fXbins[1]+fPHPlace[i])+0.5,(Float_t) i/fSf,(Float_t) fPHValue[i]);
4208 UpdateVectorPH(fXbins[1]+fPHPlace[i],i,fPHValue[i]);
4213 if ((k2+fDifference) < fTimeMax) {
4215 for (Int_t i = k2; i < fTimeMax; i++) {
4216 if (fPHValue[i] > 0.0) {
4218 fPH2d->Fill((fXbins[1]+fPHPlace[i])+0.5,(Float_t) i/fSf,(Float_t) fPHValue[i]);
4221 UpdateVectorPH(fXbins[1]+fPHPlace[i],i,fPHValue[i]);
4234 //____________Set the pad calibration variables for the detector_______________
4235 Bool_t AliTRDCalibra::LocalisationDetectorXbins(Int_t detector)
4238 // For the detector calcul the first Xbins and set the number of row
4239 // and col pads per calibration groups, the number of calibration
4240 // groups in the detector.
4243 // first Xbins of the detector
4245 CalculXBins(detector,0);
4248 CalculXBins(detector,1);
4251 CalculXBins(detector,2);
4254 // fragmentation of idect
4255 for (Int_t i = 0; i < 3; i++) {
4256 ModePadCalibration((Int_t) GetChamber(detector),i);
4257 ModePadFragmentation((Int_t) GetPlane(detector)
4258 , (Int_t) GetChamber(detector)
4259 , (Int_t) GetSector(detector),i);
4266 //____________Plot the 2D histos filled Online_________________________________
4268 //_____________________________________________________________________________
4269 void AliTRDCalibra::PlotPH2d()
4272 // Plot the 2D histo
4275 TCanvas *cph2d = new TCanvas("cph2d","",50,50,600,800);
4277 fPH2d->Draw("LEGO");
4281 //_____________________________________________________________________________
4282 void AliTRDCalibra::PlotCH2d()
4285 // Plot the 2D histos
4288 TCanvas *cch2d = new TCanvas("cch2d","",50,50,600,800);
4290 fCH2d->Draw("LEGO");
4294 //_____________________________________________________________________________
4295 void AliTRDCalibra::PlotPRF2d()
4298 // Plot the 2D histos
4301 TCanvas *cPRF2d = new TCanvas("cPRF2d","",50,50,600,800);
4303 fPRF2d->Draw("LEGO");
4307 //____________Fit______________________________________________________________
4309 //____________Create histos if fDebug == 1 or fDebug >= 3______________________
4311 //_____________________________________________________________________________
4312 void AliTRDCalibra::InitArrayFitPH()
4315 // Initialise fCoefVdrift[3] and fCoefVdriftE[2] to the right dimension
4318 Int_t nbins = fDect2[1]-fDect1[1];
4320 //Init the pointer to nbins
4321 fCoefVdrift[0] = new Double_t[nbins];
4322 fCoefVdrift[1] = new Double_t[nbins];
4323 fCoefVdrift[2] = new Double_t[nbins];
4325 fCoefVdriftE[0] = new Double_t[nbins];
4326 fCoefVdriftE[1] = new Double_t[nbins];
4328 for(Int_t k = 0; k < nbins; k++){
4329 fCoefVdriftE[0][k] = 0.0;
4330 fCoefVdriftE[1][k] = 0.0;
4335 //_____________________________________________________________________________
4336 void AliTRDCalibra::InitArrayFitT0()
4339 // Initialise fCoefT0[3] and fCoefT0E[2] to the right dimension
4342 Int_t nbins = fDect2[1]-fDect1[1];
4344 //Init the pointer to nbins
4345 fCoefT0[0] = new Double_t[nbins];
4346 fCoefT0[1] = new Double_t[nbins];
4347 fCoefT0[2] = new Double_t[nbins];
4349 fCoefT0E[0] = new Double_t[nbins];
4350 fCoefT0E[1] = new Double_t[nbins];
4352 for(Int_t k = 0; k < nbins; k++){
4353 fCoefT0E[0][k] = 0.0;
4354 fCoefT0E[1][k] = 0.0;
4359 //_____________________________________________________________________________
4360 void AliTRDCalibra::InitArrayFitCH()
4363 // Initialise fCoefCharge[4] and fCoefChargeE[3] to the right dimension
4366 Int_t nbins = fDect2[0]-fDect1[0];
4368 //Init the pointer to nbins
4369 fCoefCharge[0] = new Double_t[nbins];
4370 fCoefCharge[1] = new Double_t[nbins];
4371 fCoefCharge[2] = new Double_t[nbins];
4372 fCoefCharge[3] = new Double_t[nbins];
4374 fCoefChargeE[0] = new Double_t[nbins];
4375 fCoefChargeE[1] = new Double_t[nbins];
4376 fCoefChargeE[2] = new Double_t[nbins];
4378 for(Int_t k = 0; k < nbins; k++){
4379 fCoefChargeE[0][k] = 0.0;
4380 fCoefChargeE[1][k] = 0.0;
4381 fCoefChargeE[2][k] = 0.0;
4387 //_____________________________________________________________________________
4388 void AliTRDCalibra::InitArrayFitPRF()
4391 // Initialise fCoefPRF[2] and fCoefPRFE to the right dimension
4394 Int_t nbins = fDect2[2]-fDect1[2];
4396 //Init the pointer to nbins
4397 fCoefPRF[0] = new Double_t[nbins];
4398 fCoefPRF[1] = new Double_t[nbins];
4400 fCoefPRFE = new Double_t[nbins];
4402 for(Int_t k = 0; k < nbins; k++){
4408 //_____________________________________________________________________________
4409 void AliTRDCalibra::CreateFitHistoPRFDB(Int_t rowMax, Int_t colMax)
4412 // Create the histos for fDebug = 3 and fDebug = 4 (Fit functions)
4415 fCoefPRFDB = new TH2F("coefPRF","",rowMax,0,rowMax,colMax,0,colMax);
4417 fCoefPRFDB->SetStats(0);
4418 fCoefPRFDB->SetXTitle("row Number");
4419 fCoefPRFDB->SetYTitle("col Number");
4420 fCoefPRFDB->SetZTitle("PRF width [pad width units]");
4422 fCoefPRFDB->SetFillColor(6);
4423 fCoefPRFDB->SetLineColor(6);
4427 //_____________________________________________________________________________
4428 void AliTRDCalibra::CreateFitHistoCHDB(Int_t rowMax, Int_t colMax)
4431 // Create the histos for fDebug = 3 and fDebug = 4 (Fit functions)
4434 fCoefChargeDB[0] = new TH2F("coefchargedb0","",rowMax,0,rowMax,colMax,0,colMax);
4435 fCoefChargeDB[1] = new TH2F("coefchargedb1","",rowMax,0,rowMax,colMax,0,colMax);
4436 fCoefChargeDB[2] = new TH2F("coefchargedb2","",rowMax,0,rowMax,colMax,0,colMax);
4438 fCoefChargeDB[0]->SetStats(0);
4439 fCoefChargeDB[1]->SetStats(0);
4440 fCoefChargeDB[2]->SetStats(0);
4441 fCoefChargeDB[0]->SetXTitle("row Number");
4442 fCoefChargeDB[0]->SetYTitle("col Number");
4443 fCoefChargeDB[1]->SetXTitle("row Number");
4444 fCoefChargeDB[1]->SetYTitle("col Number");
4445 fCoefChargeDB[2]->SetXTitle("row Number");
4446 fCoefChargeDB[2]->SetYTitle("col Number");
4447 fCoefChargeDB[0]->SetZTitle("f_{g} Fit method");
4448 fCoefChargeDB[1]->SetZTitle("f_{g} Mean method");
4449 fCoefChargeDB[2]->SetZTitle("f_{g} Fitbis method");
4451 fCoefChargeDB[0]->SetFillColor(6);
4452 fCoefChargeDB[0]->SetLineColor(6);
4453 fCoefChargeDB[0]->SetLineColor(6);
4454 fCoefChargeDB[1]->SetFillColor(2);
4455 fCoefChargeDB[1]->SetLineColor(2);
4456 fCoefChargeDB[1]->SetLineColor(2);
4457 fCoefChargeDB[2]->SetFillColor(8);
4458 fCoefChargeDB[2]->SetLineColor(8);
4459 fCoefChargeDB[2]->SetLineColor(8);
4463 //_____________________________________________________________________________
4464 void AliTRDCalibra::CreateFitHistoPHDB(Int_t rowMax, Int_t colMax)
4467 // Create the histos for fDebug = 3 and fDebug = 4 (Fit functions)
4470 fCoefVdriftDB[0] = new TH2F("coefvdriftdb0","",rowMax,0,rowMax,colMax,0,colMax);
4471 fCoefVdriftDB[1] = new TH2F("coefvdriftdb1","",rowMax,0,rowMax,colMax,0,colMax);
4473 fCoefVdriftDB[0]->SetStats(0);
4474 fCoefVdriftDB[1]->SetStats(0);
4475 fCoefVdriftDB[0]->SetXTitle("row Number");
4476 fCoefVdriftDB[0]->SetYTitle("col Number");
4477 fCoefVdriftDB[1]->SetXTitle("row Number");
4478 fCoefVdriftDB[1]->SetYTitle("col Number");
4479 fCoefVdriftDB[0]->SetZTitle("v_{drift} Fit method");
4480 fCoefVdriftDB[1]->SetZTitle("v_{drift} slope method");
4482 fCoefVdriftDB[0]->SetFillColor(6);
4483 fCoefVdriftDB[0]->SetLineColor(6);
4484 fCoefVdriftDB[0]->SetLineColor(6);
4485 fCoefVdriftDB[1]->SetFillColor(2);
4486 fCoefVdriftDB[1]->SetLineColor(2);
4487 fCoefVdriftDB[1]->SetLineColor(2);
4491 //_____________________________________________________________________________
4492 void AliTRDCalibra::CreateFitHistoT0DB(Int_t rowMax, Int_t colMax)
4495 // Create the histos for fDebug = 3 and fDebug = 4 (Fit functions)
4498 fCoefT0DB[0] = new TH2F("coefT0db0","",rowMax,0,rowMax,colMax,0,colMax);
4499 fCoefT0DB[1] = new TH2F("coefT0db1","",rowMax,0,rowMax,colMax,0,colMax);
4501 fCoefT0DB[0]->SetStats(0);
4502 fCoefT0DB[1]->SetStats(0);
4503 fCoefT0DB[0]->SetXTitle("row Number");
4504 fCoefT0DB[0]->SetYTitle("col Number");
4505 fCoefT0DB[1]->SetXTitle("row Number");
4506 fCoefT0DB[1]->SetYTitle("col Number");
4507 fCoefT0DB[0]->SetZTitle("t0 Fit method");
4508 fCoefT0DB[1]->SetZTitle("t0 slope method");
4510 fCoefT0DB[0]->SetFillColor(6);
4511 fCoefT0DB[0]->SetLineColor(6);
4512 fCoefT0DB[0]->SetLineColor(6);
4513 fCoefT0DB[1]->SetFillColor(2);
4514 fCoefT0DB[1]->SetLineColor(2);
4515 fCoefT0DB[1]->SetLineColor(2);
4519 //_____________________________________________________________________________
4520 Bool_t AliTRDCalibra::FillVectorFitCH(Int_t countdet)
4523 // For the Fit functions fill the vector FitCH special for the gain calibration
4526 AliTRDFitCHInfo *fitCHInfo = new AliTRDFitCHInfo();
4529 if (GetChamber(countdet) == 2) {
4536 Float_t *coef = new Float_t[ntotal];
4537 for (Int_t i = 0; i < ntotal; i++) {
4538 coef[i] = fCoefCH[i];
4541 Int_t detector = countdet;
4543 fitCHInfo->SetCoef(coef);
4544 fitCHInfo->SetDetector(detector);
4545 fVectorFitCH->Add((TObject *) fitCHInfo);
4551 //____________Functions for initialising the AliTRDCalibra in the code_________
4552 Bool_t AliTRDCalibra::InitFit(Int_t nbins, Int_t i)
4555 // Init the calibration mode (Nz, Nrphi), the histograms for
4556 // debugging the fit methods if fDebug > 0,
4559 gStyle->SetPalette(1);
4560 gStyle->SetOptStat(1111);
4561 gStyle->SetPadBorderMode(0);
4562 gStyle->SetCanvasColor(10);
4563 gStyle->SetPadLeftMargin(0.13);
4564 gStyle->SetPadRightMargin(0.01);
4566 // Get the parameter object
4567 AliTRDCommonParam *parCom = AliTRDCommonParam::Instance();
4569 AliInfo("Could not get CommonParam");
4573 // Mode groups of pads: the total number of bins!
4574 Int_t numberofbinsexpected = 0;
4575 ModePadCalibration(2,i);
4576 ModePadFragmentation(0,2,0,i);
4577 fDetChamb2[i] = fNfragZ[i] * fNfragRphi[i];
4579 AliInfo(Form("For the chamber 2: %d",fDetChamb2[i]));
4581 numberofbinsexpected += 6 * 18 * fDetChamb2[i];
4582 ModePadCalibration(0,i);
4583 ModePadFragmentation(0,0,0,i);
4584 fDetChamb0[i] = fNfragZ[i] * fNfragRphi[i];
4586 AliInfo(Form("For the other chamber 0: %d",fDetChamb0[i]));
4588 numberofbinsexpected += 6 * 4 * 18 * fDetChamb0[i];
4590 // Quick verification that we have the good pad calibration mode if 2D histos!
4592 if (numberofbinsexpected != nbins) {
4593 AliInfo("It doesn't correspond to the mode of pad group calibration!");
4598 // Security for fDebug 3 and 4
4599 if ((fDebug >= 3) &&
4603 AliInfo("This detector doesn't exit!");
4607 // Determine fDet1 and fDet2
4611 fDect1[i] = fFitVoir;
4612 fDect2[i] = fDect1[i] +1;
4616 fDect2[i] = numberofbinsexpected;
4619 CalculXBins(AliTRDgeometry::GetDetector(fDet[0],fDet[1],fDet[2]),i);
4620 fDect1[i] = fXbins[i];
4621 CalculXBins((AliTRDgeometry::GetDetector(fDet[0],fDet[1],fDet[2])+1),i);
4622 fDect2[i] = fXbins[i];
4625 // Create the histos for debugging
4630 // Init the VectorFitCH
4631 fVectorFitCH = new TObjArray();
4632 fCoefCH = new Float_t[2304];
4633 for (Int_t k = 0; k < 2304; k++) {
4636 fScaleFitFactor = 0.0;
4638 // Number of Xbins(detectors or groups of pads) if Vector2d
4639 // Quick verification that we are not out of range!
4640 if (fVectorCH && fPlaCH) {
4642 (fVectorCH->GetEntriesFast() > 0) &&
4643 ((Int_t) fPlaCH->GetEntriesFast() > 0)) {
4644 if ((Int_t) fVectorCH->GetEntriesFast() > numberofbinsexpected) {
4645 AliInfo("ch doesn't correspond to the mode of pad group calibration!");
4648 if ((Int_t) fVectorCH->GetEntriesFast() != (Int_t) fPlaCH->GetEntriesFast()) {
4649 AliInfo("VectorCH doesn't correspond to PlaCH!");
4656 // Debugging: Create the histos
4659 // fDebug == 0 nothing
4666 // fDebug == 2 and fFitVoir no histo
4668 if (fFitVoir < numberofbinsexpected) {
4669 AliInfo(Form("We will see the fit of the object %d",fFitVoir));
4672 AliInfo("fFitVoir is out of range of the histo!");
4677 // fDebug == 3 or 4 and fDet
4679 if ((fNz[0] == 0) && (fNrphi[0] == 0)) {
4680 AliInfo("Do you really want to see one detector without pad groups?");
4684 AliInfo(Form("You will see the detector: iPlane %d, iChamb %d, iSect %d"
4685 ,fDet[0],fDet[1],fDet[2]));
4686 // A little geometry:
4687 Int_t rowMax = parCom->GetRowMax(fDet[0],fDet[1],fDet[2]);
4688 Int_t colMax = parCom->GetColMax(fDet[0]);
4689 // Create the histos to visualise
4690 CreateFitHistoCHDB(rowMax,colMax);
4702 // Number of Xbins (detectors or groups of pads) if vector2d
4703 // Quick verification that we are not out of range!
4704 if (fVectorPH && fPlaPH) {
4706 (fVectorPH->GetEntriesFast() > 0) &&
4707 ((Int_t) fPlaPH->GetEntriesFast() > 0)) {
4708 if ((Int_t) fVectorPH->GetEntriesFast() > numberofbinsexpected) {
4709 AliInfo("ph doesn't correspond to the mode of pad group calibration!");
4712 if ((Int_t) fVectorPH->GetEntriesFast() != (Int_t) fPlaPH->GetEntriesFast()) {
4713 AliInfo("VectorPH doesn't correspond to PlaPH!");
4724 // Debugging: Create the histos
4727 // fDebug == 0 nothing
4731 // Create the histos replique de ph
4736 // fDebug == 2 and fFitVoir no histo
4738 if (fFitVoir < numberofbinsexpected) {
4739 AliInfo(Form("We will see the fit of the object %d",fFitVoir));
4742 AliInfo("fFitVoir is out of range of the histo!");
4747 // fDebug == 3 or 4 and fDet
4749 if ((fNz[1] == 0) &&
4751 AliInfo("Do you really want to see one detector without pad groups?");
4755 AliInfo(Form("You will see the detector: iPlane %d, iChamb %d, iSect %d"
4756 ,fDet[0],fDet[1],fDet[2]));
4757 // A little geometry:
4758 Int_t rowMax = parCom->GetRowMax(fDet[0],fDet[1],fDet[2]);
4759 Int_t colMax = parCom->GetColMax(fDet[0]);
4760 // Create the histos to visualise
4761 CreateFitHistoPHDB(rowMax,colMax);
4762 CreateFitHistoT0DB(rowMax,colMax);
4775 // Number of Xbins(detectors or groups of pads) if vector2d
4776 if (fVectorPRF && fPlaPRF){
4778 (fVectorPRF->GetEntriesFast() > 0) &&
4779 (fPlaPRF->GetEntriesFast() > 0)) {
4780 // Quick verification that we are not out of range!
4781 if ((Int_t) fVectorPRF->GetEntriesFast() > numberofbinsexpected) {
4782 AliInfo("ch doesn't correspond to the mode of pad group calibration!");
4785 if ((Int_t) fVectorPRF->GetEntriesFast() != (Int_t) fPlaPRF->GetEntriesFast()) {
4786 AliInfo("VectorPRF doesn't correspond to PlaCH!");
4796 // Debugging: Create the histos
4799 // fDebug == 0 nothing
4803 // Create the histos replique de ch
4807 // fDebug == 2 and fFitVoir no histo
4809 if (fFitVoir < numberofbinsexpected) {
4810 AliInfo(Form("We will see the fit of the object %d",fFitVoir));
4813 AliInfo("fFitVoir is out of range of the histo!");
4818 // fDebug == 3 or 4 and fDet
4820 if ((fNz[2] == 0) &&
4822 AliInfo("Do you really want to see one detector without pad groups?");
4826 AliInfo(Form("You will see the detector: iPlane %d, iChamb %d, iSect %d"
4827 ,fDet[0],fDet[1],fDet[2]));
4828 // A little geometry:
4829 Int_t rowMax = parCom->GetRowMax(fDet[0],fDet[1],fDet[2]);
4830 Int_t colMax = parCom->GetColMax(fDet[0]);
4831 // Create the histos to visualise
4832 CreateFitHistoPRFDB(rowMax,colMax);
4845 //____________Functions for initialising the AliTRDCalibra in the code_________
4846 void AliTRDCalibra::InitfCountDetAndfCount(Int_t i)
4849 // Init the current detector where we are fCountDet and the
4850 // next fCount for the functions Fit...
4853 // Loop on the Xbins of ch!!
4854 fCountDet[i] = -1; // Current detector
4855 fCount[i] = 0; // To find the next detector
4860 // Set countdet to the detector
4861 fCountDet[i] = AliTRDgeometry::GetDetector(fDet[0],fDet[1],fDet[2]);
4863 // Determination of fNnZ, fNnRphi, fNfragZ and fNfragRphi
4864 ModePadCalibration(fDet[1],i);
4865 ModePadFragmentation(fDet[0],fDet[1],fDet[2],i);
4867 // Set counter to write at the end of the detector
4868 fCount[i] = fDect1[i] + fNfragZ[i]*fNfragRphi[i];
4874 //____________Functions for initialising the AliTRDCalibra in the code_________
4875 void AliTRDCalibra::UpdatefCountDetAndfCount(Int_t idect, Int_t i)
4878 // See if we are in a new detector and update the
4879 // variables fNfragZ and fNfragRphi if yes
4882 // Determination of fNnZ, fNnRphi, fNfragZ and fNfragRphi
4883 // If fDebug == 1 or 0
4884 if ((fDebug == 0) ||
4887 if (fCount[i] == idect) {
4889 // On en est au detector
4892 // Determination of fNnZ, fNnRphi, fNfragZ and fNfragRphi
4893 ModePadCalibration((Int_t) GetChamber(fCountDet[i]),i);
4894 ModePadFragmentation((Int_t) GetPlane(fCountDet[i])
4895 ,(Int_t) GetChamber(fCountDet[i])
4896 ,(Int_t) GetSector(fCountDet[i]),i);
4898 // Set for the next detector
4899 fCount[i] += fNfragZ[i]*fNfragRphi[i];
4907 //____________Functions for initialising the AliTRDCalibra in the code_________
4908 void AliTRDCalibra::ReconstructFitRowMinRowMax(Int_t idect, Int_t i)
4911 // Reconstruct the min pad row, max pad row, min pad col and
4912 // max pad col of the calibration group for the Fit functions
4916 ReconstructionRowPadGroup((Int_t) (idect-(fCount[i]-(fNfragZ[i]*fNfragRphi[i]))),i);
4919 ReconstructionRowPadGroup((Int_t) (idect-fDect1[i]),i);
4924 //____________Functions for initialising the AliTRDCalibra in the code_________
4925 Bool_t AliTRDCalibra::NotEnoughStatistic(Int_t idect, Int_t i)
4928 // For the case where there are not enough entries in the histograms
4929 // of the calibration group, the value present in the choosen database
4930 // will be put. A negativ sign enables to know that a fit was not possible.
4933 // Get the parameter object
4934 AliTRDCommonParam *parCom = AliTRDCommonParam::Instance();
4936 AliInfo("Could not get CommonParam Manager");
4941 AliTRDcalibDB *cal = AliTRDcalibDB::Instance();
4943 AliInfo("Could not get calibDB");
4948 AliInfo(Form("The element %d in this detector %d has not enough statistic to be fitted"
4949 ,idect-(fCount[i]-(fNfragZ[i]*fNfragRphi[i])),fCountDet[i]));
4952 AliInfo("The element has not enough statistic to be fitted");
4955 if ((i == 0) && (fDebug != 2)) {
4957 // Calcul the coef from the database choosen
4958 CalculChargeCoefMean(fCountDet[0],(Int_t) (idect-fDect1[0]),kFALSE);
4960 // Fill the coefCH[2304] with negative value to say: not fitted
4961 for (Int_t k = fRowMin[0]; k < fRowMax[0]; k++) {
4962 for (Int_t j = fColMin[0]; j < fColMax[0]; j++) {
4963 if (GetChamber(fCountDet[0]) == 2) {
4964 fCoefCH[(Int_t)(j*12+k)] = -TMath::Abs(fChargeCoef[3]);
4966 if (GetChamber(fCountDet[0]) != 2) {
4967 fCoefCH[(Int_t)(j*16+k)] = -TMath::Abs(fChargeCoef[3]);
4972 // Put the default value negative
4973 if ((fDebug == 1) ||
4976 if (fFitChargeBisOn) {
4977 fCoefCharge[2][idect-fDect1[0]]=-TMath::Abs(fChargeCoef[3]);
4979 if (fMeanChargeOn) {
4980 fCoefCharge[1][idect-fDect1[0]]=-TMath::Abs(fChargeCoef[3]);
4983 fCoefCharge[0][idect-fDect1[0]]=-TMath::Abs(fChargeCoef[3]);
4987 // End of one detector
4988 if ((idect == (fCount[0]-1))) {
4989 FillVectorFitCH((Int_t) fCountDet[0]);
4991 for (Int_t k = 0; k < 2304; k++) {
4998 if ((i == 1) && (fDebug != 2)) {
5000 CalculVdriftCoefMean(fCountDet[1],(Int_t) (idect-fDect1[1]));
5001 CalculT0CoefMean(fCountDet[1],(Int_t) (idect-fDect1[1]));
5003 // Put the default value (time0 can be negativ, so we stay with + )
5004 if ((fDebug == 1) ||
5008 fCoefVdrift[0][(idect-fDect1[1])] = -fVdriftCoef[2];
5009 fCoefT0[0][(idect-fDect1[1])] = fT0Coef[2];
5012 fCoefVdrift[1][(idect-fDect1[1])] = -fVdriftCoef[2];
5013 fCoefT0[1][(idect-fDect1[1])] = fT0Coef[2];
5017 // Put the default value
5019 fVdriftCoef[0] = fVdriftCoef[2];
5020 fVdriftCoef[1] = fVdriftCoef[2];
5022 fT0Coef[0] = fT0Coef[2];
5023 fT0Coef[1] = fT0Coef[2];
5027 // Fill the tree if end of a detector.
5028 // The pointer to the branch stays with the default value negative!!!
5030 // Pointer to the branch
5031 for (Int_t k = fRowMin[1]; k < fRowMax[1]; k++) {
5032 for (Int_t j = fColMin[1]; j < fColMax[1]; j++) {
5033 if (GetChamber(fCountDet[1]) == 2) {
5034 fVdriftPad[(Int_t)(j*12+k)] = -TMath::Abs(fVdriftCoef[2]);
5036 if (GetChamber(fCountDet[1]) != 2) {
5037 fVdriftPad[(Int_t)(j*16+k)] = -TMath::Abs(fVdriftCoef[2]);
5042 // End of one detector
5043 if ((idect == (fCount[1]-1)) && (fDebug != 2)) {
5044 FillTreeVdrift((Int_t) fCountDet[1]);
5048 // Fill the tree if end of a detector.
5049 // The pointer to the branch stays with the default value positive!!!
5050 // Pointer to the branch
5051 for (Int_t k = fRowMin[1]; k < fRowMax[1]; k++) {
5052 for (Int_t j = fColMin[1]; j < fColMax[1]; j++) {
5053 if (GetChamber(fCountDet[1]) == 2) {
5054 fT0Pad[(Int_t)(j*12+k)] = fT0Coef[2];
5056 if (GetChamber(fCountDet[1]) != 2) {
5057 fT0Pad[(Int_t)(j*16+k)] = fT0Coef[2];
5062 // End of one detector
5063 if ((idect == (fCount[1]-1)) && (fDebug != 2)) {
5064 FillTreeT0((Int_t) fCountDet[1]);
5069 if ((i == 2) && (fDebug != 2)) {
5071 CalculPRFCoefMean(fCountDet[2],(Int_t) (idect-fDect1[2]));
5073 if ((fDebug == 1) ||
5075 fCoefPRF[0][(idect-fDect1[2])] = -fPRFCoef[1];
5079 fPRFCoef[0] = fPRFCoef[1];
5083 // Fill the tree if end of a detector.
5084 // The pointer to the branch stays with the default value 1.5!!!
5085 // Pointer to the branch
5086 for (Int_t k = fRowMin[2]; k < fRowMax[2]; k++) {
5087 for (Int_t j = fColMin[2]; j < fColMax[2]; j++) {
5088 if((parCom->GetColMax(GetPlane(fCountDet[2])) != (j+1)) && (j != 0)){
5089 if (GetChamber(fCountDet[2]) == 2) {
5090 fPRFPad[(Int_t)(j*12+k)] = -fPRFCoef[1];
5092 if (GetChamber(fCountDet[2]) != 2) {
5093 fPRFPad[(Int_t)(j*16+k)] = -fPRFCoef[1];
5098 if (GetChamber(fCountDet[2]) == 2) {
5099 fPRFPad[(Int_t)(j*12+k)] = -((Float_t) cal->GetPRFWidth(fCountDet[2],j,k));
5101 if (GetChamber(fCountDet[2]) != 2) {
5102 fPRFPad[(Int_t)(j*16+k)] = -((Float_t) cal->GetPRFWidth(fCountDet[2],j,k));
5106 if (GetChamber(fCountDet[2]) == 2) {
5107 fPRFPad[(Int_t)(j*12+k)] = -((Float_t) GetPRFDefault(GetPlane(fCountDet[2])));
5109 if (GetChamber(fCountDet[2]) != 2) {
5110 fPRFPad[(Int_t)(j*16+k)] = -((Float_t) GetPRFDefault(GetPlane(fCountDet[2])));
5117 // End of one detector
5118 if ((idect == (fCount[2]-1)) && (fDebug != 2)) {
5119 FillTreePRF((Int_t) fCountDet[2]);
5128 //____________Functions for initialising the AliTRDCalibra in the code_________
5129 Bool_t AliTRDCalibra::FillInfosFit(Int_t idect, Int_t i)
5132 // Fill the coefficients found with the fits or other
5133 // methods from the Fit functions
5136 // Get the parameter object
5137 AliTRDCommonParam *parCom = AliTRDCommonParam::Instance();
5139 AliInfo("Could not get CommonParam Manager");
5144 AliTRDcalibDB *cal = AliTRDcalibDB::Instance();
5146 AliInfo("Could not get calibDB");
5150 if ((i == 0) && (fDebug != 2)) {
5151 // Fill the coefCH[2304] with fChargeCoef[0]
5152 // that would be negativ only if the fit failed totally
5153 for (Int_t k = fRowMin[0]; k < fRowMax[0]; k++) {
5154 for (Int_t j = fColMin[0]; j < fColMax[0]; j++) {
5155 if (GetChamber(fCountDet[0]) == 2) {
5156 fCoefCH[(Int_t)(j*12+k)] = fChargeCoef[0];
5158 if (GetChamber(fCountDet[0]) != 2) {
5159 fCoefCH[(Int_t)(j*16+k)] = fChargeCoef[0];
5163 // End of one detector
5164 if ((idect == (fCount[0]-1))) {
5165 FillVectorFitCH((Int_t) fCountDet[0]);
5167 for (Int_t k = 0; k < 2304; k++) {
5173 if ((i == 1) && (fDebug != 2)) {
5176 // Pointer to the branch: fVdriftCoef[1] will ne negativ only if the fit failed totally
5177 for (Int_t k = fRowMin[1]; k < fRowMax[1]; k++) {
5178 for (Int_t j = fColMin[1]; j < fColMax[1]; j++) {
5179 if (GetChamber(fCountDet[1]) == 2) {
5180 fVdriftPad[(Int_t)(j*12+k)]=fVdriftCoef[1];
5182 if (GetChamber(fCountDet[1]) != 2) {
5183 fVdriftPad[(Int_t)(j*16+k)]=fVdriftCoef[1];
5187 // End of one detector
5188 if ((idect == (fCount[1]-1)) && (fDebug != 2)) {
5189 FillTreeVdrift((Int_t) fCountDet[1]);
5193 // Pointer to the branch: fT0Coef[1] will ne negativ only if the fit failed totally
5194 for (Int_t k = fRowMin[1]; k < fRowMax[1]; k++) {
5195 for (Int_t j = fColMin[1]; j < fColMax[1]; j++) {
5196 if (GetChamber(fCountDet[1]) == 2) {
5197 fT0Pad[(Int_t)(j*12+k)]=fT0Coef[1];
5199 if (GetChamber(fCountDet[1]) != 2) {
5200 fT0Pad[(Int_t)(j*16+k)]=fT0Coef[1];
5204 // End of one detector
5205 if ((idect == (fCount[1]-1)) && (fDebug != 2)) {
5206 FillTreeT0((Int_t) fCountDet[1]);
5211 if ((i == 2) && (fDebug != 2)) {
5212 // Pointer to the branch
5213 for (Int_t k = fRowMin[2]; k < fRowMax[2]; k++) {
5214 for (Int_t j = fColMin[2]; j < fColMax[2]; j++) {
5215 if ((parCom->GetColMax(GetPlane(fCountDet[2])) != (j+1)) && (j != 0)) {
5216 if (GetChamber(fCountDet[2]) == 2) {
5217 fPRFPad[(Int_t)(j*12+k)] = fPRFCoef[0];
5219 if (GetChamber(fCountDet[2]) != 2) {
5220 fPRFPad[(Int_t)(j*16+k)] = fPRFCoef[0];
5225 if (GetChamber(fCountDet[2]) == 2) {
5226 fPRFPad[(Int_t)(j*12+k)] = (Float_t) cal->GetPRFWidth(fCountDet[2],j,k);
5228 if (GetChamber(fCountDet[2]) != 2) {
5229 fPRFPad[(Int_t)(j*16+k)] = (Float_t) cal->GetPRFWidth(fCountDet[2],j,k);
5233 if (GetChamber(fCountDet[2]) == 2) {
5234 fPRFPad[(Int_t)(j*12+k)] = (Float_t) GetPRFDefault(GetPlane(fCountDet[2]));
5236 if (GetChamber(fCountDet[2]) != 2) {
5237 fPRFPad[(Int_t)(j*16+k)] = (Float_t) GetPRFDefault(GetPlane(fCountDet[2]));
5243 // End of one detector
5244 if ((idect == (fCount[2]-1)) && (fDebug != 2)) {
5245 FillTreePRF((Int_t) fCountDet[2]);
5253 //____________Functions for initialising the AliTRDCalibra in the code_________
5254 Bool_t AliTRDCalibra::WriteFitInfos(Int_t i)
5257 // In the case the user wants to write a file with a tree of the found
5258 // coefficients for the calibration before putting them in the database
5261 TFile *fout = TFile::Open(fWriteNameCoef,"UPDATE");
5262 // Check if the file could be opened
5263 if (!fout || !fout->IsOpen()) {
5264 AliInfo("No File found!");
5268 if ((i == 0) && (fDebug != 2)) {
5270 if ((fDebug == 4) ||
5275 fout->WriteTObject(fGain,fGain->GetName(),(Option_t *) "writedelete");
5278 if ((i == 1) && (fDebug != 2)) {
5280 if ((fDebug == 4) ||
5285 fout->WriteTObject(fVdrift,fVdrift->GetName(),(Option_t *) "writedelete");
5287 if ((fDebug == 4) ||
5292 fout->WriteTObject(fT0,fT0->GetName(),(Option_t *) "writedelete");
5295 if ((i == 2) && (fDebug != 2)) {
5297 if ((fDebug == 4) ||
5302 fout->WriteTObject(fPRF,fPRF->GetName(),(Option_t *) "writedelete");
5312 //____________Fill Coef DB in case of visualisation of one detector____________
5315 //_____________________________________________________________________________
5316 void AliTRDCalibra::FillCoefVdriftDB()
5319 // Fill the histos for fDebug = 3 and fDebug = 4 to visualise the detector
5322 for (Int_t row = fRowMin[1]; row < fRowMax[1]; row++) {
5323 for (Int_t col = fColMin[1]; col < fColMax[1]; col++) {
5324 fCoefVdriftDB[1]->SetBinContent(row+1,col+1,TMath::Abs(fVdriftCoef[1]));
5326 fCoefVdriftDB[0]->SetBinContent(row+1,col+1,TMath::Abs(fVdriftCoef[0]));
5333 //_____________________________________________________________________________
5334 void AliTRDCalibra::FillCoefT0DB()
5337 // Fill the histos for fDebug = 3 and fDebug = 4 to visualise the detector
5340 for (Int_t row = fRowMin[1]; row < fRowMax[1]; row++) {
5341 for (Int_t col = fColMin[1]; col < fColMax[1]; col++) {
5342 fCoefT0DB[1]->SetBinContent(row+1,col+1,TMath::Abs(fT0Coef[1]));
5344 fCoefT0DB[0]->SetBinContent(row+1,col+1,TMath::Abs(fT0Coef[0]));
5351 //_____________________________________________________________________________
5352 void AliTRDCalibra::FillCoefChargeDB()
5355 // Fill the histos for fDebug = 3 and fDebug = 4 to visualise the detector
5358 for (Int_t row = fRowMin[0]; row < fRowMax[0]; row++) {
5359 for (Int_t col = fColMin[0]; col < fColMax[0]; col++) {
5360 if (fMeanChargeOn) {
5361 fCoefChargeDB[1]->SetBinContent(row+1,col+1,TMath::Abs(fChargeCoef[1]));
5363 if (fFitChargeBisOn) {
5364 fCoefChargeDB[2]->SetBinContent(row+1,col+1,TMath::Abs(fChargeCoef[2]));
5366 fCoefChargeDB[0]->SetBinContent(row+1,col+1,TMath::Abs(fChargeCoef[0]));
5372 //_____________________________________________________________________________
5373 void AliTRDCalibra::FillCoefPRFDB()
5376 // Fill the histos for fDebug = 3 and fDebug = 4 to visualise the detector
5379 for (Int_t row = fRowMin[2]; row < fRowMax[2]; row++) {
5380 for (Int_t col = fColMin[2]; col < fColMax[2]; col++) {
5381 fCoefPRFDB->SetBinContent(row+1,col+1,fPRFCoef[0]);
5388 //____________Plot histos CoefPRF....__________________________________________
5391 //_____________________________________________________________________________
5392 void AliTRDCalibra::PlotWriteCH()
5395 // Scale the coefficients to one, create the graph errors and write them if wanted
5398 //TObjArray of the grapherrors and so on
5399 TObjArray *listofgraphs = new TObjArray();
5401 Int_t nbins = fDect2[0]-fDect1[0];
5411 Double_t scale = 1.0;
5414 Double_t *xValuesFitted = new Double_t[nbins];
5415 Double_t *xValuesFittedMean = new Double_t[nbins];
5416 Double_t *xValuesFittedBis = new Double_t[nbins];
5417 for(Int_t k = 0; k < nbins; k ++){
5418 xValuesFitted[k] = -1;
5419 xValuesFittedMean[k] = -1;
5420 xValuesFittedBis[k] = -1;
5423 for(Int_t l = 0; l < nbins; l++){
5424 if(fCoefCharge[0][l] > 0){
5425 fCoefCharge[0][l]=fCoefCharge[0][l]*fScaleFitFactor;
5426 fCoefChargeE[0][l]=fCoefChargeE[0][l]*fScaleFitFactor;
5427 xValuesFitted[counter[0]]=l;
5433 for(Int_t l = 0; l < nbins; l++){
5434 if(fCoefCharge[1][l] > 0){
5435 sum += fCoefCharge[1][l];
5436 xValuesFittedMean[counter[1]]= l;
5441 if(sum > 0.0) scale = counter[1]/sum;
5442 for(Int_t l = 0; l < nbins; l++){
5443 if(fCoefCharge[1][l] > 0){
5444 fCoefCharge[1][l]=fCoefCharge[1][l]*scale;
5445 fCoefChargeE[1][l]=fCoefChargeE[1][l]*scale;
5449 if(fFitChargeBisOn){
5451 for(Int_t l = 0; l < nbins; l++){
5452 if(fCoefCharge[2][l] > 0){
5453 fCoefCharge[2][l]=fCoefCharge[2][l]*fScaleFitFactor;
5454 fCoefChargeE[2][l]=fCoefChargeE[2][l]*fScaleFitFactor;
5455 sum += fCoefCharge[2][l];
5456 xValuesFittedBis[counter[2]]= l;
5461 if(sum > 0.0) scale = counter[2]/sum;
5462 for(Int_t l = 0; l < nbins; l++){
5463 if(fCoefCharge[2][l] > 0){
5464 fCoefCharge[2][l]=fCoefCharge[2][l]*scale;
5465 fCoefChargeE[2][l]=fCoefChargeE[2][l]*scale;
5466 if(fCoefCharge[2][l] > 1.0) printf("for the group %d, I have the coef %f with the error %f\n",l,fCoefCharge[2][l],fCoefChargeE[2][l]);
5471 //Create the X and Xerror
5472 Double_t *xValues = new Double_t[nbins];
5473 Double_t *xValuesE = new Double_t[nbins];
5474 for(Int_t k = 0; k < nbins; k ++){
5479 //Create the graph erros and plot them
5480 TGraph *graphCharge3 = new TGraph(nbins,xValues,fCoefCharge[3]);
5481 graphCharge3->SetName("coefcharge3");
5482 graphCharge3->SetTitle("");
5483 graphCharge3->GetXaxis()->SetTitle("Det/Pad groups");
5484 graphCharge3->GetYaxis()->SetTitle("gain factor");
5485 graphCharge3->SetLineColor(4);
5486 listofgraphs->Add((TObject *)graphCharge3);
5487 TGraphErrors *graphCharge0 = new TGraphErrors(nbins,xValues,fCoefCharge[0],xValuesE,fCoefChargeE[0]);
5488 graphCharge0->SetName("coefcharge0");
5489 graphCharge0->SetTitle("");
5490 graphCharge0->GetXaxis()->SetTitle("Det/Pad groups");
5491 graphCharge0->GetYaxis()->SetTitle("gain factor");
5492 graphCharge0->SetMarkerColor(6);
5493 graphCharge0->SetLineColor(6);
5494 graphCharge0->SetMarkerStyle(26);
5495 listofgraphs->Add((TObject *)graphCharge0);
5496 TCanvas *cch1 = new TCanvas("cch1","",50,50,600,800);
5498 TLegend *legch1 = new TLegend(0.4,0.6,0.89,0.89);
5499 legch1->AddEntry(graphCharge3,"f_{g} simulated","l");
5500 legch1->AddEntry(graphCharge0,"f_{g} fit","p");
5501 graphCharge0->Draw("AP");
5502 //graphCharge3->Draw("AL");
5503 if (fMeanChargeOn) {
5504 TGraphErrors *graphCharge1 = new TGraphErrors(nbins,xValues,fCoefCharge[1],xValuesE,fCoefChargeE[1]);
5505 graphCharge1->SetName("coefcharge1");
5506 graphCharge1->GetXaxis()->SetTitle("Det/Pad groups");
5507 graphCharge1->GetYaxis()->SetTitle("gain factor");
5508 graphCharge1->SetTitle("");
5509 graphCharge1->SetMarkerColor(2);
5510 graphCharge1->SetLineColor(2);
5511 graphCharge1->SetMarkerStyle(24);
5512 legch1->AddEntry(graphCharge1,"f_{g} mean","p");
5513 graphCharge1->Draw("P");
5514 listofgraphs->Add((TObject *)graphCharge1);
5516 if (fFitChargeBisOn ) {
5517 TGraphErrors *graphCharge2 = new TGraphErrors(nbins,xValues,fCoefCharge[2],xValuesE,fCoefChargeE[2]);
5518 graphCharge2->SetName("coefcharge2");
5519 graphCharge2->SetTitle("");
5520 graphCharge2->GetXaxis()->SetTitle("Det/Pad groups");
5521 graphCharge2->GetYaxis()->SetTitle("gain factor");
5522 graphCharge2->SetMarkerColor(8);
5523 graphCharge2->SetLineColor(8);
5524 graphCharge2->SetMarkerStyle(25);
5525 legch1->AddEntry(graphCharge2,"f_{g} fitbis","p");
5526 graphCharge2->Draw("P");
5527 listofgraphs->Add((TObject *)graphCharge2);
5529 legch1->Draw("same");
5532 //Create the arrays and the graphs for the delta
5533 TCanvas *cch2 = new TCanvas("cch2","",50,50,600,800);
5536 Double_t *yValuesDelta = new Double_t[counter[0]];
5537 for(Int_t k = 0; k < counter[0]; k++){
5538 if(fCoefCharge[3][(Int_t)(xValuesFitted[k])] > 0.0) {
5539 yValuesDelta[k] = (fCoefCharge[0][(Int_t)xValuesFitted[k]]-fCoefCharge[3][(Int_t)xValuesFitted[k]])/fCoefCharge[3][(Int_t)xValuesFitted[k]];
5541 else yValuesDelta[k] = 0.0;
5543 TGraph *graphDeltaCharge0 = new TGraph(counter[0],&xValuesFitted[0],yValuesDelta);
5544 graphDeltaCharge0->SetName("deltacharge0");
5545 graphDeltaCharge0->GetXaxis()->SetTitle("Det/Pad groups");
5546 graphDeltaCharge0->GetYaxis()->SetTitle("#Deltag/g_{sim}");
5547 graphDeltaCharge0->SetMarkerColor(6);
5548 graphDeltaCharge0->SetTitle("");
5549 graphDeltaCharge0->SetLineColor(6);
5550 graphDeltaCharge0->SetMarkerStyle(26);
5551 listofgraphs->Add((TObject *)graphDeltaCharge0);
5552 TLegend *legch3 = new TLegend(0.4,0.6,0.89,0.89);
5553 legch3->AddEntry(graphDeltaCharge0,"fit","p");
5554 graphDeltaCharge0->Draw("AP");
5556 TH1I *histoErrorCharge0 = new TH1I("errorcharge0","",100 ,-0.3,0.3);
5557 histoErrorCharge0->SetXTitle("#Deltag/g_{sim}");
5558 histoErrorCharge0->SetYTitle("counts");
5559 histoErrorCharge0->SetLineColor(6);
5560 histoErrorCharge0->SetLineStyle(1);
5561 histoErrorCharge0->SetStats(0);
5562 for(Int_t k = 0; k < counter[0]; k++){
5563 histoErrorCharge0->Fill(yValuesDelta[k]);
5565 TLegend *legch2 = new TLegend(0.4,0.6,0.89,0.89);
5566 legch2->AddEntry(histoErrorCharge0,"f_{g} fit","l");
5567 histoErrorCharge0->Draw();
5568 listofgraphs->Add((TObject *)histoErrorCharge0);
5569 if (fMeanChargeOn) {
5571 Double_t *yValuesDeltaMean = new Double_t[counter[1]];
5572 for(Int_t k = 0; k < counter[1]; k++){
5573 if(fCoefCharge[3][(Int_t)xValuesFittedMean[k]] > 0.0) {
5574 yValuesDeltaMean[k] = (fCoefCharge[1][(Int_t)xValuesFittedMean[k]]-fCoefCharge[3][(Int_t)xValuesFittedMean[k]])/fCoefCharge[3][(Int_t)xValuesFittedMean[k]];
5576 else yValuesDeltaMean[k] = 0.0;
5578 TGraph *graphDeltaCharge1 = new TGraph(counter[1],&xValuesFittedMean[0],yValuesDeltaMean);
5579 graphDeltaCharge1->SetName("deltacharge1");
5580 graphDeltaCharge1->GetXaxis()->SetTitle("Det/Pad groups");
5581 graphDeltaCharge1->GetYaxis()->SetTitle("#Deltag/g_{sim}");
5582 graphDeltaCharge1->SetMarkerColor(2);
5583 graphDeltaCharge1->SetMarkerStyle(24);
5584 graphDeltaCharge1->SetLineColor(2);
5585 graphDeltaCharge1->SetTitle("");
5586 legch3->AddEntry(graphDeltaCharge1,"mean","p");
5587 graphDeltaCharge1->Draw("P");
5588 listofgraphs->Add((TObject *)graphDeltaCharge1);
5590 TH1I *histoErrorCharge1 = new TH1I("errorcharge1","",100 ,-0.3,0.3);
5591 histoErrorCharge1->SetXTitle("#Deltag/g_{sim}");
5592 histoErrorCharge1->SetYTitle("counts");
5593 histoErrorCharge1->SetLineColor(2);
5594 histoErrorCharge1->SetLineStyle(2);
5595 histoErrorCharge1->SetStats(0);
5596 for(Int_t k = 0; k < counter[1]; k++){
5597 histoErrorCharge1->Fill(yValuesDeltaMean[k]);
5599 legch2->AddEntry(histoErrorCharge1,"f_{g} mean","l");
5600 histoErrorCharge1->Draw("same");
5601 listofgraphs->Add((TObject *)histoErrorCharge1);
5604 if (fFitChargeBisOn) {
5606 Double_t *yValuesDeltaBis = new Double_t[counter[2]];
5607 for(Int_t k = 0; k < counter[2]; k++){
5608 if(fCoefCharge[3][(Int_t)xValuesFittedBis[k]] > 0.0) {
5609 yValuesDeltaBis[k] = (fCoefCharge[2][(Int_t)xValuesFittedBis[k]]-fCoefCharge[3][(Int_t)xValuesFittedBis[k]])/fCoefCharge[3][(Int_t)xValuesFittedBis[k]];
5611 else yValuesDeltaBis[k] = 0.0;
5613 TGraph *graphDeltaCharge2 = new TGraph(counter[2],&xValuesFittedBis[0],yValuesDeltaBis);
5614 graphDeltaCharge2->SetName("deltacharge2");
5615 graphDeltaCharge2->GetXaxis()->SetTitle("Det/Pad groups");
5616 graphDeltaCharge2->GetYaxis()->SetTitle("#Deltag/g_{sim}");
5617 graphDeltaCharge2->SetMarkerColor(8);
5618 graphDeltaCharge2->SetLineColor(8);
5619 graphDeltaCharge2->SetMarkerStyle(25);
5620 legch3->AddEntry(graphDeltaCharge2,"fit","p");
5621 graphDeltaCharge2->SetTitle("");
5622 graphDeltaCharge2->Draw("P");
5623 listofgraphs->Add((TObject *)graphDeltaCharge2);
5625 TH1I *histoErrorCharge2 = new TH1I("errorcharge2","",100 ,-0.3,0.3);
5626 histoErrorCharge2->SetXTitle("#Deltag/g_{sim}");
5627 histoErrorCharge2->SetYTitle("counts");
5628 histoErrorCharge2->SetLineColor(8);
5629 histoErrorCharge2->SetLineStyle(5);
5630 histoErrorCharge2->SetLineWidth(3);
5631 histoErrorCharge2->SetStats(0);
5632 for(Int_t k = 0; k < counter[2]; k++){
5633 histoErrorCharge2->Fill(yValuesDeltaBis[k]);
5635 legch2->AddEntry(histoErrorCharge2,"f_{g} fitbis","l");
5636 histoErrorCharge2->Draw("same");
5637 listofgraphs->Add((TObject *)histoErrorCharge2);
5640 legch3->Draw("same");
5642 legch2->Draw("same");
5647 TFile *fout = TFile::Open(fWriteNameCoef,"UPDATE");
5648 // Check if the file could be opened
5649 if (!fout || !fout->IsOpen()) {
5650 AliInfo("No File found!");
5654 for(Int_t k = 0; k < listofgraphs->GetEntriesFast(); k++){
5655 fout->WriteTObject((TObject *) listofgraphs->At(k),((TObject *)listofgraphs->At(k))->GetName(),(Option_t *) "OverWrite");
5663 //_____________________________________________________________________________
5664 void AliTRDCalibra::PlotWritePH()
5667 // create the graph errors and write them if wanted
5670 //TObjArray of the grapherrors and so on
5671 TObjArray *listofgraphs = new TObjArray();
5673 Int_t nbins = fDect2[1]-fDect1[1];
5675 //See the number of fitted for delta
5682 Double_t *xValuesFitted = new Double_t[nbins];
5683 Double_t *xValuesFittedPH = new Double_t[nbins];
5684 for(Int_t k = 0; k < nbins; k ++){
5685 xValuesFitted[k] = -1;
5686 xValuesFittedPH[k] = -1;
5689 for(Int_t l = 0; l < nbins; l++){
5690 if(fCoefVdrift[1][l] > 0){
5691 xValuesFitted[counter[1]]=l;
5696 for(Int_t l = 0; l < nbins; l++){
5697 if(fCoefVdrift[0][l] > 0){
5698 xValuesFittedPH[counter[0]]= l;
5705 //Create the X and Xerror
5706 Double_t *xValues = new Double_t[nbins];
5707 Double_t *xValuesE = new Double_t[nbins];
5708 for(Int_t k = 0; k < nbins; k ++){
5713 //Create the graph erros and plot them
5714 TGraph *graphVdrift2 = new TGraph(nbins,xValues,fCoefVdrift[2]);
5715 graphVdrift2->SetName("coefvdrift2");
5716 graphVdrift2->SetTitle("");
5717 graphVdrift2->GetXaxis()->SetTitle("Det/Pad groups");
5718 graphVdrift2->GetYaxis()->SetTitle("Vdrift [cm/#mus]");
5719 graphVdrift2->SetLineColor(4);
5720 listofgraphs->Add((TObject *)graphVdrift2);
5721 TGraphErrors *graphVdrift1 = new TGraphErrors(nbins,xValues,fCoefVdrift[1],xValuesE,fCoefVdriftE[1]);
5722 graphVdrift1->SetName("coefvdrift1");
5723 graphVdrift1->SetTitle("");
5724 graphVdrift1->GetXaxis()->SetTitle("Det/Pad groups");
5725 graphVdrift1->GetYaxis()->SetTitle("Vdrift [cm/#mus]");
5726 graphVdrift1->SetMarkerColor(6);
5727 graphVdrift1->SetLineColor(6);
5728 graphVdrift1->SetMarkerStyle(26);
5729 listofgraphs->Add((TObject *)graphVdrift1);
5730 TCanvas *cph1 = new TCanvas("cph1","",50,50,600,800);
5732 TLegend *legph1 = new TLegend(0.4,0.6,0.89,0.89);
5733 legph1->AddEntry(graphVdrift2,"Vdrift simulated","l");
5734 legph1->AddEntry(graphVdrift1,"Vdrift fit","p");
5735 graphVdrift2->Draw("AL");
5736 graphVdrift1->Draw("P");
5738 TGraphErrors *graphVdrift0 = new TGraphErrors(nbins,xValues,fCoefVdrift[0],xValuesE,fCoefVdriftE[0]);
5739 graphVdrift0->SetName("coefVdrift0");
5740 graphVdrift0->GetXaxis()->SetTitle("Det/Pad groups");
5741 graphVdrift0->GetYaxis()->SetTitle("Vdrift [cm/#mus]");
5742 graphVdrift0->SetTitle("");
5743 graphVdrift0->SetMarkerColor(2);
5744 graphVdrift0->SetLineColor(2);
5745 graphVdrift0->SetMarkerStyle(24);
5746 legph1->AddEntry(graphVdrift0,"v_{fit PH}","p");
5747 graphVdrift0->Draw("P");
5748 listofgraphs->Add((TObject *)graphVdrift0);
5750 legph1->Draw("same");
5753 //Create the arrays and the graphs for the delta
5754 TCanvas *cph2 = new TCanvas("cph2","",50,50,600,800);
5757 Double_t *yValuesDelta = new Double_t[counter[1]];
5758 for(Int_t k = 0; k < counter[1]; k++){
5759 if(fCoefVdrift[2][(Int_t)(xValuesFitted[k])] > 0.0) {
5760 yValuesDelta[k] = (fCoefVdrift[1][(Int_t)xValuesFitted[k]]-fCoefVdrift[2][(Int_t)xValuesFitted[k]])/fCoefVdrift[2][(Int_t)xValuesFitted[k]];
5762 else yValuesDelta[k] = 0.0;
5764 TGraph *graphDeltaVdrift1 = new TGraph(counter[1],&xValuesFitted[0],yValuesDelta);
5765 graphDeltaVdrift1->SetName("deltavdrift1");
5766 graphDeltaVdrift1->GetXaxis()->SetTitle("Det/Pad groups");
5767 graphDeltaVdrift1->GetYaxis()->SetTitle("#Deltav/v_{sim}");
5768 graphDeltaVdrift1->SetMarkerColor(6);
5769 graphDeltaVdrift1->SetTitle("");
5770 graphDeltaVdrift1->SetLineColor(6);
5771 graphDeltaVdrift1->SetMarkerStyle(26);
5772 listofgraphs->Add((TObject *)graphDeltaVdrift1);
5773 TLegend *legph3 = new TLegend(0.4,0.6,0.89,0.89);
5774 legph3->AddEntry(graphDeltaVdrift1,"v_{slope method}","p");
5775 graphDeltaVdrift1->Draw("AP");
5777 TH1I *histoErrorVdrift1 = new TH1I("errorvdrift1","",100 ,-0.05,0.05);
5778 histoErrorVdrift1->SetXTitle("#Deltav/v_{sim}");
5779 histoErrorVdrift1->SetYTitle("counts");
5780 histoErrorVdrift1->SetLineColor(6);
5781 histoErrorVdrift1->SetLineStyle(1);
5782 histoErrorVdrift1->SetStats(0);
5783 for(Int_t k = 0; k < counter[1]; k++){
5784 histoErrorVdrift1->Fill(yValuesDelta[k]);
5786 TLegend *legph2 = new TLegend(0.4,0.6,0.89,0.89);
5787 legph2->AddEntry(histoErrorVdrift1,"v_{slope method}","l");
5788 histoErrorVdrift1->Draw();
5789 listofgraphs->Add((TObject *)histoErrorVdrift1);
5792 Double_t *yValuesDeltaPH = new Double_t[counter[0]];
5793 for(Int_t k = 0; k < counter[0]; k++){
5794 if(fCoefVdrift[2][(Int_t)xValuesFittedPH[k]] > 0.0) {
5795 yValuesDeltaPH[k] = (fCoefVdrift[0][(Int_t)xValuesFittedPH[k]]-fCoefVdrift[2][(Int_t)xValuesFittedPH[k]])/fCoefVdrift[2][(Int_t)xValuesFittedPH[k]];
5797 else yValuesDeltaPH[k] = 0.0;
5799 TGraph *graphDeltaVdrift0 = new TGraph(counter[0],&xValuesFittedPH[0],yValuesDeltaPH);
5800 graphDeltaVdrift0->SetName("deltavdrift0");
5801 graphDeltaVdrift0->GetXaxis()->SetTitle("Det/Pad groups");
5802 graphDeltaVdrift0->GetYaxis()->SetTitle("#Deltav/v_{sim}");
5803 graphDeltaVdrift0->SetMarkerColor(2);
5804 graphDeltaVdrift0->SetMarkerStyle(24);
5805 graphDeltaVdrift0->SetLineColor(2);
5806 graphDeltaVdrift0->SetTitle("");
5807 legph3->AddEntry(graphDeltaVdrift0,"v_{fit PH}","p");
5808 graphDeltaVdrift0->Draw("P");
5809 listofgraphs->Add((TObject *)graphDeltaVdrift0);
5811 TH1I *histoErrorVdrift0 = new TH1I("errorvdrift0","",100 ,-0.05,0.05);
5812 histoErrorVdrift0->SetXTitle("#Deltav/v_{sim}");
5813 histoErrorVdrift0->SetYTitle("counts");
5814 histoErrorVdrift0->SetLineColor(2);
5815 histoErrorVdrift0->SetLineStyle(2);
5816 histoErrorVdrift0->SetStats(0);
5817 for(Int_t k = 0; k < counter[0]; k++){
5818 histoErrorVdrift0->Fill(yValuesDeltaPH[k]);
5820 legph2->AddEntry(histoErrorVdrift0,"v_{fit PH}","l");
5821 histoErrorVdrift0->Draw("same");
5822 listofgraphs->Add((TObject *)histoErrorVdrift0);
5825 legph3->Draw("same");
5827 legph2->Draw("same");
5832 TFile *fout = TFile::Open(fWriteNameCoef,"UPDATE");
5833 // Check if the file could be opened
5834 if (!fout || !fout->IsOpen()) {
5835 AliInfo("No File found!");
5839 for(Int_t k = 0; k < listofgraphs->GetEntriesFast(); k++){
5840 fout->WriteTObject((TObject *) listofgraphs->At(k),((TObject *)listofgraphs->At(k))->GetName(),(Option_t *) "OverWrite");
5848 //_____________________________________________________________________________
5849 void AliTRDCalibra::PlotWriteT0()
5852 // create the graph errors and write them if wanted
5855 //TObjArray of the grapherrors and so on
5856 TObjArray *listofgraphs = new TObjArray();
5858 Int_t nbins = fDect2[1]-fDect1[1];
5860 //See the number of fitted for delta: here T0 can be negative, we don't use the sign but the error
5861 //and the grapherrors of the coefficients contained the no fitted with error 0.0
5868 Double_t *xValuesFitted = new Double_t[nbins];
5869 Double_t *xValuesFittedPH = new Double_t[nbins];
5870 for(Int_t k = 0; k < nbins; k ++){
5871 xValuesFitted[k] = -1;
5872 xValuesFittedPH[k] = -1;
5875 for(Int_t l = 0; l < nbins; l++){
5876 if(fCoefT0E[1][l] != 0.0){
5877 xValuesFitted[counter[1]]=l;
5882 for(Int_t l = 0; l < nbins; l++){
5883 if(fCoefT0E[0][l] != 0.0){
5884 xValuesFittedPH[counter[0]]= l;
5891 //Create the X and Xerror
5892 Double_t *xValues = new Double_t[nbins];
5893 Double_t *xValuesE = new Double_t[nbins];
5894 for(Int_t k = 0; k < nbins; k ++){
5899 //Create the graph erros and plot them
5900 TGraph *graphT02 = new TGraph(nbins,xValues,fCoefT0[2]);
5901 graphT02->SetName("coeft02");
5902 graphT02->SetTitle("");
5903 graphT02->GetXaxis()->SetTitle("Det/Pad groups");
5904 graphT02->GetYaxis()->SetTitle("T0 [time bins]");
5905 graphT02->SetLineColor(4);
5906 listofgraphs->Add((TObject *)graphT02);
5907 TGraphErrors *graphT01 = new TGraphErrors(nbins,xValues,fCoefT0[1],xValuesE,fCoefT0E[1]);
5908 graphT01->SetName("coeft01");
5909 graphT01->SetTitle("");
5910 graphT01->GetXaxis()->SetTitle("Det/Pad groups");
5911 graphT01->GetYaxis()->SetTitle("T0 [time bins]");
5912 graphT01->SetMarkerColor(6);
5913 graphT01->SetLineColor(6);
5914 graphT01->SetMarkerStyle(26);
5915 listofgraphs->Add((TObject *)graphT01);
5916 TCanvas *ct01 = new TCanvas("ct01","",50,50,600,800);
5918 TLegend *legt01 = new TLegend(0.4,0.6,0.89,0.89);
5919 legt01->AddEntry(graphT02,"T0 simulated","l");
5920 legt01->AddEntry(graphT01,"T0 slope method","p");
5921 graphT02->Draw("AL");
5922 graphT01->Draw("P");
5924 TGraphErrors *graphT00 = new TGraphErrors(nbins,xValues,fCoefT0[0],xValuesE,fCoefT0E[0]);
5925 graphT00->SetName("coeft00");
5926 graphT00->GetXaxis()->SetTitle("Det/Pad groups");
5927 graphT00->GetYaxis()->SetTitle("T0 [time bins]");
5928 graphT00->SetTitle("");
5929 graphT00->SetMarkerColor(2);
5930 graphT00->SetLineColor(2);
5931 graphT00->SetMarkerStyle(24);
5932 legt01->AddEntry(graphT00,"T0 fit","p");
5933 graphT00->Draw("P");
5934 listofgraphs->Add((TObject *)graphT00);
5936 legt01->Draw("same");
5939 //Create the arrays and the graphs for the delta
5940 TCanvas *ct02 = new TCanvas("ct02","",50,50,600,800);
5943 Double_t *yValuesDelta = new Double_t[counter[1]];
5944 for(Int_t k = 0; k < counter[1]; k++){
5945 yValuesDelta[k] = (fCoefT0[1][(Int_t)xValuesFitted[k]]-fCoefT0[2][(Int_t)xValuesFitted[k]]);
5947 TGraph *graphDeltaT01 = new TGraph(counter[1],&xValuesFitted[0],yValuesDelta);
5948 graphDeltaT01->SetName("deltat01");
5949 graphDeltaT01->GetXaxis()->SetTitle("Det/Pad groups");
5950 graphDeltaT01->GetYaxis()->SetTitle("#Deltat0 [time bins]");
5951 graphDeltaT01->SetMarkerColor(6);
5952 graphDeltaT01->SetTitle("");
5953 graphDeltaT01->SetLineColor(6);
5954 graphDeltaT01->SetMarkerStyle(26);
5955 listofgraphs->Add((TObject *)graphDeltaT01);
5956 TLegend *legt03 = new TLegend(0.4,0.6,0.89,0.89);
5957 legt03->AddEntry(graphDeltaT01,"T0_{slope method}","p");
5958 graphDeltaT01->Draw("AP");
5960 TH1I *histoErrorT01 = new TH1I("errort01","",100 ,-0.05,0.05);
5961 histoErrorT01->SetXTitle("#Deltat0 [time bins]");
5962 histoErrorT01->SetYTitle("counts");
5963 histoErrorT01->SetLineColor(6);
5964 histoErrorT01->SetLineStyle(1);
5965 histoErrorT01->SetStats(0);
5966 for(Int_t k = 0; k < counter[1]; k++){
5967 histoErrorT01->Fill(yValuesDelta[k]);
5969 TLegend *legt02 = new TLegend(0.4,0.6,0.89,0.89);
5970 legt02->AddEntry(histoErrorT01,"T0_{slope method}","l");
5971 histoErrorT01->Draw();
5972 listofgraphs->Add((TObject *)histoErrorT01);
5975 Double_t *yValuesDeltaPH = new Double_t[counter[0]];
5976 for(Int_t k = 0; k < counter[0]; k++){
5977 yValuesDeltaPH[k] = (fCoefT0[0][(Int_t)xValuesFittedPH[k]]-fCoefT0[2][(Int_t)xValuesFittedPH[k]]);
5979 TGraph *graphDeltaT00 = new TGraph(counter[0],&xValuesFittedPH[0],yValuesDeltaPH);
5980 graphDeltaT00->SetName("deltat00");
5981 graphDeltaT00->GetXaxis()->SetTitle("Det/Pad groups");
5982 graphDeltaT00->GetYaxis()->SetTitle("#Deltat0 [time bins]");
5983 graphDeltaT00->SetMarkerColor(2);
5984 graphDeltaT00->SetMarkerStyle(24);
5985 graphDeltaT00->SetLineColor(2);
5986 graphDeltaT00->SetTitle("");
5987 legt03->AddEntry(graphDeltaT00,"T0_{fit PH}","p");
5988 graphDeltaT00->Draw("P");
5989 listofgraphs->Add((TObject *)graphDeltaT00);
5991 TH1I *histoErrorT00 = new TH1I("errort00","",100 ,-0.05,0.05);
5992 histoErrorT00->SetXTitle("#Deltat0 [time bins]");
5993 histoErrorT00->SetYTitle("counts");
5994 histoErrorT00->SetLineColor(2);
5995 histoErrorT00->SetLineStyle(2);
5996 histoErrorT00->SetStats(0);
5997 for(Int_t k = 0; k < counter[0]; k++){
5998 histoErrorT00->Fill(yValuesDeltaPH[k]);
6000 legt02->AddEntry(histoErrorT00,"T0_{fit PH}","l");
6001 histoErrorT00->Draw("same");
6002 listofgraphs->Add((TObject *)histoErrorT00);
6005 legt03->Draw("same");
6007 legt02->Draw("same");
6012 TFile *fout = TFile::Open(fWriteNameCoef,"UPDATE");
6013 // Check if the file could be opened
6014 if (!fout || !fout->IsOpen()) {
6015 AliInfo("No File found!");
6019 for(Int_t k = 0; k < listofgraphs->GetEntriesFast(); k++){
6020 fout->WriteTObject((TObject *) listofgraphs->At(k),((TObject *)listofgraphs->At(k))->GetName(),(Option_t *) "OverWrite");
6028 //_____________________________________________________________________________
6029 void AliTRDCalibra::PlotWritePRF()
6032 // create the graph errors and write them if wanted
6035 //TObjArray of the grapherrors and so on
6036 TObjArray *listofgraphs = new TObjArray();
6038 Int_t nbins = fDect2[2]-fDect1[2];
6040 //See the number of fitted for delta
6045 Double_t *xValuesFitted = new Double_t[nbins];
6046 for(Int_t k = 0; k < nbins; k ++){
6047 xValuesFitted[k] = -1;
6050 for(Int_t l = 0; l < nbins; l++){
6051 if(fCoefPRF[0][l] > 0){
6052 xValuesFitted[counter]=l;
6059 //Create the X and Xerror
6060 Double_t *xValues = new Double_t[nbins];
6061 Double_t *xValuesE = new Double_t[nbins];
6062 for(Int_t k = 0; k < nbins; k ++){
6067 //Create the graph erros and plot them
6068 TGraph *graphPRF1 = new TGraph(nbins,xValues,fCoefPRF[1]);
6069 graphPRF1->SetName("coefprf1");
6070 graphPRF1->SetTitle("");
6071 graphPRF1->GetXaxis()->SetTitle("Det/Pad groups");
6072 graphPRF1->GetYaxis()->SetTitle("PRF width [p.u]");
6073 graphPRF1->SetLineColor(4);
6074 graphPRF1->SetMarkerColor(4);
6075 graphPRF1->SetMarkerStyle(25);
6076 graphPRF1->SetMarkerSize(0.7);
6077 listofgraphs->Add((TObject *)graphPRF1);
6078 TGraphErrors *graphPRF0 = new TGraphErrors(nbins,xValues,fCoefPRF[0],xValuesE,fCoefPRFE);
6079 graphPRF0->SetName("coefprf0");
6080 graphPRF0->SetTitle("");
6081 graphPRF0->GetXaxis()->SetTitle("Det/Pad groups");
6082 graphPRF0->GetYaxis()->SetTitle("PRF Width [p.u]");
6083 graphPRF0->SetMarkerColor(6);
6084 graphPRF0->SetLineColor(6);
6085 graphPRF0->SetMarkerStyle(26);
6086 listofgraphs->Add((TObject *)graphPRF0);
6087 TCanvas *cprf1 = new TCanvas("cprf1","",50,50,600,800);
6089 TLegend *legprf1 = new TLegend(0.4,0.6,0.89,0.89);
6090 legprf1->AddEntry(graphPRF1,"PRF width simulated","p");
6091 legprf1->AddEntry(graphPRF0,"PRF fit","p");
6092 graphPRF1->Draw("AP");
6093 graphPRF0->Draw("P");
6094 legprf1->Draw("same");
6097 //Create the arrays and the graphs for the delta
6098 TCanvas *cprf2 = new TCanvas("cprf2","",50,50,600,800);
6101 Double_t *yValuesDelta = new Double_t[counter];
6102 for(Int_t k = 0; k < counter; k++){
6103 if(fCoefPRF[1][(Int_t)xValuesFitted[k]] > 0.0){
6104 yValuesDelta[k] = (fCoefPRF[0][(Int_t)xValuesFitted[k]]-fCoefPRF[1][(Int_t)xValuesFitted[k]])/(fCoefPRF[1][(Int_t)xValuesFitted[k]]);
6107 TGraph *graphDeltaPRF = new TGraph(counter,&xValuesFitted[0],yValuesDelta);
6108 graphDeltaPRF->SetName("deltaprf");
6109 graphDeltaPRF->GetXaxis()->SetTitle("Det/Pad groups");
6110 graphDeltaPRF->GetYaxis()->SetTitle("#Delta#sigma/#sigma_{sim}");
6111 graphDeltaPRF->SetMarkerColor(6);
6112 graphDeltaPRF->SetTitle("");
6113 graphDeltaPRF->SetLineColor(6);
6114 graphDeltaPRF->SetMarkerStyle(26);
6115 listofgraphs->Add((TObject *)graphDeltaPRF);
6116 TLegend *legprf3 = new TLegend(0.4,0.6,0.89,0.89);
6117 legprf3->AddEntry(graphDeltaPRF,"#sigma_{fit}","p");
6118 graphDeltaPRF->Draw("AP");
6120 TH1I *histoErrorPRF = new TH1I("errorprf1","",100 ,-0.5,0.5);
6121 histoErrorPRF->SetXTitle("#Delta#sigma/#sigma_{sim}");
6122 histoErrorPRF->SetYTitle("counts");
6123 histoErrorPRF->SetLineColor(6);
6124 histoErrorPRF->SetLineStyle(1);
6125 histoErrorPRF->SetStats(0);
6126 for(Int_t k = 0; k < counter; k++){
6127 histoErrorPRF->Fill(yValuesDelta[k]);
6129 TLegend *legprf2 = new TLegend(0.4,0.6,0.89,0.89);
6130 legprf2->AddEntry(histoErrorPRF,"#sigma_{fit}","l");
6131 histoErrorPRF->Draw();
6132 listofgraphs->Add((TObject *)histoErrorPRF);
6134 legprf3->Draw("same");
6136 legprf2->Draw("same");
6141 TFile *fout = TFile::Open(fWriteNameCoef,"UPDATE");
6142 // Check if the file could be opened
6143 if (!fout || !fout->IsOpen()) {
6144 AliInfo("No File found!");
6148 for(Int_t k = 0; k < listofgraphs->GetEntriesFast(); k++){
6149 fout->WriteTObject((TObject *) listofgraphs->At(k),((TObject *)listofgraphs->At(k))->GetName(),(Option_t *) "OverWrite");
6158 //____________Plot histos DB___________________________________________________
6161 //_____________________________________________________________________________
6162 void AliTRDCalibra::PlotCHDB()
6165 // Plot the histos for fDebug = 3 and fDebug = 4 to visualise the detector
6168 TCanvas *cchdb = new TCanvas("cchdb","",50,50,600,800);
6169 if ((fFitChargeBisOn) && (fMeanChargeOn)) {
6172 fCoefChargeDB[0]->Draw("LEGO");
6174 fCoefChargeDB[1]->Draw("LEGO");
6176 fCoefChargeDB[2]->Draw("LEGO");
6178 if ((!fFitChargeBisOn) && (fMeanChargeOn)) {
6181 fCoefChargeDB[0]->Draw("LEGO");
6183 fCoefChargeDB[1]->Draw("LEGO");
6187 fCoefChargeDB[0]->Draw("LEGO");
6192 //_____________________________________________________________________________
6193 void AliTRDCalibra::PlotPHDB()
6196 // Plot the histos for fDebug = 3 and fDebug = 4 to visualise the detector
6199 TCanvas *cphdb = new TCanvas("cphdb","",50,50,600,800);
6203 fCoefVdriftDB[0]->Draw("LEGO");
6205 fCoefVdriftDB[1]->Draw("LEGO");
6209 fCoefVdriftDB[1]->Draw("LEGO");
6214 //_____________________________________________________________________________
6215 void AliTRDCalibra::PlotT0DB()
6218 // Plot the histos for fDebug = 3 and fDebug = 4 to visualise the detector
6221 TCanvas *ct0db = new TCanvas("ct0db","",50,50,600,800);
6225 fCoefT0DB[0]->Draw("LEGO");
6227 fCoefT0DB[1]->Draw("LEGO");
6231 fCoefT0DB[1]->Draw("LEGO");
6236 //_____________________________________________________________________________
6237 void AliTRDCalibra::PlotPRFDB()
6240 // Plot the histos for fDebug = 3 and fDebug = 4 to visualise the detector
6243 TCanvas *cprfdb = new TCanvas("cprfdb","",50,50,600,800);
6245 fCoefPRFDB->Draw("LEGO");
6250 //____________Write DB Histos__________________________________________________
6253 //_____________________________________________________________________________
6254 void AliTRDCalibra::WriteCHDB(TFile *fout)
6257 // If wanted, write the debug histos for fDebug = 3 and fDebug = 4
6260 fout->WriteTObject(fCoefChargeDB[0],fCoefChargeDB[0]->GetName(),(Option_t *) "OverWrite");
6261 if (fMeanChargeOn) {
6262 fout->WriteTObject(fCoefChargeDB[1],fCoefChargeDB[1]->GetName(),(Option_t *) "OverWrite");
6264 if (fFitChargeBisOn ) {
6265 fout->WriteTObject(fCoefChargeDB[2],fCoefChargeDB[2]->GetName(),(Option_t *) "OverWrite");
6270 //_____________________________________________________________________________
6271 void AliTRDCalibra::WritePHDB(TFile *fout)
6274 // If wanted, write the debug histos for fDebug = 3 and fDebug = 4
6278 fout->WriteTObject(fCoefVdriftDB[0],fCoefVdriftDB[0]->GetName(),(Option_t *) "OverWrite");
6280 fout->WriteTObject(fCoefVdriftDB[1],fCoefVdriftDB[1]->GetName(),(Option_t *) "OverWrite");
6284 //_____________________________________________________________________________
6285 void AliTRDCalibra::WriteT0DB(TFile *fout)
6288 // If wanted, write the debug histos for fDebug = 3 and fDebug = 4
6292 fout->WriteTObject(fCoefT0DB[0],fCoefT0DB[0]->GetName(),(Option_t *) "OverWrite");
6294 fout->WriteTObject(fCoefT0DB[1],fCoefT0DB[1]->GetName(),(Option_t *) "OverWrite");
6298 //_____________________________________________________________________________
6299 void AliTRDCalibra::WritePRFDB(TFile *fout)
6302 // If wanted, write the debug histos for fDebug = 3 and fDebug = 4
6305 fout->WriteTObject(fCoefPRFDB,fCoefPRFDB->GetName(),(Option_t *) "OverWrite");
6310 //____________Calcul Coef Mean_________________________________________________
6313 //_____________________________________________________________________________
6314 Bool_t AliTRDCalibra::CalculT0CoefMean(Int_t dect, Int_t idect)
6317 // For the detector Dect calcul the mean time 0
6318 // for the calibration group idect from the choosen database
6321 AliTRDcalibDB *cal = AliTRDcalibDB::Instance();
6323 AliInfo("Could not get calibDB Manager");
6329 if ((fDebug != 2) && fAccCDB) {
6331 for (Int_t row = fRowMin[1]; row < fRowMax[1]; row++) {
6332 for (Int_t col = fColMin[1]; col < fColMax[1]; col++) {
6336 fT0Coef[2] += (Float_t) cal->GetT0(dect,col,row);
6340 fT0Coef[2] += (Float_t) cal->GetT0Average(dect);
6345 fT0Coef[2] = fT0Coef[2] / ((fColMax[1]-fColMin[1])*(fRowMax[1]-fRowMin[1]));
6346 if ((fDebug == 1) ||
6348 fCoefT0[2][idect] = fT0Coef[2];
6357 //_____________________________________________________________________________
6358 Bool_t AliTRDCalibra::CalculChargeCoefMean(Int_t dect, Int_t idect, Bool_t vrai)
6361 // For the detector Dect calcul the mean gain factor
6362 // for the calibration group idect from the choosen database
6365 AliTRDcalibDB *cal = AliTRDcalibDB::Instance();
6367 AliInfo("Could not get calibDB Manager");
6370 AliTRDCommonParam *parCom = AliTRDCommonParam::Instance();
6372 AliInfo("Could not get CommonParam Manager");
6376 fChargeCoef[3] = 0.0;
6380 for (Int_t row = fRowMin[0]; row < fRowMax[0]; row++) {
6381 for (Int_t col = fColMin[0]; col < fColMax[0]; col++) {
6386 fChargeCoef[3] += (Float_t) cal->GetGainFactor(dect,col,row);
6388 if (vrai && fAccCDB) {
6389 fScaleFitFactor += (Float_t) cal->GetGainFactor(dect,col,row);
6392 fChargeCoef[3] += 1.0;
6394 if (vrai && (!fAccCDB)) {
6395 fScaleFitFactor += 1.0;
6401 fChargeCoef[3] += (Float_t) cal->GetGainFactorAverage(dect);
6403 if (vrai && fAccCDB) {
6404 fScaleFitFactor += ((Float_t) cal->GetGainFactorAverage(dect));
6407 fChargeCoef[3] += 1.0;
6409 if (vrai && (!fAccCDB)) {
6410 fScaleFitFactor += 1.0;
6416 fChargeCoef[3] = fChargeCoef[3] / ((fColMax[0]-fColMin[0])*(fRowMax[0]-fRowMin[0]));
6417 if ((fDebug == 1) ||
6419 fCoefCharge[3][idect]=fChargeCoef[3];
6428 //_____________________________________________________________________________
6429 Bool_t AliTRDCalibra::CalculPRFCoefMean(Int_t dect, Int_t idect)
6432 // For the detector Dect calcul the mean sigma of pad response
6433 // function for the calibration group idect from the choosen database
6436 AliTRDcalibDB *cal = AliTRDcalibDB::Instance();
6438 AliInfo("Could not get calibDB Manager");
6442 AliTRDCommonParam *parCom = AliTRDCommonParam::Instance();
6444 AliInfo("Could not get CommonParam Manager");
6453 for (Int_t row = fRowMin[2]; row < fRowMax[2]; row++) {
6454 for (Int_t col = fColMin[2]; col < fColMax[2]; col++) {
6455 if ((parCom->GetColMax(GetPlane(dect)) != (col+1)) && (col != 0)) {
6458 fPRFCoef[1] += (Float_t) cal->GetPRFWidth(dect,col,row);
6461 fPRFCoef[1] += GetPRFDefault(GetPlane(dect));
6468 fPRFCoef[1] = fPRFCoef[1]/cot;
6469 if ((fDebug == 1) ||
6471 fCoefPRF[1][idect] = fPRFCoef[1];
6475 if ((fDebug == 1) ||
6478 fCoefPRF[1][idect] = cal->GetPRFWidth(dect,fColMin[2],fRowMin[2]);
6481 fCoefPRF[1][idect] = GetPRFDefault(GetPlane(dect));
6492 //_____________________________________________________________________________
6493 Bool_t AliTRDCalibra::CalculVdriftCoefMean(Int_t dect, Int_t idect)
6496 // For the detector dect calcul the mean drift velocity for the
6497 // calibration group idect from the choosen database
6500 AliTRDcalibDB *cal = AliTRDcalibDB::Instance();
6502 AliInfo("Could not get calibDB Manager");
6506 fVdriftCoef[2] = 0.0;
6509 for (Int_t row = fRowMin[1]; row < fRowMax[1]; row++) {
6510 for (Int_t col = fColMin[1]; col < fColMax[1]; col++) {
6515 fVdriftCoef[2] += (Float_t) cal->GetVdrift(dect,col,row);
6518 fVdriftCoef[2] += 1.5;
6524 fVdriftCoef[2] += (Float_t) cal->GetVdriftAverage(dect);
6527 fVdriftCoef[2] += 1.5;
6532 fVdriftCoef[2] = fVdriftCoef[2] / ((fColMax[1]-fColMin[1])*(fRowMax[1]-fRowMin[1]));
6533 if ((fDebug == 1) ||
6535 fCoefVdrift[2][idect] = fVdriftCoef[2];
6543 //_____________________________________________________________________________
6544 Float_t AliTRDCalibra::GetPRFDefault(Int_t plane) const
6547 // Default width of the PRF if there is no database as reference
6575 //____________Pad group calibration mode_______________________________________
6578 //_____________________________________________________________________________
6579 void AliTRDCalibra::ReconstructionRowPadGroup(Int_t idect, Int_t i)
6582 // For the calibration group idect in a detector calculate the
6583 // first and last row pad and col pad.
6584 // The pads in the interval will have the same calibrated coefficients
6594 if (fNfragZ[i] != 0) {
6595 posc = (Int_t) idect / fNfragZ[i];
6597 if (fNfragRphi[i] != 0) {
6598 posr = (Int_t) idect % fNfragZ[i];
6600 fRowMin[i] = posr * fNnZ[i];
6601 fRowMax[i] = (posr+1) * fNnZ[i];
6602 fColMin[i] = posc * fNnRphi[i];
6603 fColMax[i] = (posc+1) * fNnRphi[i];
6607 //_____________________________________________________________________________
6608 void AliTRDCalibra::CalculXBins(Int_t idect, Int_t i)
6611 // For the detector idect calcul the first Xbins
6616 AliInfo(Form("detector: %d", idect));
6620 Int_t sector = GetSector(idect);
6621 fXbins[i] += sector*(6*fDetChamb2[i]+6*4*fDetChamb0[i]);
6623 // In which chamber?
6624 Int_t chamber = GetChamber(idect);
6626 while (kc < chamber) {
6628 fXbins[i] += 6 * fDetChamb2[i];
6631 fXbins[i] += 6 * fDetChamb0[i];
6637 Int_t plane = GetPlane(idect);
6639 fXbins[i] += plane*fDetChamb2[i];
6642 fXbins[i] += plane*fDetChamb0[i];
6647 //_____________________________________________________________________________
6648 Int_t AliTRDCalibra::SearchInVector(Int_t group, Int_t i) const
6651 // Search if the calibration group "group" has already been
6652 // initialised by a previous track in the vector
6656 for (Int_t k = 0; k < (Int_t) fPlaCH->GetEntriesFast(); k++) {
6657 if (((AliTRDPlace *) fPlaCH->At(k))->GetPlace() == group) {
6665 for (Int_t k = 0; k < (Int_t) fPlaPH->GetEntriesFast(); k++) {
6666 if (((AliTRDPlace *) fPlaPH->At(k))->GetPlace() == group) {
6674 for (Int_t k = 0; k < (Int_t) fPlaPRF->GetEntriesFast(); k++) {
6675 if (((AliTRDPlace *) fPlaPRF->At(k))->GetPlace() == group) {
6686 //_____________________________________________________________________________
6687 Int_t AliTRDCalibra::SearchInTreeVector(TObjArray *vectorplace, Int_t group) const
6690 // Search if the calibration group "group" is present in the tree
6693 for (Int_t k = 0; k < (Int_t) vectorplace->GetEntriesFast(); k++) {
6694 if (((AliTRDPlace *) vectorplace->At(k))->GetPlace() == group) {
6703 //_____________________________________________________________________________
6704 Int_t AliTRDCalibra::SearchBin(Float_t value, Int_t i) const
6712 Int_t fbinmax = (Int_t) value;
6713 Int_t fNumberOfBin = -1;
6719 fNumberOfBin = fNumberBinCharge;
6726 fNumberOfBin = fNumberBinPRF;
6730 if ((value >= fbinmax) ||
6731 (value < fbinmin)) {
6736 reponse = (Int_t) ((fNumberOfBin*(value-fbinmin)) / (fbinmax-fbinmin));
6743 //_____________________________________________________________________________
6744 Bool_t AliTRDCalibra::UpdateVectorCH(Int_t group, Float_t value)
6747 // Fill the vector if a new calibration group "group" or update the
6748 // values of the calibration group "group" if already here
6752 Int_t bin = SearchBin(value,0);
6754 if ((bin < 0) || (bin >= fNumberBinCharge)) {
6759 Int_t place = SearchInVector(group,0);
6763 AliTRDPlace *placegroup = new AliTRDPlace();
6764 placegroup->SetPlace(group);
6765 fPlaCH->Add((TObject *) placegroup);
6767 AliTRDCTInfo *fCHInfo = new AliTRDCTInfo();
6768 UShort_t *entries = new UShort_t[fNumberBinCharge];
6770 for(Int_t k = 0; k < fNumberBinCharge; k++) {
6776 fCHInfo->SetEntries(entries);
6777 // Set in the vector
6778 fVectorCH->Add((TObject *) fCHInfo);
6780 // Group already exits
6783 AliTRDCTInfo *fCHInfo = new AliTRDCTInfo();
6785 fCHInfo = ((AliTRDCTInfo *) fVectorCH->At(place));
6786 UShort_t *entries = fCHInfo->GetEntries();
6790 fCHInfo->SetEntries(entries);
6791 // Update the vector
6792 fVectorCH->AddAt((TObject *) fCHInfo,place);
6799 //_____________________________________________________________________________
6800 Bool_t AliTRDCalibra::UpdateVectorPRF(Int_t group, Float_t x, Float_t y)
6803 // Fill the vector if a new calibration group "group" or update the
6804 // values of the calibration group "group" if already here
6808 Int_t bin = SearchBin(x,2);
6810 if ((bin < 0) || (bin >= fNumberBinPRF)) {
6815 Int_t place = SearchInVector(group,2);
6820 AliTRDPlace *placegroup = new AliTRDPlace();
6821 placegroup->SetPlace(group);
6822 fPlaPRF->Add((TObject *) placegroup);
6823 AliTRDPInfo *fPRFInfo = new AliTRDPInfo();
6825 Float_t *sum = new Float_t[fNumberBinPRF];
6826 Float_t *sumsquare = new Float_t[fNumberBinPRF];
6827 UShort_t *entries = new UShort_t[fNumberBinPRF];
6830 for (Int_t k = 0; k < fNumberBinPRF; k++) {
6838 sumsquare[bin] += y*y;
6842 fPRFInfo->SetSum(sum);
6843 fPRFInfo->SetSumSquare(sumsquare);
6844 fPRFInfo->SetEntries(entries);
6846 // Set in the vector
6847 fVectorPRF->Add((TObject *) fPRFInfo);
6850 // Group already exits
6853 AliTRDPInfo *fPRFInfo = new AliTRDPInfo();
6855 fPRFInfo = (AliTRDPInfo *) fVectorPRF->At(place);
6857 Float_t *sum = fPRFInfo->GetSum();
6858 Float_t *sumsquare = fPRFInfo->GetSumSquare();
6859 UShort_t *entries = fPRFInfo->GetEntries();
6862 Double_t calcul = (((Double_t) fPRFInfo->GetEntries()[bin])
6863 * ((Double_t) fPRFInfo->GetSum()[bin]) + (Double_t) y)
6864 / (((Double_t) fPRFInfo->GetEntries()[bin]) + 1);
6865 sum[bin] = (Float_t) calcul;
6866 Double_t calculsquare = (((Double_t) fPRFInfo->GetSumSquare()[bin])
6867 * ((Double_t) fPRFInfo->GetEntries()[bin]) + ((Double_t) y)*((Double_t) y))
6868 / (((Double_t) fPRFInfo->GetEntries()[bin]) + 1);
6869 sumsquare[bin] = (Float_t) calculsquare;
6873 fPRFInfo->SetSum(sum);
6874 fPRFInfo->SetSumSquare(sumsquare);
6875 fPRFInfo->SetEntries(entries);
6877 // Update the vector
6878 fVectorPRF->AddAt((TObject *) fPRFInfo,place);
6886 //_____________________________________________________________________________
6887 Bool_t AliTRDCalibra::UpdateVectorPH(Int_t group, Int_t time, Float_t value)
6890 // Fill the vector if a new calibration group "group" or update
6891 // the values of the calibration group "group" if already here
6898 (bin >= fTimeMax)) {
6903 Int_t place = SearchInVector(group,1);
6908 AliTRDPlace *placegroup = new AliTRDPlace();
6909 placegroup->SetPlace(group);
6910 fPlaPH->Add((TObject *) placegroup);
6911 AliTRDPInfo *fPHInfo = new AliTRDPInfo();
6913 Float_t *sum = new Float_t[fTimeMax];
6914 Float_t *sumsquare = new Float_t[fTimeMax];
6915 UShort_t *entries = new UShort_t[fTimeMax];
6918 for (Int_t k = 0; k < fTimeMax; k++) {
6926 sumsquare[bin] += value*value;
6930 fPHInfo->SetSum(sum);
6931 fPHInfo->SetSumSquare(sumsquare);
6932 fPHInfo->SetEntries(entries);
6934 // Set in the vector
6935 fVectorPH->Add((TObject *) fPHInfo);
6938 // Group already exits
6941 AliTRDPInfo *fPHInfo = new AliTRDPInfo();
6943 fPHInfo = (AliTRDPInfo *) fVectorPH->At(place);
6945 Float_t *sum = fPHInfo->GetSum();
6946 Float_t *sumsquare = fPHInfo->GetSumSquare();
6947 UShort_t *entries = fPHInfo->GetEntries();
6950 Double_t calcul = (((Double_t) fPHInfo->GetEntries()[bin])
6951 * ((Double_t) fPHInfo->GetSum()[bin]) + (Double_t) value)
6952 / (((Double_t) fPHInfo->GetEntries()[bin]) + 1);
6953 sum[bin] = (Float_t) calcul;
6954 Double_t calculsquare = ((((Double_t) fPHInfo->GetSumSquare()[bin])
6955 * ((Double_t) fPHInfo->GetEntries()[bin]))
6956 + (((Double_t) value) * ((Double_t)value)))
6957 / (((Double_t) fPHInfo->GetEntries()[bin]) + 1);
6958 sumsquare[bin] = (Float_t) calculsquare;
6962 fPHInfo->SetSum(sum);
6963 fPHInfo->SetSumSquare(sumsquare);
6964 fPHInfo->SetEntries(entries);
6966 // Update the vector
6967 fVectorPH->AddAt((TObject *) fPHInfo,place);
6975 //_____________________________________________________________________________
6976 TGraphErrors *AliTRDCalibra::ConvertVectorPHisto(AliTRDPInfo *pInfo
6977 , const Char_t *name) const
6980 // Convert the PInfo in a 1D grapherror, name must contains "PRF"
6981 // if PRF calibration and not "PRF" for Vdrift calibration
6984 TGraphErrors *histo;
6985 const Char_t *pattern1 = "PRF";
6992 Double_t step = 0.0;
6997 if (strstr(name,pattern1)) {
6998 ntimes = fNumberBinPRF;
7003 x = new Double_t[ntimes]; // Xaxis
7004 y = new Double_t[ntimes]; // Mean
7005 ex = new Double_t[ntimes]; // Nentries
7006 ey = new Double_t[ntimes]; // Sum of square/nentries
7009 if (!strstr(name,pattern1)) {
7014 step = (1.0 - (-1.0)) / fNumberBinPRF;
7015 min = -1.0 + step / 2.0;
7019 for (Int_t k = 0; k < ntimes; k++) {
7020 x[k] = min + k*step;
7024 // Fill only if there is more than 0 something
7025 if (pInfo->GetEntries()[k] > 0) {
7026 ex[k] = pInfo->GetEntries()[k];
7027 y[k] = pInfo->GetSum()[k];
7028 ey[k] = pInfo->GetSumSquare()[k];
7032 // Define the TGraphErrors
7033 histo = new TGraphErrors(ntimes,x,y,ex,ey);
7034 histo->SetTitle(name);
7039 //_____________________________________________________________________________
7040 TH1F *AliTRDCalibra::ConvertVectorCTHisto(AliTRDCTInfo *cTInfo
7041 , const Char_t * name) const
7044 // Convert the CTInfo in a 1D histo
7049 Int_t ntimes = fNumberBinCharge;
7050 UShort_t *entries = cTInfo->GetEntries();
7053 histo = new TH1F(name,name,fNumberBinCharge,0,300);
7056 for (Int_t k = 0; k < ntimes; k++) {
7057 histo->SetBinContent(k+1,entries[k]);
7058 histo->SetBinError(k+1,TMath::Sqrt(TMath::Abs(entries[k])));
7065 //_____________________________________________________________________________
7066 TTree *AliTRDCalibra::ConvertVectorCTTreeHisto(TObjArray *vVectorCT
7068 , const Char_t *name
7069 , const Char_t *nametitle) const
7072 // Convert the vector in a tree with two branchs: the group number
7073 // and the TH1F histo reconstructed from the vector
7076 // Size of the things
7077 Int_t ntotal = (Int_t) pPlaCT->GetEntriesFast();
7079 AliInfo("nothing to write!");
7080 TTree *treeCT = new TTree(name,nametitle);
7084 // Variable of the tree
7085 Int_t groupnumber = -1; // Group calibration
7087 TObjArray vectorCT = *vVectorCT;
7088 TObjArray plaCT = *pPlaCT;
7091 TTree *treeCT = new TTree(name,nametitle);
7092 treeCT->Branch("groupnumber",&groupnumber,"groupnumber/I");
7093 treeCT->Branch("histo","TH1F",&histo,32000,0);
7097 while (k < ntotal) {
7099 groupnumber = ((AliTRDPlace *) plaCT.At(0))->GetPlace();
7100 nome += groupnumber;
7101 histo = ConvertVectorCTHisto(((AliTRDCTInfo *) vectorCT.At(0)),nome);
7103 vectorCT.RemoveAt(0);
7104 vectorCT.Compress();
7114 //_____________________________________________________________________________
7115 TTree *AliTRDCalibra::ConvertVectorPTreeHisto(TObjArray *vVectorP
7117 , const Char_t *name
7118 , const Char_t *nametitle) const
7121 // Convert the vector in a tree with two branchs: the group number
7122 // and the TGraphErrors histo reconstructed from the vector.
7123 // The name must contain "PRF" for PRF calibration and not "PRF"
7124 // for Vdrift calibration
7127 // Size of the things
7128 Int_t ntotal = (Int_t) pPlaP->GetEntriesFast();
7130 AliInfo("nothing to write!");
7131 TTree *treeP = new TTree(name,nametitle);
7135 // Variable of the tree
7136 Int_t groupnumber = -1; // Group calibration
7137 TGraphErrors *histo = 0x0;
7138 TObjArray vectorP = *vVectorP;
7139 TObjArray plaP = *pPlaP;
7142 TTree *treeP = new TTree(name,nametitle);
7143 treeP->Branch("groupnumber",&groupnumber,"groupnumber/I");
7144 treeP->Branch("histo","TGraphErrors",&histo,32000,0);
7148 while (k < ntotal) {
7150 groupnumber = ((AliTRDPlace *) plaP.At(0))->GetPlace();
7151 nome += groupnumber;
7152 histo = ConvertVectorPHisto((AliTRDPInfo *) vectorP.At(0),nome);
7154 vectorP.RemoveAt(0);
7165 //_____________________________________________________________________________
7166 TObjArray *AliTRDCalibra::ConvertTreeVector(TTree *tree) const
7169 // Convert the branch groupnumber of the tree taken from
7170 // TRD.calibration.root in case of vector method in a std::vector
7175 TObjArray *vectorplace = new TObjArray();
7177 // Variable of the tree
7178 Int_t groupnumber = -1; // Group calibration
7181 tree->SetBranchAddress("groupnumber",&groupnumber);
7184 Int_t ntotal = tree->GetEntries();
7185 for (Int_t k = 0; k < ntotal; k++) {
7187 AliTRDPlace *placegroupnumber = new AliTRDPlace();
7188 placegroupnumber->SetPlace(groupnumber);
7189 vectorplace->Add((TObject *) placegroupnumber);
7196 //_____________________________________________________________________________
7197 Bool_t AliTRDCalibra::MergeVectorCT(TObjArray *vVectorCT2, TObjArray *pPlaCT2)
7200 // Add the two vectors and place the result in the first
7203 if (((Int_t) pPlaCT2->GetEntriesFast()) != ((Int_t) vVectorCT2->GetEntriesFast())){
7204 AliInfo("VectorCT2 doesn't correspond to PlaCT2!");
7209 for (Int_t k = 0; k < (Int_t) fPlaCH->GetEntriesFast(); k++) {
7211 // Look if PlaCT1[k] it is also in the second vector
7213 for (Int_t j = 0; j < (Int_t) pPlaCT2->GetEntriesFast(); j++) {
7214 if (((AliTRDPlace *) pPlaCT2->At(j))->GetPlace() ==
7215 ((AliTRDPlace *) fPlaCH->At(k))->GetPlace()) {
7221 // If not in the second vector nothing to do
7223 // If in the second vector
7226 AliTRDCTInfo *fCTInfo = new AliTRDCTInfo();
7227 UShort_t *entries = new UShort_t[fNumberBinCharge];
7229 for (Int_t nu = 0; nu < fNumberBinCharge; nu++) {
7230 entries[nu] = ((AliTRDCTInfo *) fVectorCH->At(((AliTRDPlace *) fPlaCH->At(k))->GetPlace()))->GetEntries()[nu]
7231 + ((AliTRDCTInfo *) vVectorCT2->At(((AliTRDPlace *) fPlaCH->At(k))->GetPlace()))->GetEntries()[nu];
7235 fCTInfo->SetEntries(entries);
7237 // Nothing to do on PlaCT1
7239 // Update the vector
7240 fVectorCH->AddAt((TObject *) fCTInfo,((AliTRDPlace *) fPlaCH->At(k))->GetPlace());
7246 // And at the end the vector in CT2 but not in CH1
7247 for (Int_t k = 0; k < (Int_t) pPlaCT2->GetEntriesFast(); k++) {
7249 // Look if pPlaCT2[k] it is also in the second vector
7251 for (Int_t j = 0; j < (Int_t) fPlaCH->GetEntriesFast(); j++) {
7252 if (((AliTRDPlace *) fPlaCH->At(j))->GetPlace() == ((AliTRDPlace *) pPlaCT2->At(k))->GetPlace()) {
7258 // If not in the first vector
7261 AliTRDCTInfo *fCTInfo = new AliTRDCTInfo();
7262 fCTInfo = ((AliTRDCTInfo *) vVectorCT2->At(((AliTRDPlace *) pPlaCT2->At(k))->GetPlace()));
7265 fPlaCH->Add((TObject *) (pPlaCT2->At(k)));
7266 fVectorCH->Add((TObject *) fCTInfo);
7276 //_____________________________________________________________________________
7277 Bool_t AliTRDCalibra::MergeVectorP(TObjArray *vVectorP2
7282 // Add the two vectors and place the result in the first
7285 if (((Int_t) pPlaP2->GetEntriesFast()) != ((Int_t) vVectorP2->GetEntriesFast())) {
7286 AliInfo("VectorP2 doesn't correspond to PlaP2!");
7293 for (Int_t k = 0; k < (Int_t) fPlaPH->GetEntriesFast(); k++) {
7295 // Look if fPlaPH[k] it is also in the second vector
7297 for (Int_t j = 0; j < (Int_t) pPlaP2->GetEntriesFast(); j++) {
7298 if (((AliTRDPlace *) pPlaP2->At(j))->GetPlace() == ((AliTRDPlace *) fPlaPH->At(k))->GetPlace()) {
7304 // If not in the second vector nothing to do
7306 // If in the second vector
7309 AliTRDPInfo *fPInfo = new AliTRDPInfo();
7310 UShort_t *entries = new UShort_t[fTimeMax];
7311 Float_t *sum = new Float_t[fTimeMax];
7312 Float_t *sumsquare = new Float_t[fTimeMax];
7314 for (Int_t nu = 0; nu < fTimeMax; nu++) {
7316 entries[nu] = ((AliTRDPInfo *) fVectorPH->At(((AliTRDPlace *) fPlaPH->At(k))->GetPlace()))->GetEntries()[nu]
7317 + ((AliTRDPInfo *) vVectorP2->At(((AliTRDPlace *) fPlaPH->At(k))->GetPlace()))->GetEntries()[nu];
7319 Double_t calcul = ((((Double_t) ((AliTRDPInfo *) fVectorPH->At(((AliTRDPlace *) fPlaPH->At(k))->GetPlace()))->GetSum()[nu])
7320 * ((Double_t) ((AliTRDPInfo *) fVectorPH->At(((AliTRDPlace *) fPlaPH->At(k))->GetPlace()))->GetEntries()[nu]))
7321 + (((Double_t) ((AliTRDPInfo *) vVectorP2->At(((AliTRDPlace *) fPlaPH->At(k))->GetPlace()))->GetSum()[nu])
7322 * ((Double_t) ((AliTRDPInfo *) vVectorP2->At(((AliTRDPlace *) fPlaPH->At(k))->GetPlace()))->GetEntries()[nu])))
7323 / ((Double_t) fPInfo->GetEntries()[nu]);
7325 sum[nu] = (Float_t) calcul;
7327 Double_t calculsquare = ((((Double_t) ((AliTRDPInfo *) fVectorPH->At(((AliTRDPlace *) fPlaPH->At(k))->GetPlace()))->GetSumSquare()[nu])
7328 * ((Double_t) ((AliTRDPInfo *) fVectorPH->At(((AliTRDPlace *) fPlaPH->At(k))->GetPlace()))->GetEntries()[nu]))
7329 + (((Double_t) ((AliTRDPInfo *) vVectorP2->At(((AliTRDPlace *) fPlaPH->At(k))->GetPlace()))->GetSumSquare()[nu])
7330 * ((Double_t) ((AliTRDPInfo *) vVectorP2->At(((AliTRDPlace *) fPlaPH->At(k))->GetPlace()))->GetEntries()[nu])))
7331 / ((Double_t) fPInfo->GetEntries()[nu]);
7334 sumsquare[nu] = calculsquare;
7339 fPInfo->SetSum(sum);
7340 fPInfo->SetSumSquare(sumsquare);
7341 fPInfo->SetEntries(entries);
7343 // Nothing to do on PlaCT1
7345 // Update the vector VectorCT1
7346 fVectorPH->AddAt((TObject *) fPInfo,((AliTRDPlace *) fPlaPH->At(k))->GetPlace());
7352 // And at the end the vector in P2 but not in CH1
7353 for (Int_t k = 0; k < (Int_t) pPlaP2->GetEntriesFast(); k++) {
7355 // Look if PlaCT2[k] it is also in the second vector
7357 for (Int_t j = 0; j < (Int_t) fPlaPH->GetEntriesFast(); j++) {
7358 if (((AliTRDPlace *) fPlaPH->At(j))->GetPlace() == ((AliTRDPlace *) pPlaP2->At(k))->GetPlace()) {
7364 // If not in the first vector
7367 AliTRDPInfo *fPInfo = new AliTRDPInfo();
7368 fPInfo = (AliTRDPInfo *) vVectorP2->At(((AliTRDPlace *) pPlaP2->At(k))->GetPlace());
7370 // Add at the end of CH1
7371 fPlaPH->Add(((TObject *) pPlaP2->At(k)));
7372 fVectorPH->Add((TObject *) fPInfo);
7384 for (Int_t k = 0; k < (Int_t) fPlaPRF->GetEntriesFast(); k++) {
7386 // Look if fPlaPRF[k] it is also in the second vector
7388 for (Int_t j = 0; j < (Int_t) pPlaP2->GetEntriesFast(); j++) {
7389 if (((AliTRDPlace *) pPlaP2->At(j))->GetPlace() == ((AliTRDPlace *) fPlaPRF->At(k))->GetPlace()) {
7395 // If not in the second vector nothing to do
7397 // If in the second vector
7400 AliTRDPInfo *fPInfo = new AliTRDPInfo();
7401 UShort_t *entries = new UShort_t[fNumberBinPRF];
7402 Float_t *sum = new Float_t[fNumberBinPRF];
7403 Float_t *sumsquare = new Float_t[fNumberBinPRF];
7405 for (Int_t nu = 0; nu < fNumberBinPRF; nu++) {
7407 entries[nu] = ((AliTRDPInfo *) fVectorPRF->At(((AliTRDPlace *) fPlaPRF->At(k))->GetPlace()))->GetEntries()[nu]
7408 + ((AliTRDPInfo *) vVectorP2->At(((AliTRDPlace *) fPlaPRF->At(k))->GetPlace()))->GetEntries()[nu];
7410 Double_t calcul = ((((Double_t) ((AliTRDPInfo *) fVectorPRF->At(((AliTRDPlace *) fPlaPRF->At(k))->GetPlace()))->GetSum()[nu])
7411 * ((Double_t) ((AliTRDPInfo *) fVectorPRF->At(((AliTRDPlace *) fPlaPRF->At(k))->GetPlace()))->GetEntries()[nu]))
7412 + (((Double_t) ((AliTRDPInfo *) vVectorP2->At(((AliTRDPlace *) fPlaPRF->At(k))->GetPlace()))->GetSum()[nu])
7413 * ((Double_t) ((AliTRDPInfo *) vVectorP2->At(((AliTRDPlace *) fPlaPRF->At(k))->GetPlace()))->GetEntries()[nu])))
7414 / ((Double_t) fPInfo->GetEntries()[nu]);
7416 sum[nu] = (Float_t) calcul;
7418 Double_t calculsquare = ((((Double_t) ((AliTRDPInfo *) fVectorPRF->At(((AliTRDPlace *) fPlaPRF->At(k))->GetPlace()))->GetSumSquare()[nu])
7419 * ((Double_t) ((AliTRDPInfo *) fVectorPRF->At(((AliTRDPlace *) fPlaPRF->At(k))->GetPlace()))->GetEntries()[nu]))
7420 + (((Double_t) ((AliTRDPInfo *) vVectorP2->At(((AliTRDPlace *) fPlaPRF->At(k))->GetPlace()))->GetSumSquare()[nu])
7421 * ((Double_t) ((AliTRDPInfo *) vVectorP2->At(((AliTRDPlace *) fPlaPRF->At(k))->GetPlace()))->GetEntries()[nu])))
7422 / ((Double_t) fPInfo->GetEntries()[nu]);
7424 sumsquare[nu] = calculsquare;
7429 fPInfo->SetSum(sum);
7430 fPInfo->SetSumSquare(sumsquare);
7431 fPInfo->SetEntries(entries);
7433 // Nothing to do on PlaCT1
7435 // Update the vector VectorCT1
7436 fVectorPRF->AddAt((TObject *) fPInfo,((AliTRDPlace *) fPlaPRF->At(k))->GetPlace());
7442 // And at the end the vector in P2 but not in CH1
7443 for (Int_t k = 0; k < (Int_t) pPlaP2->GetEntriesFast(); k++) {
7445 // Look if PlaCT2[k] it is also in the second vector
7447 for (Int_t j = 0; j < (Int_t) fPlaPRF->GetEntriesFast(); j++) {
7448 if (((AliTRDPlace *) fPlaPRF->At(j))->GetPlace() == ((AliTRDPlace *) pPlaP2->At(k))->GetPlace()) {
7454 // If not in the first vector
7457 AliTRDPInfo *fPInfo = new AliTRDPInfo();
7458 fPInfo = (AliTRDPInfo *) vVectorP2->At(((AliTRDPlace *) pPlaP2->At(k))->GetPlace());
7460 // Add at the end of CH1
7461 fPlaPRF->Add(((TObject *) pPlaP2->At(k)));
7462 fVectorPRF->Add((TObject *) fPInfo);
7474 //____________Fit Methods______________________________________________________
7476 //_____________________________________________________________________________
7477 void AliTRDCalibra::FitPente(TH1* projPH, Int_t idect)
7480 // Slope methode for the drift velocity
7484 const Float_t kDrWidth = AliTRDgeometry::DrThick();
7491 Double_t vdriftCoefE = 0.0;
7492 Double_t t0CoefE = 0.0;
7493 fVdriftCoef[1] = 0.0;
7495 TLine *line = new TLine();
7498 TAxis *xpph = projPH->GetXaxis();
7499 Int_t nbins = xpph->GetNbins();
7500 Double_t lowedge = xpph->GetBinLowEdge(1);
7501 Double_t upedge = xpph->GetBinUpEdge(xpph->GetNbins());
7502 Double_t widbins = (upedge - lowedge) / nbins;
7503 Double_t limit = upedge + 0.5 * widbins;
7505 // Beginning of the signal
7506 TH1D *pentea = new TH1D("pentea","pentea",projPH->GetNbinsX(),0,(Float_t) limit);
7507 for (Int_t k = 1; k < projPH->GetNbinsX(); k++) {
7508 pentea->SetBinContent(k,(Double_t) (projPH->GetBinContent(k+1) - projPH->GetBinContent(k)));
7511 binmax = (Int_t) pentea->GetMaximumBin();
7514 AliInfo("Put the binmax from 1 to 2 to enable the fit");
7516 pentea->Fit("pol2","0MR","",TMath::Max(pentea->GetBinCenter(binmax-1),0.0),pentea->GetBinCenter(binmax+1));
7517 Float_t l3P1am = pentea->GetFunction("pol2")->GetParameter(1);
7518 Float_t l3P2am = pentea->GetFunction("pol2")->GetParameter(2);
7519 Float_t l3P1amE = pentea->GetFunction("pol2")->GetParError(1);
7520 Float_t l3P2amE = pentea->GetFunction("pol2")->GetParError(2);
7522 fPhd[0] = -(l3P1am / (2 * l3P2am));
7524 if((l3P1am != 0.0) && (l3P2am != 0.0)){
7525 t0CoefE = (l3P1amE/l3P1am + l3P2amE/l3P2am)*fPhd[0];
7528 // Amplification region
7531 for (Int_t kbin = 1; kbin < projPH->GetNbinsX(); kbin ++) {
7532 if (((projPH->GetBinContent(kbin+1) - projPH->GetBinContent(kbin)) <= 0.0) && (ju == 0) && (kbin > (fPhd[0]/widbins))) {
7539 AliInfo("Put the binmax from 1 to 2 to enable the fit");
7541 projPH->Fit("pol2","0MR","",TMath::Max(projPH->GetBinCenter(binmax-1),0.0),projPH->GetBinCenter(binmax+1));
7542 Float_t l3P1amf = projPH->GetFunction("pol2")->GetParameter(1);
7543 Float_t l3P2amf = projPH->GetFunction("pol2")->GetParameter(2);
7544 Float_t l3P1amfE = projPH->GetFunction("pol2")->GetParError(1);
7545 Float_t l3P2amfE = projPH->GetFunction("pol2")->GetParError(2);
7548 fPhd[1] = -(l3P1amf / (2 * l3P2amf));
7550 if((l3P1amf != 0.0) && (l3P2amf != 0.0)){
7551 vdriftCoefE = (l3P1amfE/l3P1amf + l3P2amfE/l3P2amf)*fPhd[1];
7555 TH1D *pente = new TH1D("pente","pente",projPH->GetNbinsX(),0,(Float_t) limit);
7556 for (Int_t k = binmax+4; k < projPH->GetNbinsX(); k++) {
7557 pente->SetBinContent(k,(Double_t) (projPH->GetBinContent(k+1) - projPH->GetBinContent(k)));
7559 binmin = (Int_t) pente->GetMinimumBin();
7562 AliInfo("Put the binmax from 1 to 2 to enable the fit");
7567 ,TMath::Max(pente->GetBinCenter(binmin-1), 0.0)
7568 ,TMath::Min(pente->GetBinCenter(binmin+2),(Double_t) limit));
7569 Float_t l3P1dr = pente->GetFunction("pol2")->GetParameter(1);
7570 Float_t l3P2dr = pente->GetFunction("pol2")->GetParameter(2);
7571 Float_t l3P1drE = pente->GetFunction("pol2")->GetParError(1);
7572 Float_t l3P2drE = pente->GetFunction("pol2")->GetParError(2);
7574 fPhd[2] = -(l3P1dr / (2 * l3P2dr));
7576 if((l3P1dr != 0.0) && (l3P2dr != 0.0)){
7577 vdriftCoefE += (l3P1drE/l3P1dr + l3P2drE/l3P2dr)*fPhd[2];
7580 if ((fPhd[2] > fPhd[0]) &&
7581 (fPhd[2] > fPhd[1]) &&
7582 (fPhd[1] > fPhd[0])) {
7583 fVdriftCoef[1] = (kDrWidth) / (fPhd[2]-fPhd[1]);
7584 if (fPhd[0] >= 0.0) {
7585 fT0Coef[1] = (fPhd[0] - fT0Shift) / widbins;
7586 if (fT0Coef[1] < -1.0) {
7587 fT0Coef[1] = fT0Coef[2];
7591 fT0Coef[1] = fT0Coef[2];
7595 fVdriftCoef[1] = -TMath::Abs(fVdriftCoef[2]);
7596 fT0Coef[1] = fT0Coef[2];
7599 if ((fDebug == 1) ||
7601 fCoefVdrift[1][idect] = fVdriftCoef[1];
7602 fCoefVdriftE[1] [idect] = vdriftCoefE;
7603 fCoefT0[1][idect] = fT0Coef[1];
7604 fCoefT0E[1][idect] = t0CoefE;
7608 TCanvas *cpentei = new TCanvas("cpentei","cpentei",50,50,600,800);
7611 line->SetLineColor(2);
7612 line->DrawLine(fPhd[0],0,fPhd[0],projPH->GetMaximum());
7613 line->DrawLine(fPhd[1],0,fPhd[1],projPH->GetMaximum());
7614 line->DrawLine(fPhd[2],0,fPhd[2],projPH->GetMaximum());
7615 AliInfo(Form("fPhd[0] (beginning of the signal): %f" ,(Float_t) fPhd[0]));
7616 AliInfo(Form("fPhd[1] (end of the amplification region): %f" ,(Float_t) fPhd[1]));
7617 AliInfo(Form("fPhd[2] (end of the drift region): %f" ,(Float_t) fPhd[2]));
7618 AliInfo(Form("fVriftCoef[1] (with only the drift region(default)): %f",(Float_t) fVdriftCoef[1]));
7630 //_____________________________________________________________________________
7631 void AliTRDCalibra::FitPH(TH1* projPH, Int_t idect)
7634 // Fit methode for the drift velocity
7638 const Float_t kDrWidth = AliTRDgeometry::DrThick();
7641 TAxis *xpph = projPH->GetXaxis();
7642 Double_t upedge = xpph->GetBinUpEdge(xpph->GetNbins());
7644 TF1 *fPH = new TF1("fPH",AliTRDCalibra::PH,-0.05,3.2,6);
7645 fPH->SetParameter(0,0.469); // Scaling
7646 fPH->SetParameter(1,0.18); // Start
7647 fPH->SetParameter(2,0.0857325); // AR
7648 fPH->SetParameter(3,1.89); // DR
7649 fPH->SetParameter(4,0.08); // QA/QD
7650 fPH->SetParameter(5,0.0); // Baseline
7652 TLine *line = new TLine();
7654 fVdriftCoef[0] = 0.0;
7656 Double_t vdriftCoefE = 0.0;
7657 Double_t t0CoefE = 0.0;
7659 if (idect%fFitPHPeriode == 0) {
7661 AliInfo(Form("<AliTRDCalibra::FitPH> The detector %d will be fitted",idect));
7662 fPH->SetParameter(0,(projPH->Integral()-(projPH->GetBinContent(1)*projPH->GetNbinsX())) * 0.00028); // Scaling
7663 fPH->SetParameter(1,fPhd[0] - 0.1); // Start
7664 fPH->SetParameter(2,fPhd[1] - fPhd[0]); // AR
7665 fPH->SetParameter(3,fPhd[2] - fPhd[1]); // DR
7666 fPH->SetParameter(4,0.225); // QA/QD
7667 fPH->SetParameter(5,(Float_t) projPH->GetBinContent(1));
7670 projPH->Fit(fPH,"0M","",0.0,upedge);
7674 TCanvas *cpente = new TCanvas("cpente","cpente",50,50,600,800);
7676 projPH->Fit(fPH,"M+","",0.0,upedge);
7678 line->SetLineColor(4);
7679 line->DrawLine(fPH->GetParameter(1)
7681 ,fPH->GetParameter(1)
7682 ,projPH->GetMaximum());
7683 line->DrawLine(fPH->GetParameter(1)+fPH->GetParameter(2)
7685 ,fPH->GetParameter(1)+fPH->GetParameter(2)
7686 ,projPH->GetMaximum());
7687 line->DrawLine(fPH->GetParameter(1)+fPH->GetParameter(2)+fPH->GetParameter(3)
7689 ,fPH->GetParameter(1)+fPH->GetParameter(2)+fPH->GetParameter(3)
7690 ,projPH->GetMaximum());
7693 if (fPH->GetParameter(3) != 0) {
7694 fVdriftCoef[0] = kDrWidth / (fPH->GetParameter(3));
7695 vdriftCoefE = (fPH->GetParError(3)/fPH->GetParameter(3))*fVdriftCoef[0];
7696 fT0Coef[0] = fPH->GetParameter(1);
7697 t0CoefE = fPH->GetParError(1);
7700 fVdriftCoef[0] = -TMath::Abs(fVdriftCoef[2]);
7701 fT0Coef[0] = fT0Coef[2];
7704 if ((fDebug == 1) ||
7706 fCoefVdrift[0][idect] = fVdriftCoef[0];
7707 fCoefVdriftE[0][idect] = vdriftCoefE;
7708 fCoefT0[0][idect] = fT0Coef[0];
7709 fCoefT0E[0][idect] = t0CoefE;
7712 AliInfo(Form("fVdriftCoef[0]: %f",(Float_t) fVdriftCoef[0]));
7719 // Put the default value
7720 if ((fDebug <= 1) ||
7722 fCoefVdrift[0][idect] = -TMath::Abs(fVdriftCoef[2]);
7723 fCoefT0[0][idect] = -TMath::Abs(fT0Coef[2]);
7734 //_____________________________________________________________________________
7735 void AliTRDCalibra::FitPRF(TH1 *projPRF, Int_t idect)
7738 // Fit methode for the sigma of the pad response function
7742 Double_t prfCoefE = 0.0;
7745 projPRF->Fit("gaus","0M","",-fRangeFitPRF,fRangeFitPRF);
7749 TCanvas *cfit = new TCanvas("cfit","cfit",50,50,600,800);
7751 projPRF->Fit("gaus","M+","",-fRangeFitPRF,fRangeFitPRF);
7755 fPRFCoef[0] = projPRF->GetFunction("gaus")->GetParameter(2);
7756 prfCoefE = projPRF->GetFunction("gaus")->GetParError(2);
7757 if(fPRFCoef[0] <= 0.0) fPRFCoef[0] = -fPRFCoef[1];
7759 if ((fDebug == 1) ||
7761 fCoefPRF[0][idect] = fPRFCoef[0];
7762 fCoefPRFE[idect] = prfCoefE;
7765 AliInfo(Form("fPRFCoef[0]: %f",(Float_t) fPRFCoef[0]));
7770 //_____________________________________________________________________________
7771 void AliTRDCalibra::FitCH(TH1 *projch, Int_t idect)
7774 // Fit methode for the gain factor
7777 fChargeCoef[0] = 0.0;
7778 fChargeCoef[1] = 0.0;
7779 Double_t chargeCoefE0 = 0.0;
7780 Double_t chargeCoefE1 = 0.0;
7781 TF1 *fLandauGaus = new TF1("fLandauGaus",FuncLandauGaus,0,300,5);
7783 fChargeCoef[1] = projch->GetMean();
7784 chargeCoefE1 = projch->GetMeanError();
7785 projch->Fit("landau","0",""
7786 ,(Float_t) fChargeCoef[1]/fBeginFitCharge
7787 ,projch->GetBinCenter(projch->GetNbinsX()));
7788 fL3P0 = projch->GetFunction("landau")->GetParameter(0);
7789 Double_t l3P1 = projch->GetFunction("landau")->GetParameter(1);
7790 fL3P2 = projch->GetFunction("landau")->GetParameter(2);
7792 projch->Fit("gaus","0",""
7793 ,(Float_t) fChargeCoef[1]/fBeginFitCharge
7794 ,projch->GetBinCenter(projch->GetNbinsX()));
7795 Double_t g3P0 = projch->GetFunction("gaus")->GetParameter(0);
7796 fG3P2 = projch->GetFunction("gaus")->GetParameter(2);
7798 fLandauGaus->SetParameters(fL3P0,l3P1,fL3P2,g3P0,fG3P2);
7799 if ((fDebug <= 1) ||
7801 projch->Fit("fLandauGaus","0",""
7802 ,(Float_t) fChargeCoef[1]/fBeginFitCharge
7803 ,projch->GetBinCenter(projch->GetNbinsX()));
7807 TCanvas *cp = new TCanvas("cp","cp",50,50,600,800);
7809 projch->Fit("fLandauGaus","+",""
7810 ,(Float_t) fChargeCoef[1]/fBeginFitCharge
7811 ,projch->GetBinCenter(projch->GetNbinsX()));
7813 fLandauGaus->Draw("same");
7816 if (projch->GetFunction("fLandauGaus")->GetParameter(1) > 0) {
7817 // Calcul of "real" coef
7818 CalculChargeCoefMean(fCountDet[0],(Int_t) idect,kTRUE);
7819 fChargeCoef[0] = projch->GetFunction("fLandauGaus")->GetParameter(1);
7820 chargeCoefE0 = projch->GetFunction("fLandauGaus")->GetParError(1);
7823 // Calcul of "real" coef
7824 CalculChargeCoefMean(fCountDet[0],(Int_t) idect,kFALSE);
7825 fChargeCoef[0] = -TMath::Abs(fChargeCoef[3]);
7829 AliInfo(Form("fChargeCoef[0]: %f",(Float_t) fChargeCoef[0]));
7830 AliInfo(Form("fChargeCoef[1]: %f",(Float_t) fChargeCoef[1]));
7833 if ((fDebug == 1) ||
7835 if (fChargeCoef[0] > 0.0) {
7836 fCoefCharge[0][idect]= fChargeCoef[0];
7837 fCoefChargeE[0][idect]= chargeCoefE0;
7838 fCoefCharge[1][idect]= fChargeCoef[1];
7839 fCoefChargeE[1][idect]= chargeCoefE1;
7842 fL3P0 = fLandauGaus->Integral(0.3*projch->GetMean(),3*projch->GetMean());
7843 fG3P2 = fLandauGaus->GetParameter(2);
7844 fL3P2 = fLandauGaus->GetParameter(4);
7852 //_____________________________________________________________________________
7853 void AliTRDCalibra::FitBisCH(TH1* projch, Int_t idect)
7856 // Fit methode for the gain factor more time consuming
7859 // Setting fit range and start values
7861 //Double_t sv[4] = { l3P2, fChargeCoef[1], projch->Integral("width"), fG3P2 };
7862 Double_t sv[4] = { fL3P2, fChargeCoef[1], fL3P0, fG3P2 };
7863 Double_t pllo[4] = { 0.001, 0.001, 0.001, 0.001 };
7864 Double_t plhi[4] = { 300.0, 300.0, 100000000.0, 300.0 };
7865 Double_t fp[4] = { 1.0, 1.0, 1.0, 1.0 };
7866 Double_t fpe[4] = { 1.0, 1.0, 1.0, 1.0 };
7867 fr[0] = 0.3 * projch->GetMean();
7868 fr[1] = 3.0 * projch->GetMean();
7869 fChargeCoef[2] = 0.0;
7870 Double_t chargeCoefE2 = 0.0;
7874 TF1 *fitsnr = LanGauFit(projch,&fr[0],&sv[0]
7879 Double_t projchPeak;
7880 Double_t projchFWHM;
7881 LanGauPro(fp,projchPeak,projchFWHM);
7884 fChargeCoef[2] = fp[1];
7885 chargeCoefE2 = fpe[1];
7886 //chargeCoefE2 = chisqr;
7889 fChargeCoef[2] = -TMath::Abs(fChargeCoef[3]);
7893 AliInfo(Form("fChargeCoef[2]: %f",(Float_t) fChargeCoef[2]));
7894 TCanvas *cpy = new TCanvas("cpy","cpy",50,50,600,800);
7897 fitsnr->Draw("same");
7900 if ((fDebug == 1) ||
7902 if (fChargeCoef[2] > 0.0) {
7903 fCoefCharge[2][idect]= fChargeCoef[2];
7904 fCoefChargeE[2][idect]= chargeCoefE2;
7914 //_____________________________________________________________________________
7915 void AliTRDCalibra::NormierungCharge()
7918 // Normalisation of the gain factor resulting for the fits
7921 // Calcul of the mean of the fit
7923 for (Int_t k = 0; k < (Int_t) fVectorFitCH->GetEntriesFast(); k++) {
7925 Int_t detector = ((AliTRDFitCHInfo *) fVectorFitCH->At(k))->GetDetector();
7926 Float_t *coef = ((AliTRDFitCHInfo *) fVectorFitCH->At(k))->GetCoef();
7927 if (GetChamber(detector) == 2) {
7930 if (GetChamber(detector) != 2) {
7933 for (Int_t j = 0; j < total; j++) {
7941 fScaleFitFactor = fScaleFitFactor / sum;
7944 fScaleFitFactor = 1.0;
7947 if ((fDebug == 3) ||
7949 if ((fCoefChargeDB[0]->GetEntries() > 0.0) &&
7950 (fCoefChargeDB[0]->GetSumOfWeights() > 0.0)) {
7951 fCoefChargeDB[0]->Scale(fCoefChargeDB[0]->GetEntries() / fCoefChargeDB[0]->GetSumOfWeights());
7953 if ((fMeanChargeOn) &&
7954 (fCoefChargeDB[1]->GetEntries() > 0.0) &&
7955 (fCoefChargeDB[1]->GetSumOfWeights() > 0.0)) {
7956 fCoefChargeDB[1]->Scale(fCoefChargeDB[1]->GetEntries() / fCoefChargeDB[1]->GetSumOfWeights());
7958 if ((fFitChargeBisOn) &&
7959 (fCoefChargeDB[2]->GetEntries() > 0.0) &&
7960 (fCoefChargeDB[2]->GetSumOfWeights() > 0.0)) {
7961 fCoefChargeDB[2]->Scale(fCoefChargeDB[2]->GetEntries() / fCoefChargeDB[2]->GetSumOfWeights());
7967 //_____________________________________________________________________________
7968 TH1I *AliTRDCalibra::ReBin(TH1I *hist) const
7971 // Rebin of the 1D histo for the gain calibration if needed.
7972 // you have to choose fRebin, divider of fNumberBinCharge
7975 TAxis *xhist = hist->GetXaxis();
7976 TH1I *rehist = new TH1I("projrebin","",(Int_t) xhist->GetNbins()/fRebin
7977 ,xhist->GetBinLowEdge(1)
7978 ,xhist->GetBinUpEdge(xhist->GetNbins()));
7980 AliInfo(Form("fRebin: %d",fRebin));
7982 for (Int_t k = 1; k <= (Int_t) xhist->GetNbins()/fRebin; k++) {
7984 for (Int_t ji = i; ji < i+fRebin; ji++) {
7985 sum += hist->GetBinContent(ji);
7988 rehist->SetBinContent(k,sum);
7993 TCanvas *crebin = new TCanvas("crebin","",50,50,600,800);
8002 //_____________________________________________________________________________
8003 TH1F *AliTRDCalibra::ReBin(TH1F *hist) const
8006 // Rebin of the 1D histo for the gain calibration if needed
8007 // you have to choose fRebin divider of fNumberBinCharge
8010 TAxis *xhist = hist->GetXaxis();
8011 TH1F *rehist = new TH1F("projrebin","",(Int_t) xhist->GetNbins()/fRebin
8012 ,xhist->GetBinLowEdge(1)
8013 ,xhist->GetBinUpEdge(xhist->GetNbins()));
8015 AliInfo(Form("fRebin: %d",fRebin));
8017 for (Int_t k = 1; k <= (Int_t) xhist->GetNbins()/fRebin; k++) {
8019 for (Int_t ji = i; ji < i+fRebin; ji++) {
8020 sum += hist->GetBinContent(ji);
8023 rehist->SetBinContent(k,sum);
8028 TCanvas *crebin = new TCanvas("crebin","",50,50,600,800);
8037 //_____________________________________________________________________________
8038 TH1F *AliTRDCalibra::CorrectTheError(TGraphErrors *hist)
8041 // In the case of the vectors method the trees contains TGraphErrors for PH and PRF
8042 // to be able to add them after
8043 // We convert it to a TH1F to be able to applied the same fit function method
8044 // After having called this function you can not add the statistics anymore
8049 Int_t nbins = hist->GetN();
8050 Double_t *x = hist->GetX();
8051 Double_t *entries = hist->GetEX();
8052 Double_t *mean = hist->GetY();
8053 Double_t *square = hist->GetEY();
8054 fEntriesCurrent = 0;
8060 Double_t step = x[1] - x[0];
8061 Double_t minvalue = x[0] - step/2;
8062 Double_t maxvalue = x[(nbins-1)] + step/2;
8064 rehist = new TH1F("projcorrecterror","",nbins,minvalue,maxvalue);
8066 for (Int_t k = 0; k < nbins; k++) {
8067 rehist->SetBinContent(k+1,mean[k]);
8068 if (entries[k] > 0.0) {
8069 fEntriesCurrent += (Int_t) entries[k];
8070 Double_t d = TMath::Abs(square[k] - (mean[k]*mean[k]));
8071 rehist->SetBinError(k+1,TMath::Sqrt(d/entries[k]));
8074 rehist->SetBinError(k+1,0.0);
8082 //_____________________________________________________________________________
8083 TGraphErrors *AliTRDCalibra::AddProfiles(TGraphErrors *hist1
8084 , TGraphErrors *hist2) const
8087 // In the case of the vectors method we use TGraphErrors for PH and PRF
8088 // to be able to add the them after
8089 // Here we add the TGraphErrors
8092 // First TGraphErrors
8093 Int_t nbins1 = hist1->GetN();
8094 Double_t *x1 = hist1->GetX();
8095 Double_t *ex1 = hist1->GetEX();
8096 Double_t *y1 = hist1->GetY();
8097 Double_t *ey1 = hist1->GetEY();
8099 TGraphErrors *rehist = new TGraphErrors(nbins1);
8101 // Second TGraphErrors
8102 Double_t *ex2 = hist2->GetEX();
8103 Double_t *y2 = hist2->GetY();
8104 Double_t *ey2 = hist2->GetEY();
8106 // Define the Variables for the new TGraphErrors
8112 for (Int_t k = 0; k < nbins1; k++) {
8113 Double_t nentries = 0.0;
8118 if ((ex2[k] == 0.0) &&
8122 if ((ex2[k] == 0.0) &&
8129 if ((ex2[k] > 0.0) &&
8136 if ((ex2[k] > 0.0) &&
8138 nentries = ex1[k] + ex2[k];
8139 y = ( y1[k]*ex1[k]+ y2[k]*ex2[k]) / nentries;
8140 ey = (ey1[k]*ex1[k]+ey2[k]*ex2[k]) / nentries;
8143 rehist->SetPoint(k,x,y);
8144 rehist->SetPointError(k,ex,ey);
8152 //____________Some basic geometry function_____________________________________
8155 //_____________________________________________________________________________
8156 Int_t AliTRDCalibra::GetPlane(Int_t d) const
8159 // Reconstruct the plane number from the detector number
8162 return ((Int_t) (d % 6));
8166 //_____________________________________________________________________________
8167 Int_t AliTRDCalibra::GetChamber(Int_t d) const
8170 // Reconstruct the chamber number from the detector number
8174 return ((Int_t) (d % 30) / fgkNplan);
8178 //_____________________________________________________________________________
8179 Int_t AliTRDCalibra::GetSector(Int_t d) const
8182 // Reconstruct the sector number from the detector number
8186 return ((Int_t) (d / fg));
8191 //____________Fill and Init tree Gain, PRF, Vdrift and T0______________________
8194 //_____________________________________________________________________________
8195 void AliTRDCalibra::InitTreePRF()
8198 // Init the tree where the coefficients from the fit methods can be stored
8202 fPRFPad = new Float_t[2304];
8203 fPRF = new TTree("PRF","PRF");
8204 fPRF->Branch("detector",&fPRFDetector,"detector/I");
8205 fPRF->Branch("width" ,fPRFPad ,"width[2304]/F");
8207 // Set to default value for the plane 0 supposed to be the first one
8208 for (Int_t k = 0; k < 2304; k++) {
8215 //_____________________________________________________________________________
8216 void AliTRDCalibra::FillTreePRF(Int_t countdet)
8219 // Fill the tree with the sigma of the pad response function for the detector countdet
8222 Int_t numberofgroup = 0;
8223 fPRFDetector = countdet;
8226 if (GetChamber((Int_t)(countdet+1)) == 2) {
8227 numberofgroup = 1728;
8230 numberofgroup = 2304;
8233 // Reset to default value for the next
8234 for (Int_t k = 0; k < numberofgroup; k++) {
8235 if (GetPlane((Int_t) (countdet+1)) == 0) {
8238 if (GetPlane((Int_t) (countdet+1)) == 1) {
8241 if (GetPlane((Int_t) (countdet+1)) == 2) {
8244 if (GetPlane((Int_t) (countdet+1)) == 3) {
8247 if (GetPlane((Int_t) (countdet+1)) == 4) {
8250 if (GetPlane((Int_t) (countdet+1)) == 5) {
8259 //_____________________________________________________________________________
8260 void AliTRDCalibra::ConvertVectorFitCHTree()
8263 // Convert the vector stuff to a tree of 1D histos if the user
8264 // want to write it after the fill functions
8267 Int_t detector = -1;
8268 Int_t numberofgroup = 1;
8269 Float_t gainPad[2304];
8271 fGain = new TTree("Gain","Gain");
8272 fGain->Branch("detector",&detector,"detector/I");
8273 fGain->Branch("gainPad" ,gainPad ,"gainPad[2304]/F");
8275 Int_t loop = (Int_t) fVectorFitCH->GetEntriesFast();
8276 for (Int_t k = 0; k < loop; k++) {
8277 detector = ((AliTRDFitCHInfo *) fVectorFitCH->At(k))->GetDetector();
8278 if (GetChamber((Int_t) ((AliTRDFitCHInfo *) fVectorFitCH->At(k))->GetDetector()) == 2) {
8279 numberofgroup = 1728;
8282 numberofgroup = 2304;
8284 for (Int_t i = 0; i < numberofgroup; i++) {
8285 if (((AliTRDFitCHInfo *) fVectorFitCH->At(k))->GetCoef()[i] >= 0) {
8286 gainPad[i] = ((AliTRDFitCHInfo *) fVectorFitCH->At(k))->GetCoef()[i] * fScaleFitFactor;
8289 gainPad[i] = (Float_t) ((AliTRDFitCHInfo *) fVectorFitCH->At(k))->GetCoef()[i];
8297 //_____________________________________________________________________________
8298 void AliTRDCalibra::FillTreeVdrift(Int_t countdet)
8301 // Fill the tree with the drift velocities for the detector countdet
8304 Int_t numberofgroup = 0;
8305 fVdriftDetector = countdet;
8308 if (GetChamber((Int_t)(countdet+1)) == 2) {
8309 numberofgroup = 1728;
8312 numberofgroup = 2304;
8314 // Reset to default value the gain coef
8315 for (Int_t k = 0; k < numberofgroup; k++) {
8316 fVdriftPad[k] = -1.5;
8318 fVdriftDetector = -1;
8322 //_____________________________________________________________________________
8323 void AliTRDCalibra::InitTreePH()
8326 // Init the tree where the coefficients from the fit methods can be stored
8330 fVdriftPad = new Float_t[2304];
8331 fVdrift = new TTree("Vdrift","Vdrift");
8332 fVdrift->Branch("detector",&fVdriftDetector,"detector/I");
8333 fVdrift->Branch("vdrift" ,fVdriftPad ,"vdrift[2304]/F");
8334 // Set to default value for the plane 0 supposed to be the first one
8335 for (Int_t k = 0; k < 2304; k++) {
8336 fVdriftPad[k] = -1.5;
8338 fVdriftDetector = -1;
8342 //_____________________________________________________________________________
8343 void AliTRDCalibra::FillTreeT0(Int_t countdet)
8346 // Fill the tree with the t0 value for the detector countdet
8349 Int_t numberofgroup = 0;
8351 fT0Detector = countdet;
8354 if (GetChamber((Int_t) (countdet+1)) == 2) {
8355 numberofgroup = 1728;
8358 numberofgroup = 2304;
8360 // Reset to default value
8361 for (Int_t k = 0; k < numberofgroup; k++) {
8368 //_____________________________________________________________________________
8369 void AliTRDCalibra::InitTreeT0()
8372 // Init the tree where the coefficients from the fit methods can be stored
8376 fT0Pad = new Float_t[2304];
8377 fT0 = new TTree("T0","T0");
8378 fT0->Branch("detector",&fT0Detector,"detector/I");
8379 fT0->Branch("t0",fT0Pad,"t0[2304]/F");
8380 //Set to default value for the plane 0 supposed to be the first one
8381 for(Int_t k = 0; k < 2304; k++){
8389 //____________Private Functions________________________________________________
8392 //_____________________________________________________________________________
8393 Double_t AliTRDCalibra::PH(Double_t *x, Double_t *par)
8396 // Function for the fit
8399 //TF1 *fAsymmGauss = new TF1("fAsymmGauss",AsymmGauss,0,4,6);
8401 //PARAMETERS FOR FIT PH
8403 //fAsymmGauss->SetParameter(0,0.113755);
8404 //fAsymmGauss->SetParameter(1,0.350706);
8405 //fAsymmGauss->SetParameter(2,0.0604244);
8406 //fAsymmGauss->SetParameter(3,7.65596);
8407 //fAsymmGauss->SetParameter(4,1.00124);
8408 //fAsymmGauss->SetParameter(5,0.870597); // No tail cancelation
8416 Double_t dx = 0.005;
8417 Double_t xs = par[1];
8419 Double_t paras[2] = { 0.0, 0.0 };
8422 if ((xs >= par[1]) &&
8423 (xs < (par[1]+par[2]))) {
8424 //fAsymmGauss->SetParameter(0,par[0]);
8425 //fAsymmGauss->SetParameter(1,xs);
8426 //ss += fAsymmGauss->Eval(xx);
8429 ss += AsymmGauss(&xx,paras);
8431 if ((xs >= (par[1]+par[2])) &&
8432 (xs < (par[1]+par[2]+par[3]))) {
8433 //fAsymmGauss->SetParameter(0,par[0]*par[4]);
8434 //fAsymmGauss->SetParameter(1,xs);
8435 //ss += fAsymmGauss->Eval(xx);
8436 paras[0] = par[0]*par[4];
8438 ss += AsymmGauss(&xx,paras);
8447 //_____________________________________________________________________________
8448 Double_t AliTRDCalibra::AsymmGauss(Double_t *x, Double_t *par)
8451 // Function for the fit
8454 //par[0] = normalization
8462 Double_t par1save = par[1];
8463 //Double_t par2save = par[2];
8464 Double_t par2save = 0.0604244;
8465 //Double_t par3save = par[3];
8466 Double_t par3save = 7.65596;
8467 //Double_t par5save = par[5];
8468 Double_t par5save = 0.870597;
8469 Double_t dx = x[0] - par1save;
8471 Double_t sigma2 = par2save*par2save;
8472 Double_t sqrt2 = TMath::Sqrt(2.0);
8473 Double_t exp1 = par3save * TMath::Exp(-par3save * (dx - 0.5 * par3save * sigma2))
8474 * (1.0 - TMath::Erf((par3save * sigma2 - dx) / (sqrt2 * par2save)));
8475 Double_t exp2 = par5save * TMath::Exp(-par5save * (dx - 0.5 * par5save * sigma2))
8476 * (1.0 - TMath::Erf((par5save * sigma2 - dx) / (sqrt2 * par2save)));
8478 //return par[0]*(exp1+par[4]*exp2);
8479 return par[0] * (exp1 + 1.00124 * exp2);
8483 //_____________________________________________________________________________
8484 Double_t AliTRDCalibra::FuncLandauGaus(Double_t *x, Double_t *par)
8487 // Sum Landau + Gaus with identical mean
8490 Double_t valLandau = par[0] * TMath::Landau(x[0],par[1],par[2]);
8491 //Double_t valGaus = par[3] * TMath::Gaus(x[0],par[4],par[5]);
8492 Double_t valGaus = par[3] * TMath::Gaus(x[0],par[1],par[4]);
8493 Double_t val = valLandau + valGaus;
8499 //_____________________________________________________________________________
8500 Double_t AliTRDCalibra::LanGauFun(Double_t *x, Double_t *par)
8503 // Function for the fit
8506 // par[0]=Width (scale) parameter of Landau density
8507 // par[1]=Most Probable (MP, location) parameter of Landau density
8508 // par[2]=Total area (integral -inf to inf, normalization constant)
8509 // par[3]=Width (sigma) of convoluted Gaussian function
8511 // In the Landau distribution (represented by the CERNLIB approximation),
8512 // the maximum is located at x=-0.22278298 with the location parameter=0.
8513 // This shift is corrected within this function, so that the actual
8514 // maximum is identical to the MP parameter.
8517 // Numeric constants
8518 Double_t invsq2pi = 0.3989422804014; // (2 pi)^(-1/2)
8519 Double_t mpshift = -0.22278298; // Landau maximum location
8521 // Control constants
8522 Double_t np = 100.0; // Number of convolution steps
8523 Double_t sc = 5.0; // Convolution extends to +-sc Gaussian sigmas
8535 // MP shift correction
8536 mpc = par[1] - mpshift * par[0];
8538 // Range of convolution integral
8539 xlow = x[0] - sc * par[3];
8540 xupp = x[0] + sc * par[3];
8542 step = (xupp - xlow) / np;
8544 // Convolution integral of Landau and Gaussian by sum
8545 for (i = 1.0; i <= np/2; i++) {
8547 xx = xlow + (i-.5) * step;
8548 fland = TMath::Landau(xx,mpc,par[0]) / par[0];
8549 sum += fland * TMath::Gaus(x[0],xx,par[3]);
8551 xx = xupp - (i-.5) * step;
8552 fland = TMath::Landau(xx,mpc,par[0]) / par[0];
8553 sum += fland * TMath::Gaus(x[0],xx,par[3]);
8557 return (par[2] * step * sum * invsq2pi / par[3]);
8561 //_____________________________________________________________________________
8562 TF1 *AliTRDCalibra::LanGauFit(TH1 *his, Double_t *fitrange, Double_t *startvalues
8563 , Double_t *parlimitslo, Double_t *parlimitshi
8564 , Double_t *fitparams, Double_t *fiterrors
8565 , Double_t *chiSqr, Int_t *ndf)
8568 // Function for the fit
8572 Char_t funname[100];
8574 AliInfo(Form(funname,"Fitfcn_%s",his->GetName()));
8576 TF1 *ffitold = (TF1 *) gROOT->GetListOfFunctions()->FindObject(funname);
8581 TF1 *ffit = new TF1(funname,LanGauFun,fitrange[0],fitrange[1],4);
8582 ffit->SetParameters(startvalues);
8583 ffit->SetParNames("Width","MP","Area","GSigma");
8585 for (i = 0; i < 4; i++) {
8586 ffit->SetParLimits(i,parlimitslo[i],parlimitshi[i]);
8589 his->Fit(funname,"RB0"); // Fit within specified range, use ParLimits, do not plot
8591 ffit->GetParameters(fitparams); // Obtain fit parameters
8592 for (i = 0; i < 4; i++) {
8593 fiterrors[i] = ffit->GetParError(i); // Obtain fit parameter errors
8595 chiSqr[0] = ffit->GetChisquare(); // Obtain chi^2
8596 ndf[0] = ffit->GetNDF(); // Obtain ndf
8598 return (ffit); // Return fit function
8602 //_____________________________________________________________________________
8603 Int_t AliTRDCalibra::LanGauPro(Double_t *params, Double_t &maxx, Double_t &fwhm)
8606 // Function for the fit
8619 Int_t maxcalls = 10000;
8621 // Search for maximum
8622 p = params[1] - 0.1 * params[0];
8623 step = 0.05 * params[0];
8627 while ((l != lold) && (i < maxcalls)) {
8631 l = LanGauFun(&x,params);
8633 step = -step / 10.0;
8638 if (i == maxcalls) {
8644 // Search for right x location of fy
8645 p = maxx + params[0];
8651 while ( (l != lold) && (i < maxcalls) ) {
8656 l = TMath::Abs(LanGauFun(&x,params) - fy);
8670 // Search for left x location of fy
8672 p = maxx - 0.5 * params[0];
8678 while ((l != lold) && (i < maxcalls)) {
8682 l = TMath::Abs(LanGauFun(&x,params) - fy);
8684 step = -step / 10.0;
8689 if (i == maxcalls) {