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