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