- Channels (or pads) with bad addressing of manus or occupancy rate > 1 are not calib...
[u/mrichter/AliRoot.git] / MUON / AliMUONGain.cxx
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 //AMORE
39 //
40 #ifdef ALI_AMORE
41 #include <AmoreDA.h>
42 #include "TSystem.h"
43 #endif
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 //______________________________________________________________________________
65 Double_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 //______________________________________________________________________________
72 Double_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 //______________________________________________________________________________
79 Double_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
92 ClassImp(AliMUONGain)
93 /// \endcond
94
95 //______________________________________________________________________________
96 AliMUONGain::AliMUONGain()
97 : AliMUONPedestal(),
98 fInjCharge(0), 
99 fnInit(1),
100 fnEntries(11),
101 fnbpf1(6),
102 fPrintLevel(0), 
103 fPlotLevel(0)
104 {
105 /// Default constructor
106
107   sprintf(fRootDataFileName," "); //Gain
108 }
109 //  AliMUONPedestal& operator=(const AliMUONPedestal& other); Copy ctor
110   
111 //______________________________________________________________________________
112 AliMUONGain::~AliMUONGain()
113 {
114 /// Destructor
115 }
116
117 //______________________________________________________________________________
118 TString 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 //______________________________________________________________________________
136 void 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)
173       sprintf(flatFile,"%s_%d.ped",fprefixDA,fRunNumber);
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 //______________________________________________________________________________
214 TString AliMUONGain::WriteGainHeader(Int_t nInit, Int_t nEntries, Int_t nbpf2, Int_t *numrun, Double_t *injCharge) 
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
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;
235
236   stream<<"//   RUN     DAC   " << endl;
237   stream<<"//-----------------" << endl;
238   for (Int_t i = 0; i < nEntries; ++i) {
239         stream<<Form("//   %d    %5.0f \n",numrun[i],injCharge[i]);
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 //______________________________________________________________________________
249 TString 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("");
254   stream << Form("%4i %5i %2i %7.4f %10.3e %4i %2x\n",BP,Manu,ch,p1,p2,threshold,q);
255   return TString(stream.str().c_str());
256
257 }
258
259 //__________
260 void 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       {
345         sprintf(filename,"%s_%d.param",fprefixDA,fRunNumber);
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
356
357
358
359   //  Write Ascii file -> Shuttle
360
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
368
369   // print mean and sigma values in file
370   FILE *pfilep = 0;
371   if(fPrintLevel>1)
372     {
373       sprintf(filename,"%s_%d.peak",fprefixDA,fRunNumber);
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     {
408       sprintf(fHistoFileName,"%s_%d.root",fprefixDA,fRunNumber);
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];
505               if(x[j]==0.)fitproceed=0;
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
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];
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             }
632           gain=1./(par[1]*capa); // mv/fC
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)
639                             if(p1>10 && p2>10  and  nplot < 100)
640                   //    if(p1>=1 and p1<=2  and  nplot < 100)
641 //                if((p1==1 || p2==1) and  nplot < 100)
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
673           tempstring = WriteGainData(busPatchId,manuId,channelId,par[1],par[2],threshold,q);
674           pfilew << tempstring;
675 #ifdef ALI_AMORE
676           stringout << tempstring;
677 #endif
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
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
754   if(fPlotLevel>0){tg->Write();histoFile->Close();}
755   if(fPrintLevel>1){fclose(pfilep); fclose(pfilen);}
756 }