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:
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
11 For more details see below.
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)
25 AliTPCCalPad * final = CreateGainMap(kryptonRaw, kryptonRMS)
27 TFile *h = new TFile("GainMap.root", "RECREATE")
33 AliTPCCalPad * CreateGainMap(AliTPCCalPad *krypFitMean, AliTPCCalPad *krypFitRMS, AliTPCCalPad *noiseMap = 0, AliTPCCalPad *krypSpectrMean = 0,
34 AliTPCCalPad *krypChi2 = 0, AliTPCCalPad *pulser = 0, AliTPCCalPad *electrode = 0) {
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);
44 const Double_t kryptonMean = 0.05012;
45 const Double_t kryptonSigma = 0.00386;
47 TObjArray arrayRocFinal(72);
50 // Loop over all sectors
52 for(Int_t isector=0; isector < 72; isector++){
54 AliTPCCalROC *rocKrypFitMean = krypFitMean->GetCalROC(isector);
55 AliTPCCalROC *rocKrypFitRMS = krypFitRMS->GetCalROC(isector);
57 AliTPCCalROC *rocOutlier = new AliTPCCalROC(*rocKrypFitMean);
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.
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);
75 //TCanvas *DEBUG = new TCanvas();
76 //rocOutlier->MakeHisto2D()->Draw("colz");
78 Double_t rocMedian = rocKrypFitMean->GetMedian(rocOutlier);
81 // 2. make a parabolic fit for the whole chamber with excluded outliers
86 rocKrypFitMean->GlobalFit(rocOutlier,0,params,cov,chi2);
87 AliTPCCalROC *rocParabolicFit=AliTPCCalROC::CreateGlobalFitCalROC(params,isector);
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);
97 // VERY IMPORTANT **** VERY IMPORTANT **** VERY IMPORTANT **** VERY IMPORTANT **** VERY IMPORTANT **** VERY IMPORTANT
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) !!!!!!!!!!
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);
110 // 3. replace outliers with fitted values
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));
120 if (TMath::Abs(fitRMS/fitMean - kryptonMean) < 4*kryptonSigma) {
121 rocFinal->SetValue(iChannel,rocKrypFitMean->GetValue(iChannel));
126 // 4. Postprocessing: Set dead channels and very noisy channels (time dependent) to 0
128 const Double_t noiseMin = 0.01;
129 const Double_t noiseMax = 2;
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);
138 // Fill an array of ROCs
140 arrayRocFinal.AddAt(rocFinal,isector);
144 AliTPCCalPad *final = new AliTPCCalPad(&arrayRocFinal);
146 // 4. normalize separately IROCs and OROCs to the mean of their medians
148 Double_t meanMedianIROC;
149 Double_t meanMedianOROC;
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();
159 meanMedianIROC = meanMedianIROC/n;
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();
169 meanMedianOROC = meanMedianOROC/n;
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);
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();
186 final->SetName("GainMap");
193 void MakeCalibTree(char * inputKr="calibKr.root", char * inputCE ="fitCE.root", char * inputPulser=0){
197 AliTPCPreprocessorOnline * preprocesor = new AliTPCPreprocessorOnline;
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");
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());
219 preprocesor->AddComponent(ceqIn->Clone());
220 preprocesor->AddComponent(ceqF1->Clone());
221 preprocesor->AddComponent(ceqF2->Clone());
223 preprocesor->DumpToFile("gainTree.root");
227 AliTPCCalibViewerGUI*viewer =0;
234 TObjArray * array = AliTPCCalibViewerGUI::ShowGUI("gainTree.root");
235 AliTPCCalibViewerGUI* viewer = (AliTPCCalibViewerGUI*)array->At(0);
236 makePad = viewer->GetViewer();
237 tree = viewer->GetViewer()->GetTree();
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");
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");