o add Reset function to CalPad and CalROC o Add functionality to AliTPCdataQA - Reset...
[u/mrichter/AliRoot.git] / MUON / AliMUONGain.cxx
CommitLineData
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 57namespace {
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 87using std::endl;
88using std::cout;
89using std::istringstream;
90using std::ostringstream;
5253c153 91/// \cond CLASSIMP
92ClassImp(AliMUONGain)
93/// \endcond
94
95//______________________________________________________________________________
96AliMUONGain::AliMUONGain()
97: AliMUONPedestal(),
98fInjCharge(0),
f3cfa63e 99fRootDataFileName(),
5253c153 100fnInit(1),
101fnEntries(11),
102fnbpf1(6),
103fPrintLevel(0),
104fPlotLevel(0)
105{
106/// Default constructor
107
5253c153 108}
5253c153 109
110//______________________________________________________________________________
111AliMUONGain::~AliMUONGain()
112{
113/// Destructor
114}
115
116//______________________________________________________________________________
117TString 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//______________________________________________________________________________
135void 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 215TString 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//______________________________________________________________________________
254TString 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 265void 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}