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