AliRoot Core  3dc7879 (3dc7879)
CreateGainMap.C
Go to the documentation of this file.
1 
26 AliTPCCalPad * 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
31 
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 
185 
186 
187 void MakeCalibTree(char * inputKr="calibKr.root", char * inputCE ="fitCE.root", char * inputPulser=0){
188 
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 
220 TTree * tree =0;
221 
222 void LoadViewer(){
224 
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 
237 void 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 }
AliTPCCalROC * GetCalROC(Int_t sector) const
Definition: AliTPCCalPad.h:39
AliTPCCalPad * CreateGainMap(AliTPCCalPad *krypFitMean, AliTPCCalPad *krypFitRMS, AliTPCCalPad *noiseMap=0, AliTPCCalPad *krypSpectrMean=0, AliTPCCalPad *krypChi2=0, AliTPCCalPad *pulser=0, AliTPCCalPad *electrode=0)
Definition: CreateGainMap.C:26
TTree * tree
AliTPCCalibViewerGUI * viewer
Double_t GetMedian(AliTPCCalROC *const outlierROC=0, EPadType padType=kAll) const
#define TObjArray
GUI for the AliTPCCalibViewer used for the calibration monitor All functionalities of the AliTPCCalib...
void LoadViewer()
Float_t GetValue(UInt_t row, UInt_t pad) const
Definition: AliTPCCalROC.h:38
void Multiply(Float_t c1)
void Fit()
void GlobalFit(const AliTPCCalROC *ROCoutliers, Bool_t robust, TVectorD &fitParam, TMatrixD &covMatrix, Float_t &chi2, Int_t fitType=1, Double_t chi2Threshold=5, Double_t robustFraction=0.7, Double_t err=1, EPadType padType=kAll)
AliTPCCalibViewer * GetViewer()
TObjArray * array
Definition: AnalyzeLaser.C:12
AliTPCCalibViewer * makePad
void MakeCalibTree(char *inputKr="calibKr.root", char *inputCE="fitCE.root", char *inputPulser=0)
Double_t chi2
Definition: AnalyzeLaser.C:7
void test3()
void SetValue(UInt_t row, UInt_t pad, Float_t vd)
Definition: AliTPCCalROC.h:40
TPC calibration base class for one ROC.
Definition: AliTPCCalROC.h:20
TH2F * MakeHisto2D(Int_t side=0)
static TObjArray * ShowGUI(const char *fileName=0)
TF1 * f
Definition: interpolTest.C:21
UInt_t GetNchannels() const
Definition: AliTPCCalROC.h:34
void test()
Definition: interpolTest.C:100
void DumpToFile(const char *fileName)
Preprocessor class for HLT and DAQ.
static AliTPCCalROC * CreateGlobalFitCalROC(TVectorD &fitParam, Int_t sector, EPadType padType=kAll, AliTPCCalROC *oldTPCCalROC=0)
void test2()
TTree * GetTree() const
void test4()