Update timestamps for new AMANDA simulation (17/02/2015)
[u/mrichter/AliRoot.git] / TPC / CalibMacros / CreateGainMap.C
CommitLineData
35d81915 1/// \file CreateGainMap.C
2/// This macro creates a gain map for the TPC based on the results of the Krypton calibration.
3/// The main steps are the following:
4///
5/// 1. Define outlier-pads where the krypton calibration was not succesful
6/// 2. A parabolic fit for the whole chamber is performed
7/// 3. replace outliers with fitted values
8/// 4. normalize separately IROCs and OROCs
9///
10/// For more details see below.
11///
12/// ~~~{.cpp}
13/// TFile f("calibKr.root")
14/// AliTPCCalPad * kryptonRaw = new AliTPCCalPad(*fitMean)
15/// AliTPCCalPad * kryptonMean = new AliTPCCalPad(*spectrMean)
16/// AliTPCCalPad * kryptonChi2 = new AliTPCCalPad(*fitNormChi2)
17/// AliTPCCalPad * kryptonRMS = new AliTPCCalPad(*fitRMS)
18///
19/// .L CreateGainMap.C
20/// AliTPCCalPad * final = CreateGainMap(kryptonRaw, kryptonRMS)
21///
22/// TFile *h = new TFile("GainMap.root", "RECREATE")
23/// final.Write()
24/// ~~~
6ba6d650 25
26AliTPCCalPad * CreateGainMap(AliTPCCalPad *krypFitMean, AliTPCCalPad *krypFitRMS, AliTPCCalPad *noiseMap = 0, AliTPCCalPad *krypSpectrMean = 0,
27 AliTPCCalPad *krypChi2 = 0, AliTPCCalPad *pulser = 0, AliTPCCalPad *electrode = 0) {
28
29
30 // Draw input map
35d81915 31
6ba6d650 32 TCanvas *test3 = new TCanvas("ASIDE3", "Aoriginal");
33 krypFitMean->MakeHisto2D()->Draw("colz");
34 TCanvas *test4 = new TCanvas("CSIDE4", "Coriginal");
35 krypFitMean->MakeHisto2D(1)->Draw("colz");
36 TH1F * hDEBUG = new TH1F("bla", "fitRMS/fitMean",100,0,0.3);
37
38 const Double_t kryptonMean = 0.05012;
39 const Double_t kryptonSigma = 0.00386;
40
41 TObjArray arrayRocFinal(72);
42
43 //
44 // Loop over all sectors
45 //
46 for(Int_t isector=0; isector < 72; isector++){
47
48 AliTPCCalROC *rocKrypFitMean = krypFitMean->GetCalROC(isector);
49 AliTPCCalROC *rocKrypFitRMS = krypFitRMS->GetCalROC(isector);
50
51 AliTPCCalROC *rocOutlier = new AliTPCCalROC(*rocKrypFitMean);
52
53 // 1. define map of outliers: the 41keV peak is fitted and krypFitMean represents the peak position and krypFitRMS its sigma.
54 // In order to control the fit qualtiy krypFitRMS/krypFitMean is computed and plotted; it is roughly gaussian distributed
55 // with the mean resolution kryptonMean = 0.05012 and sigma kryptonSigma = 0.00386. If the deviation is larger than 4sigma
56 // the fit to the 41keV peak is considered as failed.
57
58 Int_t nPointsFit = 0;
59 for(Int_t iChannel=0; iChannel < rocOutlier->GetNchannels(); iChannel++){
60 Double_t fitMean = rocKrypFitMean->GetValue(iChannel);
61 Double_t fitRMS = rocKrypFitRMS->GetValue(iChannel);
62 if (fitRMS < 0.001 || fitMean < 0.001) continue;
63 hDEBUG->Fill(fitRMS/fitMean);
64 if (TMath::Abs(fitRMS/fitMean - kryptonMean) < 4*kryptonSigma || fitRMS < 0.001 || fitMean < 0.001) {
65 rocOutlier->SetValue(iChannel,0);
66 nPointsFit++;
67 }
68 }
69 //TCanvas *DEBUG = new TCanvas();
70 //rocOutlier->MakeHisto2D()->Draw("colz");
71
72 Double_t rocMedian = rocKrypFitMean->GetMedian(rocOutlier);
73
74
75 // 2. make a parabolic fit for the whole chamber with excluded outliers
76
77 TVectorD params;
78 TMatrixD cov;
79 Float_t chi2;
80 rocKrypFitMean->GlobalFit(rocOutlier,0,params,cov,chi2);
81 AliTPCCalROC *rocParabolicFit=AliTPCCalROC::CreateGlobalFitCalROC(params,isector);
82
83 if (nPointsFit != 0) cout << "sector: "<< isector << " chi2: " << chi2/nPointsFit <<" median: "<< rocMedian <<" n: " << nPointsFit << endl;
84 if (nPointsFit == 0) {
85 AliTPCCalROC *rocFinal = new AliTPCCalROC(*rocParabolicFit); // What to do with sectors being switched off??
86 arrayRocFinal.AddAt(rocFinal,isector);
87 continue;
88 }
89
90 //
91 // VERY IMPORTANT **** VERY IMPORTANT **** VERY IMPORTANT **** VERY IMPORTANT **** VERY IMPORTANT **** VERY IMPORTANT
92 //
93 // if the fit is considered as failed, the mean value is taken for the whole chamber (TO AVOID THIS, RELEASE THIS CUT)
94 // if cosmic data can be used for gain calibration, this cut should be strong (e.g. 5) !!!!!!!!!!
95
96 if (chi2/nPointsFit > 5) {
97 AliTPCCalROC *rocFinal = new AliTPCCalROC(*rocParabolicFit);
98 if (rocMedian == 0 && isector == 51) rocMedian = 1075; // manual patch to be removed for sector 51!
99 for(Int_t iChannel=0; iChannel < rocFinal->GetNchannels(); iChannel++) rocFinal->SetValue(iChannel,rocMedian);
100 arrayRocFinal.AddAt(rocFinal,isector);
101 continue;
102 }
103
104 // 3. replace outliers with fitted values
105
106 AliTPCCalROC *rocFinal = new AliTPCCalROC(*rocParabolicFit);
107 for(Int_t iChannel=0; iChannel < rocFinal->GetNchannels(); iChannel++){
108 Double_t fitMean = rocKrypFitMean->GetValue(iChannel);
109 Double_t fitRMS = rocKrypFitRMS->GetValue(iChannel);
110 if (fitRMS < 0.001 || fitMean < 0.001 || TMath::Abs(rocParabolicFit->GetValue(iChannel)/fitMean - 1) > 0.35) {
111 rocFinal->SetValue(iChannel,rocParabolicFit->GetValue(iChannel));
112 continue;
113 }
114 if (TMath::Abs(fitRMS/fitMean - kryptonMean) < 4*kryptonSigma) {
115 rocFinal->SetValue(iChannel,rocKrypFitMean->GetValue(iChannel));
116 }
117 }
118
119
120 // 4. Postprocessing: Set dead channels and very noisy channels (time dependent) to 0
121
122 const Double_t noiseMin = 0.01;
123 const Double_t noiseMax = 2;
124
125 if (noiseMap) {
126 AliTPCCalROC *rocNoise = noiseMap->GetCalROC(isector);
127 for(Int_t iChannel=0; iChannel < rocFinal->GetNchannels(); iChannel++){
128 Double_t noise = rocNoise->GetValue(iChannel);
129 if (noise < noiseMin || noise > noiseMax) rocFinal->SetValue(iChannel, 0);
130 }
131 }
132 // Fill an array of ROCs
133
134 arrayRocFinal.AddAt(rocFinal,isector);
135
136 }
137
138 AliTPCCalPad *final = new AliTPCCalPad(&arrayRocFinal);
139
140 // 4. normalize separately IROCs and OROCs to the mean of their medians
141
142 Double_t meanMedianIROC;
143 Double_t meanMedianOROC;
144 Int_t n = 0;
145
146 for(Int_t isector=0; isector < 36; isector++){ // IROCs
147 AliTPCCalROC *rocFinal = final->GetCalROC(isector);
148 if (rocFinal->GetMedian() != 0) {
149 meanMedianIROC += rocFinal->GetMedian();
150 n++;
151 }
152 }
153 meanMedianIROC = meanMedianIROC/n;
154
155 n = 0;
156 for(Int_t isector=36; isector < 72; isector++){ // OROCs
157 AliTPCCalROC *rocFinal = final->GetCalROC(isector);
158 if (rocFinal->GetMedian() != 0) {
159 meanMedianOROC += rocFinal->GetMedian();
160 n++;
161 }
162 }
163 meanMedianOROC = meanMedianOROC/n;
164
165 for(Int_t isector=0; isector < 72; isector++){ // OROCs
166 AliTPCCalROC *rocFinal = final->GetCalROC(isector);
167 if (isector<36) rocFinal->Multiply(1./meanMedianIROC);
168 if (isector>35) rocFinal->Multiply(1./meanMedianOROC);
169 }
170
171 // Draw results
172 TCanvas *test = new TCanvas("ASIDE", "A");
173 final->MakeHisto2D()->Draw("colz");
174 TCanvas *test2 = new TCanvas("CSIDE", "C");
175 final->MakeHisto2D(1)->Draw("colz");
176 TCanvas *cDEBUG = new TCanvas();
177 hDEBUG->Draw();
178
179 //return results
180 final->SetName("GainMap");
181 return final;
182
183}
184
1b0bf200 185
186
187void MakeCalibTree(char * inputKr="calibKr.root", char * inputCE ="fitCE.root", char * inputPulser=0){
35d81915 188
1b0bf200 189 AliTPCPreprocessorOnline * preprocesor = new AliTPCPreprocessorOnline;
190 TFile f(inputKr);
191 TFile fce(inputCE);
192 AliTPCCalPad * kryptonMean = (AliTPCCalPad*)f.Get("spectrMean");
193 AliTPCCalPad * kryptonRMS = (AliTPCCalPad*)f.Get("spectrRMS");
194 AliTPCCalPad * kryptonMeanG = (AliTPCCalPad*)f.Get("fitMean");
195 AliTPCCalPad * kryptonRMSG = (AliTPCCalPad*)f.Get("fitRMS");
196 AliTPCCalPad * kryptonNormChi2 = (AliTPCCalPad*)f.Get("fitNormChi2");
197 AliTPCCalPad * kryptonEntries = (AliTPCCalPad*)f.Get("entries");
198 AliTPCCalPad * ceqIn = (AliTPCCalPad*)fce.Get("qIn");
199 AliTPCCalPad * ceqF1 = (AliTPCCalPad*)fce.Get("qF1");
200 AliTPCCalPad * ceqF2 = (AliTPCCalPad*)fce.Get("qF2");
201
202
203
204 preprocesor->AddComponent(kryptonMean->Clone());
205 preprocesor->AddComponent(kryptonRMS->Clone());
206 preprocesor->AddComponent(kryptonMeanG->Clone());
207 preprocesor->AddComponent(kryptonRMSG->Clone());
208 preprocesor->AddComponent(kryptonNormChi2->Clone());
209 preprocesor->AddComponent(kryptonEntries->Clone());
210 //
211 preprocesor->AddComponent(ceqIn->Clone());
212 preprocesor->AddComponent(ceqF1->Clone());
213 preprocesor->AddComponent(ceqF2->Clone());
214
215 preprocesor->DumpToFile("gainTree.root");
216 //
217}
218
219AliTPCCalibViewerGUI*viewer =0;
220TTree * tree =0;
221
222void LoadViewer(){
35d81915 223 /// Load calib Viewer
224
1b0bf200 225 TObjArray * array = AliTPCCalibViewerGUI::ShowGUI("gainTree.root");
226 AliTPCCalibViewerGUI* viewer = (AliTPCCalibViewerGUI*)array->At(0);
227 makePad = viewer->GetViewer();
228 tree = viewer->GetViewer()->GetTree();
229
230 tree->SetAlias("krAccept0","abs(fitRMS.fElements/fitMean.fElements-0.06)<0.04");
231 tree->SetAlias("krAccept1","abs(fitRMS.fElements)>30");
232 tree->SetAlias("yedge","tan(10*pi/180.)*lx.fElements");
233 tree->SetAlias("ceAccept1","abs(qIn.fElements/qIn_Median.fElements-1.5)<1.4&&qIn.fElements>3&&qIn_Median.fElements>3");
234
235}
236
237void Fit(){
238 TF1 f1("f1","[0]*exp(-[1]*x)+[2]");
239 f1.SetParameters(1,1,0.2);
240 tree->Draw("1-qIn.fElements/qF1.fElements:yedge-abs(ly.fElements)>>his(50,1.5,5,100,-0.5,1.5)","ceAccept1&&sector>36");
241 his->FitSlicesY();
242 his_1->Fit(&f1);
243
244
245}