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