]>
Commit | Line | Data |
---|---|---|
5253c153 | 1 | /************************************************************************** |
2 | * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. * | |
3 | * | |
4 | * Author: The ALICE Off-line Project. * | |
5 | * Contributors are mentioned in the code where appropriate. * | |
6 | * * | |
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 | **************************************************************************/ | |
15 | ||
16 | // $Id$ | |
17 | ||
18 | #include "AliMUONGain.h" | |
19 | #include "AliMUONPedestal.h" | |
20 | #include "AliMUONErrorCounter.h" | |
21 | #include "AliMUON2DMap.h" | |
22 | #include "AliMUONVCalibParam.h" | |
23 | #include "AliMUONVStore.h" | |
24 | #include "AliMpIntPair.h" | |
25 | ||
26 | #include <TString.h> | |
27 | #include <THashList.h> | |
28 | #include <TTimeStamp.h> | |
29 | #include <TTree.h> | |
30 | #include <TFile.h> | |
31 | #include <TF1.h> | |
32 | #include <TGraphErrors.h> | |
33 | #include <TMath.h> | |
34 | #include <Riostream.h> | |
35 | ||
36 | #include <sstream> | |
b785c4bd | 37 | #include <cstdio> |
5253c153 | 38 | |
39 | #define NFITPARAMS 4 | |
40 | ||
41 | //----------------------------------------------------------------------------- | |
42 | /// \class AliMUONGain | |
43 | /// | |
44 | /// Implementation of calibration computing | |
45 | /// | |
46 | /// add | |
47 | /// | |
48 | /// | |
49 | /// \author Alberto Baldisseri, JL Charvet | |
50 | //----------------------------------------------------------------------------- | |
51 | ||
52 | ||
53 | ||
54 | ||
55 | // functions | |
56 | ||
f3cfa63e | 57 | namespace { |
58 | ||
59 | //______________________________________________________________________________ | |
60 | Double_t funcLin (const Double_t *x, const Double_t *par) | |
61 | { | |
62 | /// Linear function | |
63 | return par[0] + par[1]*x[0]; | |
64 | } | |
65 | ||
66 | //______________________________________________________________________________ | |
67 | Double_t funcParabolic (const Double_t *x, const Double_t *par) | |
68 | { | |
69 | /// Parabolic function | |
70 | return par[0]*x[0]*x[0]; | |
71 | } | |
72 | ||
73 | //______________________________________________________________________________ | |
74 | Double_t funcCalib (const Double_t *x, const Double_t *par) | |
75 | { | |
76 | /// Calibration function | |
77 | Double_t xLim= par[3]; | |
78 | ||
79 | if(x[0] <= xLim) return par[0] + par[1]*x[0]; | |
80 | ||
81 | Double_t yLim = par[0]+ par[1]*xLim; | |
82 | return yLim + par[1]*(x[0] - xLim) + par[2]*(x[0] - xLim)*(x[0] - xLim); | |
83 | } | |
84 | ||
5253c153 | 85 | } |
86 | ||
b80faac0 | 87 | using std::endl; |
88 | using std::cout; | |
89 | using std::istringstream; | |
90 | using std::ostringstream; | |
5253c153 | 91 | /// \cond CLASSIMP |
92 | ClassImp(AliMUONGain) | |
93 | /// \endcond | |
94 | ||
95 | //______________________________________________________________________________ | |
96 | AliMUONGain::AliMUONGain() | |
97 | : AliMUONPedestal(), | |
98 | fInjCharge(0), | |
f3cfa63e | 99 | fRootDataFileName(), |
5253c153 | 100 | fnInit(1), |
101 | fnEntries(11), | |
102 | fnbpf1(6), | |
103 | fPrintLevel(0), | |
104 | fPlotLevel(0) | |
105 | { | |
106 | /// Default constructor | |
107 | ||
5253c153 | 108 | } |
5253c153 | 109 | |
110 | //______________________________________________________________________________ | |
111 | AliMUONGain::~AliMUONGain() | |
112 | { | |
113 | /// Destructor | |
114 | } | |
115 | ||
116 | //______________________________________________________________________________ | |
117 | TString AliMUONGain::WriteDummyHeader(void) | |
118 | { | |
119 | /// | |
120 | ||
121 | ostringstream stream; | |
122 | stream<<"//DUMMY FILE (to prevent Shuttle failure)"<< endl; | |
123 | stream<<"//================================================" << endl; | |
f3cfa63e | 124 | stream<<"// MUONTRKGAINda: Calibration run " << endl; |
5253c153 | 125 | stream<<"//================================================" << endl; |
126 | stream<<"// * Run : " << fRunNumber << endl; | |
127 | stream<<"// * Date : " << fDate->AsString("l") <<endl; | |
128 | stream<<"// * DAC : " << fInjCharge << endl; | |
129 | stream<<"//-------------------------------------------------" << endl; | |
130 | ||
131 | return TString(stream.str().c_str()); | |
132 | } | |
133 | ||
134 | //______________________________________________________________________________ | |
135 | void AliMUONGain::MakePedStoreForGain(TString shuttleFile) | |
136 | { | |
f3cfa63e | 137 | /// Store Pedmean and sigma to pedestal-like ascii file |
5253c153 | 138 | |
139 | ofstream fileout; | |
140 | TString tempstring; | |
f3cfa63e | 141 | TString flatFile; |
5253c153 | 142 | TString outputFile; |
143 | ||
144 | // Store pedestal map in root file | |
145 | TTree* tree = 0x0; | |
146 | ||
147 | // write dummy ascii file -> Shuttle | |
148 | if(fIndex<fnEntries) | |
149 | { | |
5253c153 | 150 | FILE *pfilew=0; |
151 | pfilew = fopen (shuttleFile,"w"); | |
152 | ||
153 | fprintf(pfilew,"//DUMMY FILE (to prevent Shuttle failure)\n"); | |
154 | fprintf(pfilew,"//================================================\n"); | |
f3cfa63e | 155 | fprintf(pfilew,"// MUONTRKGAINda: Calibration run \n"); |
5253c153 | 156 | fprintf(pfilew,"//=================================================\n"); |
157 | fprintf(pfilew,"// * Run : %d \n",fRunNumber); | |
158 | fprintf(pfilew,"// * Date : %s \n",fDate->AsString("l")); | |
159 | fprintf(pfilew,"// * DAC : %d \n",fInjCharge); | |
160 | fprintf(pfilew,"//-------------------------------------------------\n"); | |
161 | fclose(pfilew); | |
162 | } | |
163 | ||
164 | ||
165 | ||
dbbb2c64 | 166 | Finalize(); |
167 | MakeControlHistos(); | |
5253c153 | 168 | if(fPrintLevel>0) |
169 | { | |
170 | // compute and store mean DAC values (like pedestals) | |
f3cfa63e | 171 | flatFile = Form("%s.ped",fPrefixDA.Data()); |
5253c153 | 172 | outputFile=flatFile; |
9c2d4e89 | 173 | cout << "\n" << fPrefixLDC.Data() << " : Flat file generated : " << flatFile.Data() << "\n"; |
dbbb2c64 | 174 | if (!outputFile.IsNull()) |
175 | { | |
176 | ofstream out(outputFile.Data()); | |
177 | MakeASCIIoutput(out); | |
178 | out.close(); | |
179 | } | |
5253c153 | 180 | } |
5253c153 | 181 | |
182 | TString mode("UPDATE"); | |
183 | ||
184 | if (fIndex==1) { | |
185 | mode = "RECREATE"; | |
186 | } | |
f3cfa63e | 187 | TFile* histoFile = new TFile(fRootDataFileName.Data(), mode.Data(), "MUON Tracking Gains"); |
5253c153 | 188 | |
189 | // second argument should be the injected charge, taken from config crocus file | |
190 | // put also info about run number could be usefull | |
191 | AliMpIntPair* pair = new AliMpIntPair(fRunNumber,fInjCharge ); | |
192 | ||
193 | if (mode.CompareTo("UPDATE") == 0) { | |
194 | tree = (TTree*)histoFile->Get("t"); | |
195 | tree->SetBranchAddress("run",&pair); | |
196 | tree->SetBranchAddress("ped",&fPedestalStore); | |
197 | ||
198 | } else { | |
199 | tree = new TTree("t","Pedestal tree"); | |
200 | tree->Branch("run", "AliMpIntPair",&pair); | |
201 | tree->Branch("ped", "AliMUON2DMap",&fPedestalStore); | |
202 | tree->SetBranchAddress("run",&pair); | |
203 | tree->SetBranchAddress("ped",&fPedestalStore); | |
204 | ||
205 | } | |
206 | ||
207 | tree->Fill(); | |
208 | tree->Write("t", TObject::kOverwrite); // overwrite the tree | |
209 | histoFile->Close(); | |
210 | ||
211 | delete pair; | |
212 | } | |
213 | ||
214 | //______________________________________________________________________________ | |
5431405e | 215 | TString AliMUONGain::WriteGainHeader(Int_t nInit, Int_t nEntries, Int_t nbpf2, Int_t *numrun, Double_t *injCharge) |
5253c153 | 216 | { |
1b15b395 | 217 | /// Header of the calibration output file |
5253c153 | 218 | |
219 | ostringstream stream; | |
220 | ||
221 | ||
1b15b395 | 222 | stream<<"//=======================================================" << endl; |
f3cfa63e | 223 | stream<<"// Calibration file calculated by " << fPrefixDA.Data() <<endl; |
1b15b395 | 224 | stream<<"//=======================================================" << endl; |
5253c153 | 225 | stream<<"// * Run : " << fRunNumber << endl; |
226 | stream<<"// * Date : " << fDate->AsString("l") <<endl; | |
227 | stream<<"// * Statictics : " << fNEvents << endl; | |
1b15b395 | 228 | if(fConfig) |
f3cfa63e | 229 | stream<<"// * # of MANUS : " << fNManuConfig << " read in the Det. config. " << endl; |
1b15b395 | 230 | stream<<"// * # of MANUS : " << fNManu << " read in raw data " << endl; |
231 | stream<<"// * # of MANUS : " << fNChannel/64 << " written in calibration file " << endl; | |
5253c153 | 232 | stream<<"// * # of channels : " << fNChannel << endl; |
1b15b395 | 233 | stream<<"//-------------------------------------------------------" << endl; |
5253c153 | 234 | |
5431405e | 235 | if(nInit==0) |
236 | stream<<"// "<< nEntries <<" DAC values fit: "<< fnbpf1 << " pts (1st order) " << nbpf2 << " pts (2nd order)" << endl; | |
237 | if(nInit==1) | |
238 | stream<<"// "<< nEntries <<" DAC values fit: "<< fnbpf1 << " pts (1st order) " << nbpf2 << " pts (2nd order) DAC=0 excluded" << endl; | |
f3cfa63e | 239 | stream<<"// * nInit = " << nInit << " * f1nbp = " << fnbpf1 << " * f2nbp = " << nbpf2 << endl; |
5253c153 | 240 | |
241 | stream<<"// RUN DAC " << endl; | |
242 | stream<<"//-----------------" << endl; | |
5431405e | 243 | for (Int_t i = 0; i < nEntries; ++i) { |
244 | stream<<Form("// %d %5.0f \n",numrun[i],injCharge[i]); | |
5253c153 | 245 | } |
246 | stream<<"//=======================================" << endl; | |
247 | stream<<"// BP MANU CH. a1 a2 thres. q " << endl; | |
248 | stream<<"//=======================================" << endl; | |
249 | ||
250 | return TString(stream.str().c_str()); | |
251 | } | |
252 | ||
253 | //______________________________________________________________________________ | |
254 | TString AliMUONGain::WriteGainData(Int_t BP, Int_t Manu, Int_t ch, Double_t p1, Double_t p2, Int_t threshold, Int_t q) | |
255 | { | |
1b15b395 | 256 | /// Write calibration parameters per channel |
5253c153 | 257 | |
258 | ostringstream stream(""); | |
5431405e | 259 | stream << Form("%4i %5i %2i %7.4f %10.3e %4i %2x\n",BP,Manu,ch,p1,p2,threshold,q); |
5253c153 | 260 | return TString(stream.str().c_str()); |
261 | ||
262 | } | |
263 | ||
da329aff | 264 | //_______________________________________________________________________________ |
5253c153 | 265 | void AliMUONGain::MakeGainStore(TString shuttleFile) |
266 | { | |
9c2d4e89 | 267 | Int_t status=0; |
5253c153 | 268 | /// Store gains in ASCII files |
269 | ofstream fileout; | |
270 | ofstream filcouc; | |
271 | TString tempstring; | |
f3cfa63e | 272 | TString filename; |
9c2d4e89 | 273 | char* detail; |
5253c153 | 274 | |
275 | Double_t goodA1Min = 0.5; | |
276 | Double_t goodA1Max = 2.; | |
f3cfa63e | 277 | // Double_t goodA2Min = -0.5E-03; |
278 | // Double_t goodA2Max = 1.E-03; | |
f9e83f1d | 279 | Double_t goodA2Min = -0.5E-01; // changed 28/10/2009 (JLC) <=> enlarged condition on a2 |
f3cfa63e | 280 | Double_t goodA2Max = 1.E-01; |
5253c153 | 281 | // Table for uncalibrated buspatches and manus |
282 | THashList* uncalBuspatchManuTable = new THashList(1000,2); | |
283 | ||
284 | Int_t numrun[11]; | |
285 | ||
286 | // open file MUONTRKda_gain_data.root | |
287 | // read again the pedestal for the calibration runs (11 runs) | |
288 | // need the injection charge from config file (to be done) | |
289 | // For each channel make a TGraphErrors (mean, sigma) vs injected charge | |
290 | // Fit with a polynomial fct | |
291 | // store the result in a flat file. | |
292 | ||
f3cfa63e | 293 | if(fIndex==0)cout << " Root data file = " << fRootDataFileName.Data() << endl; |
294 | TFile* histoFile = new TFile(fRootDataFileName.Data()); | |
5253c153 | 295 | |
296 | AliMUON2DMap* map[11]; | |
297 | AliMUONVCalibParam* ped[11]; | |
298 | AliMpIntPair* run[11]; | |
299 | ||
300 | //read back from root file | |
301 | TTree* tree = (TTree*)histoFile->Get("t"); | |
302 | Int_t nEntries = tree->GetEntries(); | |
303 | ||
304 | Int_t nbpf2 = nEntries - (fnInit + fnbpf1) + 1; // nb pts used for 2nd order fit | |
305 | ||
306 | // read back info | |
307 | for (Int_t i = 0; i < nEntries; ++i) { | |
308 | map[i] = 0x0; | |
309 | run[i] = 0x0; | |
310 | tree->SetBranchAddress("ped",&map[i]); | |
311 | tree->SetBranchAddress("run",&run[i]); | |
312 | tree->GetEvent(i); | |
313 | // std::cout << map[i] << " " << run[i] << std::endl; | |
314 | } | |
315 | // RunNumber extracted from Root data fil | |
316 | if(fIndex==0)fRunNumber=(UInt_t)run[nEntries-1]->GetFirst(); | |
317 | // sscanf(getenv("DATE_RUN_NUMBER"),"%d",&gAliRunNumber); | |
318 | ||
319 | Double_t pedMean[11]; | |
320 | Double_t pedSigma[11]; | |
321 | for ( Int_t i=0 ; i<11 ; i++) {pedMean[i]=0.;pedSigma[i]=1.;}; | |
322 | Double_t injCharge[11]; | |
323 | Double_t injChargeErr[11]; | |
324 | for ( Int_t i=0 ; i<11 ; i++) {injCharge[i]=0.;injChargeErr[i]=1.;}; | |
325 | ||
5253c153 | 326 | // print out in .log file |
9c2d4e89 | 327 | detail=Form("\n%s : ------ MUONTRKGAINda for Gain computing (Last Run = %d) ------\n",fPrefixLDC.Data(),fRunNumber); printf("%s",detail); |
5253c153 | 328 | (*fFilcout)<<"\n\n//=================================================" << endl; |
f3cfa63e | 329 | (*fFilcout)<<"// MUONTRKGAINda: Gain Computing Run = " << fRunNumber << endl; |
330 | (*fFilcout)<<"// RootDataFile = "<< fRootDataFileName.Data() << endl; | |
5253c153 | 331 | (*fFilcout)<<"//=================================================" << endl; |
332 | (*fFilcout)<<"//* Date : " << fDate->AsString("l") << "\n" << endl; | |
333 | ||
9c2d4e89 | 334 | (*fFilcout) << " Entries = " << nEntries << " DAC values \n" << endl; |
335 | for (Int_t i = 0; i < nEntries; ++i) { | |
336 | (*fFilcout) << " Run = " << run[i]->GetFirst() << " DAC = " << run[i]->GetSecond() << endl; | |
337 | numrun[i] = run[i]->GetFirst(); | |
338 | injCharge[i] = run[i]->GetSecond(); | |
339 | injChargeErr[i] = 0.01*injCharge[i]; | |
340 | if(injChargeErr[i] <= 1.) injChargeErr[i]=1.; | |
341 | } | |
342 | detail=Form("%s : .... Fitting .... \n",fPrefixLDC.Data()); printf("%s",detail); | |
5253c153 | 343 | |
344 | ||
345 | // why 2 files ? (Ch. F.) => second file contains detailed results | |
346 | FILE *pfilen = 0; | |
347 | if(fPrintLevel>1) | |
348 | { | |
f3cfa63e | 349 | filename=Form("%s.param",fPrefixDA.Data()); |
9c2d4e89 | 350 | detail=Form("%s : Second fit parameter file = %s\n",fPrefixLDC.Data(),filename.Data()); printf("%s",detail); |
351 | // cout << " Second fit parameter file = " << filename.Data() << "\n"; | |
f3cfa63e | 352 | pfilen = fopen (filename.Data(),"w"); |
5253c153 | 353 | |
354 | fprintf(pfilen,"//===================================================================\n"); | |
355 | fprintf(pfilen,"// BP MANU CH. par[0] [1] [2] [3] xlim P(chi2) p1 P(chi2)2 p2\n"); | |
356 | fprintf(pfilen,"//===================================================================\n"); | |
357 | fprintf(pfilen,"// * Run : %d \n",fRunNumber); | |
358 | fprintf(pfilen,"//===================================================================\n"); | |
359 | } | |
360 | ||
5253c153 | 361 | |
dbbb2c64 | 362 | // file outputs for gain |
5253c153 | 363 | |
5431405e | 364 | ofstream pfilew; |
365 | pfilew.open(shuttleFile.Data()); | |
dbbb2c64 | 366 | // Write Header Data of the .par file |
5431405e | 367 | pfilew << WriteGainHeader(fnInit,nEntries,nbpf2,numrun,injCharge); |
5253c153 | 368 | |
369 | // print mean and sigma values in file | |
370 | FILE *pfilep = 0; | |
371 | if(fPrintLevel>1) | |
372 | { | |
f3cfa63e | 373 | filename=Form("%s.peak",fPrefixDA.Data()); |
9c2d4e89 | 374 | detail=Form("%s : File containing Peak mean values = %s\n",fPrefixLDC.Data(),filename.Data()); printf("%s",detail); |
375 | // cout << " File containing Peak mean values = " << filename << "\n"; | |
5253c153 | 376 | pfilep = fopen (filename,"w"); |
377 | ||
378 | fprintf(pfilep,"//==============================================================================================================================\n"); | |
379 | fprintf(pfilep,"// * Run : %d \n",fRunNumber); | |
380 | fprintf(pfilep,"//==============================================================================================================================\n"); | |
381 | fprintf(pfilep,"// BP MANU CH. Ped. <0> <1> <2> <3> <4> <5> <6> <7> <8> <9> <10> \n"); | |
382 | fprintf(pfilep,"//==============================================================================================================================\n"); | |
383 | fprintf(pfilep,"// DAC= %9.1f%9.1f%9.1f%9.1f%9.1f%9.1f%9.1f%9.1f%9.1f%9.1f%9.1f fC\n",injCharge[0],injCharge[1],injCharge[2],injCharge[3],injCharge[4],injCharge[5],injCharge[6],injCharge[7],injCharge[8],injCharge[9],injCharge[10]); | |
384 | fprintf(pfilep,"// %9.1f%9.1f%9.1f%9.1f%9.1f%9.1f%9.1f%9.1f%9.1f%9.1f%9.1f\n",injChargeErr[0],injChargeErr[1],injChargeErr[2],injChargeErr[3],injChargeErr[4],injChargeErr[5],injChargeErr[6],injChargeErr[7],injChargeErr[8],injChargeErr[9],injChargeErr[10]); | |
385 | fprintf(pfilep,"//==============================================================================================================================\n"); | |
386 | } | |
387 | ||
388 | Double_t chi2 = 0.; | |
389 | Double_t chi2P2 = 0.; | |
390 | Double_t prChi2 = 0; | |
391 | Double_t prChi2P2 =0; | |
392 | Double_t a0=0.,a1=1.,a2=0.; | |
393 | Int_t busPatchId ; | |
394 | Int_t manuId ; | |
395 | Int_t channelId ; | |
396 | Int_t threshold = 0; | |
397 | Int_t q = 0; | |
398 | Int_t p1 =0; | |
399 | Int_t p2 =0; | |
f3cfa63e | 400 | Double_t gain=10.; // max value (= bad value) |
5253c153 | 401 | Double_t capa=0.2; // internal capacitor (pF) |
402 | ||
403 | // plot out | |
404 | ||
5253c153 | 405 | TTree* tg = 0x0; |
406 | if(fPlotLevel>0) | |
407 | { | |
f3cfa63e | 408 | fHistoFileName=Form("%s.root",fPrefixDA.Data()); |
065a2d7d | 409 | new TFile(fHistoFileName.Data(),"RECREATE","MUON Tracking gains"); |
5253c153 | 410 | tg = new TTree("tg","TTree avec class Manu_DiMu"); |
411 | ||
412 | tg->Branch("bp",&busPatchId, "busPatchId/I"); | |
413 | tg->Branch("manu",&manuId, "manuId/I"); | |
414 | tg->Branch("channel",&channelId, "channelId/I"); | |
415 | ||
416 | tg->Branch("a0",&a0, "a0/D"); | |
417 | tg->Branch("a1",&a1, "a1/D"); | |
418 | tg->Branch("a2",&a2, "a2/D"); | |
419 | tg->Branch("Pchi2",&prChi2, "prChi2/D"); | |
420 | tg->Branch("Pchi2_2",&prChi2P2, "prChi2P2/D"); | |
421 | tg->Branch("Threshold",&threshold, "threshold/I"); | |
422 | tg->Branch("q",&q, "q/I"); | |
423 | tg->Branch("p1",&p1, "p1/I"); | |
424 | tg->Branch("p2",&p2, "p2/I"); | |
425 | tg->Branch("gain",&gain, "gain/D"); | |
426 | } | |
427 | ||
428 | char graphName[256]; | |
429 | ||
430 | // iterates over the first pedestal run | |
431 | TIter next(map[0]->CreateIterator()); | |
432 | AliMUONVCalibParam* p; | |
433 | ||
434 | Int_t nmanu = 0; | |
435 | Int_t nGoodChannel = 0; | |
436 | Int_t nBadChannel = 0; | |
437 | Int_t noFitChannel = 0; | |
438 | Int_t nplot=0; | |
439 | Double_t sumProbChi2 = 0.; | |
440 | Double_t sumA1 = 0.; | |
441 | Double_t sumProbChi2P2 = 0.; | |
442 | Double_t sumA2 = 0.; | |
443 | ||
444 | Double_t x[11], xErr[11], y[11], yErr[11]; | |
445 | Double_t xp[11], xpErr[11], yp[11], ypErr[11]; | |
446 | ||
447 | Int_t uncalcountertotal=0 ; | |
1ccd531d | 448 | Int_t unparabolicfit=0; |
5253c153 | 449 | |
450 | while ( ( p = dynamic_cast<AliMUONVCalibParam*>(next() ) ) ) | |
451 | { | |
452 | ped[0] = p; | |
1ccd531d | 453 | unparabolicfit=0; |
5253c153 | 454 | |
455 | busPatchId = p->ID0(); | |
456 | manuId = p->ID1(); | |
457 | ||
458 | // read back pedestal from the other runs for the given (bupatch, manu) | |
459 | for (Int_t i = 1; i < nEntries; ++i) { | |
460 | ped[i] = static_cast<AliMUONVCalibParam*>(map[i]->FindObject(busPatchId, manuId)); | |
461 | } | |
462 | ||
463 | // compute for each channel the gain parameters | |
464 | for ( channelId = 0; channelId < ped[0]->Size() ; ++channelId ) | |
465 | { | |
466 | ||
467 | Int_t n = 0; | |
468 | for (Int_t i = 0; i < nEntries; ++i) { | |
469 | ||
470 | if (!ped[i]) continue; //shouldn't happen. | |
471 | pedMean[i] = ped[i]->ValueAsDouble(channelId, 0); | |
472 | pedSigma[i] = ped[i]->ValueAsDouble(channelId, 1); | |
473 | ||
474 | if (pedMean[i] < 0) continue; // not connected | |
475 | ||
476 | if (pedSigma[i] <= 0) pedSigma[i] = 1.; // should not happen. | |
477 | n++; | |
478 | } | |
479 | ||
5253c153 | 480 | // print_peak_mean_values |
481 | if(fPrintLevel>1) | |
482 | { | |
483 | ||
484 | fprintf(pfilep,"%4i%5i%5i%10.3f",busPatchId,manuId,channelId,0.); | |
485 | fprintf(pfilep,"%9.1f%9.1f%9.1f%9.1f%9.1f%9.1f%9.1f%9.1f%9.1f%9.1f%9.1f \n",pedMean[0],pedMean[1],pedMean[2],pedMean[3],pedMean[4],pedMean[5],pedMean[6],pedMean[7],pedMean[8],pedMean[9],pedMean[10]); | |
486 | fprintf(pfilep," sig= %9.3f%9.3f%9.3f%9.3f%9.3f%9.3f%9.3f%9.3f%9.3f%9.3f%9.3f \n",pedSigma[0],pedSigma[1],pedSigma[2],pedSigma[3],pedSigma[4],pedSigma[5],pedSigma[6],pedSigma[7],pedSigma[8],pedSigma[9],pedSigma[10]); | |
487 | } | |
488 | ||
9c2d4e89 | 489 | // makegain |
5253c153 | 490 | // Fit Method: Linear fit over gAlinbpf1 points + parabolic fit over nbpf2 points) |
491 | // nInit=1 : 1st pt DAC=0 excluded | |
492 | ||
493 | // 1. - linear fit over gAlinbpf1 points | |
494 | ||
2942f542 | 495 | Double_t par[4] = {0.,0.5,0.,static_cast<Double_t>(ADCMax())}; |
5253c153 | 496 | Int_t nbs = nEntries - fnInit; |
497 | if(nbs < fnbpf1)fnbpf1=nbs; | |
498 | ||
499 | Int_t fitproceed=1; | |
f3cfa63e | 500 | Int_t nbpf2Dynamic=nbpf2; |
9ee1d6ff | 501 | Int_t adcLimit=4090; // when RMS < 0.5 (in other cases mean values forced to 4095, see DA_PED) |
5253c153 | 502 | for (Int_t j = 0; j < nbs; ++j) |
503 | { | |
504 | Int_t k = j + fnInit; | |
505 | x[j] = pedMean[k]; | |
f3cfa63e | 506 | if(x[j]<=0.){fitproceed=0; break;} |
f9e83f1d | 507 | // if(x[j]>= ADCMax()) |
9ee1d6ff | 508 | if(x[j]>= adcLimit) |
f3cfa63e | 509 | { |
510 | if(j < nbs-1){fitproceed=0; break;} | |
f9e83f1d | 511 | else { nbpf2Dynamic=nbpf2-1; break;} |
f3cfa63e | 512 | } |
5253c153 | 513 | xErr[j] = pedSigma[k]; |
514 | y[j] = injCharge[k]; | |
515 | yErr[j] = injChargeErr[k]; | |
516 | ||
517 | } | |
518 | ||
519 | TGraphErrors *graphErr; | |
520 | if(!fitproceed) { p1=0; p2=0; noFitChannel++;} | |
521 | ||
522 | if(fitproceed) | |
523 | { | |
524 | ||
f3cfa63e | 525 | TF1 *f1 = new TF1("f1",funcLin,0.,ADCMax(),2); |
5253c153 | 526 | graphErr = new TGraphErrors(fnbpf1, x, y, xErr, yErr); |
527 | ||
528 | f1->SetParameters(0,0); | |
529 | ||
530 | graphErr->Fit("f1","RQ"); | |
531 | ||
532 | chi2 = f1->GetChisquare(); | |
533 | f1->GetParameters(par); | |
534 | ||
535 | delete graphErr; | |
536 | graphErr=0; | |
537 | delete f1; | |
538 | ||
539 | prChi2 = TMath::Prob(chi2, fnbpf1 - 2); | |
540 | ||
541 | Double_t xLim = pedMean[fnInit + fnbpf1 - 1]; | |
542 | Double_t yLim = par[0]+par[1] * xLim; | |
543 | ||
544 | a0 = par[0]; | |
545 | a1 = par[1]; | |
546 | ||
547 | // 2. - Translation : new origin (xLim, yLim) + parabolic fit over nbf2 points | |
1ccd531d | 548 | //checking: if(busPatchId ==1841 && manuId==4)nbpf2Dynamic=2; |
549 | if(nbpf2Dynamic > 2) | |
5253c153 | 550 | { |
f3cfa63e | 551 | for (Int_t j = 0; j < nbpf2Dynamic; j++) |
5253c153 | 552 | { |
553 | Int_t k = j + (fnInit + fnbpf1) - 1; | |
554 | xp[j] = pedMean[k] - xLim; | |
555 | xpErr[j] = pedSigma[k]; | |
556 | ||
557 | yp[j] = injCharge[k] - yLim - par[1]*xp[j]; | |
558 | ypErr[j] = injChargeErr[k]; | |
559 | } | |
560 | ||
f3cfa63e | 561 | TF1 *f2 = new TF1("f2",funcParabolic,0.,ADCMax(),1); |
562 | graphErr = new TGraphErrors(nbpf2Dynamic, xp, yp, xpErr, ypErr); | |
5253c153 | 563 | |
564 | graphErr->Fit(f2,"RQ"); | |
565 | chi2P2 = f2->GetChisquare(); | |
566 | f2->GetParameters(par); | |
567 | ||
568 | delete graphErr; | |
569 | graphErr=0; | |
570 | delete f2; | |
571 | ||
f3cfa63e | 572 | prChi2P2 = TMath::Prob(chi2P2, nbpf2Dynamic-1); |
5253c153 | 573 | a2 = par[0]; |
574 | } | |
1ccd531d | 575 | else |
576 | { | |
577 | unparabolicfit++; | |
578 | (*fFilcout) << " Warning : BP = " << busPatchId << " Manu = " << manuId << " Channel = " << channelId <<": parabolic fit not possible (nbpf2=" << nbpf2Dynamic << ") => a2=0 and linear fit OK" << std::endl; | |
579 | if(unparabolicfit==1) std::cout << " Warning : BP = " << busPatchId << " Manu = " << manuId << ": no parabolic fit for some channels (nbpf2=" << nbpf2Dynamic << "), linear fit is OK (see .log for details)" << std::endl; | |
580 | a2=0. ; prChi2P2=0. ; | |
581 | } | |
5253c153 | 582 | |
583 | par[0] = a0; | |
584 | par[1] = a1; | |
585 | par[2] = a2; | |
586 | par[3] = xLim; | |
587 | ||
5253c153 | 588 | if(prChi2>0.999999)prChi2=0.999999 ; if(prChi2P2>0.999999)prChi2P2=0.9999999; // avoiding Pr(Chi2)=1 value |
589 | p1 = TMath::Nint(floor(prChi2*15))+1; // round down value : floor(2.8)=2. | |
590 | p2 = TMath::Nint(floor(prChi2P2*15))+1; | |
591 | q = p1*16 + p2; // fit quality | |
592 | ||
593 | Double_t x0 = -par[0]/par[1]; // value of x corresponding to à 0 fC | |
594 | threshold = TMath::Nint(ceil(par[3]-x0)); // linear if x < threshold | |
595 | ||
596 | if(fPrintLevel>1) | |
597 | { | |
598 | fprintf(pfilen,"%4i %4i %2i",busPatchId,manuId,channelId); | |
599 | fprintf(pfilen," %6.2f %6.4f %10.3e %4.2f %4i %8.6f %8.6f %x %8.6f %8.6f %x\n", | |
600 | par[0], par[1], par[2], par[3], threshold, prChi2, floor(prChi2*15), p1, prChi2P2, floor(prChi2P2*15),p2); | |
601 | } | |
602 | ||
603 | ||
604 | // tests | |
605 | if(par[1]< goodA1Min || par[1]> goodA1Max) p1=0; | |
606 | if(par[2]< goodA2Min || par[2]> goodA2Max) p2=0; | |
607 | ||
608 | } // fitproceed | |
609 | ||
610 | if(fitproceed && p1>0 && p2>0) | |
611 | { | |
612 | nGoodChannel++; | |
613 | sumProbChi2 += prChi2; | |
614 | sumA1 += par[1]; | |
5253c153 | 615 | sumProbChi2P2 += prChi2P2; |
616 | sumA2 += par[2]; | |
617 | } | |
618 | else // bad calibration | |
619 | { | |
620 | nBadChannel++; | |
621 | q=0; | |
622 | par[1]=0.5; a1=0.5; p1=0; | |
623 | par[2]=0.; a2=0.; p2=0; | |
f3cfa63e | 624 | threshold=ADCMax(); |
5253c153 | 625 | |
626 | // bad calibration counter | |
627 | char bpmanuname[256]; | |
628 | AliMUONErrorCounter* uncalcounter; | |
629 | ||
b785c4bd | 630 | snprintf(bpmanuname,256,"bp%dmanu%d",busPatchId,manuId); |
5253c153 | 631 | if (!(uncalcounter = (AliMUONErrorCounter*)uncalBuspatchManuTable->FindObject(bpmanuname))) |
632 | { | |
633 | // New buspatch_manu name | |
634 | uncalcounter= new AliMUONErrorCounter (busPatchId,manuId); | |
635 | uncalcounter->SetName(bpmanuname); | |
636 | uncalBuspatchManuTable->Add(uncalcounter); | |
637 | } | |
638 | else | |
639 | { | |
640 | // Existing buspatch_manu name | |
641 | uncalcounter->Increment(); | |
642 | } | |
643 | // uncalcounter->Print_uncal() | |
644 | uncalcountertotal ++; | |
645 | } | |
5431405e | 646 | gain=1./(par[1]*capa); // mv/fC |
5253c153 | 647 | |
648 | if(fPlotLevel>0) | |
649 | {if(fPlotLevel>1) | |
650 | { | |
651 | // if(q==0 and nplot < 100) | |
652 | // if(p1>1 && p2==0 and nplot < 100) | |
f3cfa63e | 653 | // if(p1>10 && p2>10 and nplot < 100) |
5253c153 | 654 | // if(p1>=1 and p1<=2 and nplot < 100) |
5431405e | 655 | // if((p1==1 || p2==1) and nplot < 100) |
f3cfa63e | 656 | if(nbpf2Dynamic<nbpf2 and nplot < 100) |
5253c153 | 657 | { |
658 | nplot++; | |
659 | // cout << " nplot = " << nplot << endl; | |
f3cfa63e | 660 | TF1 *f2Calib = new TF1("f2Calib",funcCalib,0.,ADCMax(),NFITPARAMS); |
5253c153 | 661 | |
662 | graphErr = new TGraphErrors(nEntries,pedMean,injCharge,pedSigma,injChargeErr); | |
663 | ||
b785c4bd | 664 | snprintf(graphName,256,"BusPatch_%d_Manu_%d_Ch_%d",busPatchId, manuId,channelId); |
5253c153 | 665 | |
666 | graphErr->SetTitle(graphName); | |
667 | graphErr->SetMarkerColor(3); | |
668 | graphErr->SetMarkerStyle(12); | |
669 | graphErr->Write(graphName); | |
670 | ||
b785c4bd | 671 | snprintf(graphName,256,"f2_BusPatch_%d_Manu_%d_Ch_%d",busPatchId, manuId,channelId); |
5253c153 | 672 | f2Calib->SetTitle(graphName); |
673 | f2Calib->SetLineColor(4); | |
674 | f2Calib->SetParameters(par); | |
675 | f2Calib->Write(graphName); | |
676 | ||
677 | delete graphErr; | |
678 | graphErr=0; | |
679 | delete f2Calib; | |
680 | } | |
681 | } | |
5253c153 | 682 | tg->Fill(); |
683 | } | |
dbbb2c64 | 684 | pfilew << WriteGainData(busPatchId,manuId,channelId,par[1],par[2],threshold,q); |
5253c153 | 685 | } |
686 | nmanu++; | |
1ccd531d | 687 | Int_t step=500; |
9c2d4e89 | 688 | if(nmanu % step == 0)printf("%s : Nb manu = %d\n",fPrefixLDC.Data(),nmanu); |
5253c153 | 689 | } |
690 | ||
691 | // print in logfile | |
692 | if (uncalBuspatchManuTable->GetSize()) | |
693 | { | |
694 | uncalBuspatchManuTable->Sort(); // use compare | |
695 | TIterator* iter = uncalBuspatchManuTable->MakeIterator(); | |
696 | AliMUONErrorCounter* uncalcounter; | |
697 | (*fFilcout) << "\n List of problematic BusPatch and Manu " << endl; | |
698 | (*fFilcout) << " ========================================" << endl; | |
699 | (*fFilcout) << " BP Manu Nb Channel " << endl ; | |
700 | (*fFilcout) << " ========================================" << endl; | |
701 | while((uncalcounter = (AliMUONErrorCounter*) iter->Next())) | |
702 | { | |
703 | (*fFilcout) << "\t" << uncalcounter->BusPatch() << "\t " << uncalcounter->ManuId() << "\t\t" << uncalcounter->Events() << endl; | |
704 | } | |
705 | (*fFilcout) << " ========================================" << endl; | |
706 | ||
707 | (*fFilcout) << " Number of bad calibrated Manu = " << uncalBuspatchManuTable->GetSize() << endl ; | |
708 | (*fFilcout) << " Number of bad calibrated channel = " << uncalcountertotal << endl; | |
709 | ||
710 | } | |
9c2d4e89 | 711 | if(nmanu && nGoodChannel) |
0e5fb403 | 712 | { |
9c2d4e89 | 713 | Double_t ratio_limit=0.25; |
714 | Double_t ratio=float (nBadChannel)/float (nmanu*64); | |
715 | ||
716 | detail=Form("\n%s : Nb of channels in raw data = %d (%d Manu)",fPrefixLDC.Data(),nmanu*64,nmanu); | |
717 | (*fFilcout) << detail ; printf("%s",detail); | |
718 | detail=Form("\n%s : Nb of calibrated channel = %d (%4.2f <a1< %4.2f and %4.2f <a2< %4.2f)",fPrefixLDC.Data(),nGoodChannel,goodA1Min,goodA1Max, goodA2Min,goodA2Max); | |
719 | (*fFilcout) << detail ; printf("%s",detail); | |
720 | detail=Form("\n%s : Nb of uncalibrated channel = %d [%6.4f] (%d unfitted channels) ",fPrefixLDC.Data(),nBadChannel,ratio,noFitChannel); | |
721 | (*fFilcout) << detail ; printf("%s",detail); | |
722 | ||
723 | if(ratio > ratio_limit) { status=-1; | |
724 | detail=Form("\n%s : !!!!! WARNING : Nb of uncalibrated channels very large : %6.4f > %6.4f (status= %d) ",fPrefixLDC.Data(),ratio,ratio_limit,status); | |
725 | (*fFilcout) << detail ; printf("%s",detail); | |
726 | } | |
727 | ||
0e5fb403 | 728 | Double_t meanA1 = sumA1/(nGoodChannel); |
729 | Double_t meanProbChi2 = sumProbChi2/(nGoodChannel); | |
730 | Double_t meanA2 = sumA2/(nGoodChannel); | |
731 | Double_t meanProbChi2P2 = sumProbChi2P2/(nGoodChannel); | |
0e5fb403 | 732 | Double_t capaManu = 0.2; // pF |
9c2d4e89 | 733 | (*fFilcout) << "\n linear fit : <a1> = " << meanA1 << " Prob(chi2)> = " << meanProbChi2 << " <gain> = " << 1./(meanA1*capaManu) |
0e5fb403 | 734 | << " mV/fC (capa= " << capaManu << " pF)" << endl; |
9c2d4e89 | 735 | (*fFilcout) <<" parabolic fit: <a2> = " << meanA2 << " Prob(chi2)> = " << meanProbChi2P2 << "\n" << endl; |
736 | ||
737 | detail=Form("\n%s : <gain> = %7.5f mV/fC (capa= %3.1f pF) Prob(chi2) = %5.3f\n" ,fPrefixLDC.Data(),1./(meanA1*capaManu),capaManu,meanProbChi2); | |
738 | printf("%s",detail); | |
0e5fb403 | 739 | } |
740 | else | |
9c2d4e89 | 741 | { status=-1; |
742 | detail=Form("\n%s : !!!!! ERROR : Nb of Manu = %d or Nb calibrated channel = %d !!!!! (status= %d)\n",fPrefixLDC.Data(),nmanu,nGoodChannel,status); | |
743 | (*fFilcout) << detail ; printf("%s",detail); | |
0e5fb403 | 744 | } |
5431405e | 745 | pfilew.close(); |
5431405e | 746 | |
5253c153 | 747 | if(fPlotLevel>0){tg->Write();histoFile->Close();} |
748 | if(fPrintLevel>1){fclose(pfilep); fclose(pfilen);} | |
9c2d4e89 | 749 | SetStatusDA(status); |
5253c153 | 750 | } |