Adding PiKP-only histograms and eliminating a number of switches where histograms...
[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 #include <cstdio>
38
39 #define  NFITPARAMS 4
40
41 //-----------------------------------------------------------------------------
42 /// \class AliMUONGain
43 ///
44 /// Implementation of calibration computing
45 ///
46 /// add
47 /// 
48 ///
49 /// \author Alberto Baldisseri, JL Charvet 
50 //-----------------------------------------------------------------------------
51
52
53
54
55 // functions
56
57 namespace {
58   
59   //______________________________________________________________________________
60   Double_t funcLin (const Double_t *x, const Double_t *par)
61   {
62     /// Linear function
63     return par[0] + par[1]*x[0];
64   }
65   
66   //______________________________________________________________________________
67   Double_t funcParabolic (const Double_t *x, const Double_t *par)
68   {
69     /// Parabolic function
70     return par[0]*x[0]*x[0];
71   }
72   
73   //______________________________________________________________________________
74   Double_t funcCalib (const Double_t *x, const Double_t *par)  
75   {
76     /// Calibration function
77     Double_t xLim= par[3];
78     
79     if(x[0] <= xLim) return par[0] + par[1]*x[0];
80     
81     Double_t yLim = par[0]+ par[1]*xLim;
82     return yLim + par[1]*(x[0] - xLim) + par[2]*(x[0] - xLim)*(x[0] - xLim);
83   }
84   
85 }
86
87 /// \cond CLASSIMP
88 ClassImp(AliMUONGain)
89 /// \endcond
90
91 //______________________________________________________________________________
92 AliMUONGain::AliMUONGain()
93 : AliMUONPedestal(),
94 fInjCharge(0), 
95 fRootDataFileName(),
96 fnInit(1),
97 fnEntries(11),
98 fnbpf1(6),
99 fPrintLevel(0), 
100 fPlotLevel(0)
101 {
102 /// Default constructor
103
104 }
105   
106 //______________________________________________________________________________
107 AliMUONGain::~AliMUONGain()
108 {
109 /// Destructor
110 }
111
112 //______________________________________________________________________________
113 TString AliMUONGain::WriteDummyHeader(void) 
114 {
115 ///
116
117   ostringstream stream;
118   stream<<"//DUMMY FILE (to prevent Shuttle failure)"<< endl;
119   stream<<"//================================================" << endl;
120   stream<<"//       MUONTRKGAINda: Calibration run  " << endl;
121   stream<<"//================================================" << endl;
122   stream<<"//   * Run           : " << fRunNumber << endl; 
123   stream<<"//   * Date          : " << fDate->AsString("l") <<endl;
124   stream<<"//   * DAC           : " << fInjCharge << endl;
125   stream<<"//-------------------------------------------------" << endl;
126
127   return TString(stream.str().c_str());
128 }
129
130 //______________________________________________________________________________
131 void AliMUONGain::MakePedStoreForGain(TString shuttleFile)
132 {
133 /// Store Pedmean and sigma to pedestal-like ascii file
134
135   ofstream fileout;
136   TString tempstring;
137   TString flatFile;
138   TString outputFile;
139   
140   // Store pedestal map in root file
141   TTree* tree = 0x0;
142
143   // write dummy ascii file -> Shuttle
144   if(fIndex<fnEntries)
145     {  
146       FILE *pfilew=0;
147       pfilew = fopen (shuttleFile,"w");
148
149       fprintf(pfilew,"//DUMMY FILE (to prevent Shuttle failure)\n");
150       fprintf(pfilew,"//================================================\n");
151       fprintf(pfilew,"//       MUONTRKGAINda: Calibration run  \n");
152       fprintf(pfilew,"//=================================================\n");
153       fprintf(pfilew,"//   * Run           : %d \n",fRunNumber); 
154       fprintf(pfilew,"//   * Date          : %s \n",fDate->AsString("l"));
155       fprintf(pfilew,"//   * DAC           : %d \n",fInjCharge);
156       fprintf(pfilew,"//-------------------------------------------------\n");
157       fclose(pfilew);
158     }
159
160
161
162   Finalize();
163   MakeControlHistos();
164   if(fPrintLevel>0)
165     {
166       // compute and store mean DAC values (like pedestals)
167       flatFile = Form("%s.ped",fPrefixDA.Data());
168       outputFile=flatFile;
169       cout << "\n" << fPrefixLDC.Data() << " : Flat file  generated  : " << flatFile.Data() << "\n";
170       if (!outputFile.IsNull())  
171       {
172         ofstream out(outputFile.Data());
173         MakeASCIIoutput(out);
174         out.close();
175       }      
176     }
177
178   TString mode("UPDATE");
179
180   if (fIndex==1) {
181     mode = "RECREATE";
182   }
183   TFile* histoFile = new TFile(fRootDataFileName.Data(), mode.Data(), "MUON Tracking Gains");
184
185   // second argument should be the injected charge, taken from config crocus file
186   // put also info about run number could be usefull
187   AliMpIntPair* pair   = new AliMpIntPair(fRunNumber,fInjCharge );
188
189   if (mode.CompareTo("UPDATE") == 0) {
190     tree = (TTree*)histoFile->Get("t");
191     tree->SetBranchAddress("run",&pair);
192     tree->SetBranchAddress("ped",&fPedestalStore);
193
194   } else {
195     tree = new TTree("t","Pedestal tree");
196     tree->Branch("run", "AliMpIntPair",&pair);
197     tree->Branch("ped", "AliMUON2DMap",&fPedestalStore);
198     tree->SetBranchAddress("run",&pair);
199     tree->SetBranchAddress("ped",&fPedestalStore);
200
201   }
202
203   tree->Fill();
204   tree->Write("t", TObject::kOverwrite); // overwrite the tree
205   histoFile->Close();
206
207   delete pair;
208 }
209
210 //______________________________________________________________________________
211 TString AliMUONGain::WriteGainHeader(Int_t nInit, Int_t nEntries, Int_t nbpf2, Int_t *numrun, Double_t *injCharge) 
212 {
213 /// Header of the calibration output file
214
215   ostringstream stream;
216
217
218   stream<<"//=======================================================" << endl;
219   stream<<"//      Calibration file calculated by " << fPrefixDA.Data() <<endl;
220   stream<<"//=======================================================" << endl;
221   stream<<"//   * Run           : " << fRunNumber << endl; 
222   stream<<"//   * Date          : " << fDate->AsString("l") <<endl;
223   stream<<"//   * Statictics    : " << fNEvents << endl;
224   if(fConfig)
225   stream<<"//   * # of MANUS    : " << fNManuConfig << " read in the Det. config. " << endl;
226   stream<<"//   * # of MANUS    : " << fNManu << " read in raw data " << endl;
227   stream<<"//   * # of MANUS    : " << fNChannel/64 << " written in calibration file " << 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   stream<<"//   *  nInit = " << nInit << "  *  f1nbp = " << fnbpf1 << "  *  f2nbp = " <<  nbpf2 << endl; 
236
237   stream<<"//   RUN     DAC   " << endl;
238   stream<<"//-----------------" << endl;
239   for (Int_t i = 0; i < nEntries; ++i) {
240         stream<<Form("//   %d    %5.0f \n",numrun[i],injCharge[i]);
241   }
242   stream<<"//=======================================" << endl;
243   stream<<"// BP MANU CH.   a1      a2     thres. q " << endl;
244   stream<<"//=======================================" << endl;
245
246   return TString(stream.str().c_str());
247 }
248
249 //______________________________________________________________________________
250 TString AliMUONGain::WriteGainData(Int_t BP, Int_t Manu, Int_t ch, Double_t p1, Double_t p2, Int_t threshold, Int_t q) 
251 {
252 /// Write calibration parameters per channel
253
254   ostringstream stream("");
255   stream << Form("%4i %5i %2i %7.4f %10.3e %4i %2x\n",BP,Manu,ch,p1,p2,threshold,q);
256   return TString(stream.str().c_str());
257
258 }
259
260 //_______________________________________________________________________________
261 void AliMUONGain::MakeGainStore(TString shuttleFile)
262 {
263   Int_t status=0;
264   /// Store gains in ASCII files
265   ofstream fileout;
266   ofstream filcouc;
267   TString tempstring;   
268   TString filename; 
269   char* detail;
270
271   Double_t goodA1Min =  0.5;
272   Double_t goodA1Max =  2.;
273 //   Double_t goodA2Min = -0.5E-03;
274 //   Double_t goodA2Max =  1.E-03;
275   Double_t goodA2Min = -0.5E-01; // changed 28/10/2009 (JLC) <=> enlarged condition on a2
276   Double_t goodA2Max =  1.E-01;
277   // Table for uncalibrated  buspatches and manus
278   THashList* uncalBuspatchManuTable = new THashList(1000,2);
279
280   Int_t numrun[11];
281
282   // open file MUONTRKda_gain_data.root
283   // read again the pedestal for the calibration runs (11 runs)
284   // need the injection charge from config file (to be done)
285   // For each channel make a TGraphErrors (mean, sigma) vs injected charge
286   // Fit with a polynomial fct
287   // store the result in a flat file.
288
289   if(fIndex==0)cout << " Root data file = " << fRootDataFileName.Data() << endl;  
290   TFile*  histoFile = new TFile(fRootDataFileName.Data());
291
292   AliMUON2DMap* map[11];
293   AliMUONVCalibParam* ped[11];
294   AliMpIntPair* run[11];
295
296   //read back from root file
297   TTree* tree = (TTree*)histoFile->Get("t");
298   Int_t nEntries = tree->GetEntries();
299
300   Int_t nbpf2 = nEntries - (fnInit + fnbpf1) + 1; // nb pts used for 2nd order fit
301
302   // read back info
303   for (Int_t i = 0; i < nEntries; ++i) {
304     map[i] = 0x0;
305     run[i] = 0x0;
306     tree->SetBranchAddress("ped",&map[i]);
307     tree->SetBranchAddress("run",&run[i]);
308     tree->GetEvent(i);
309     //        std::cout << map[i] << " " << run[i] << std::endl;
310   }
311   // RunNumber extracted from Root data fil
312   if(fIndex==0)fRunNumber=(UInt_t)run[nEntries-1]->GetFirst();
313   //     sscanf(getenv("DATE_RUN_NUMBER"),"%d",&gAliRunNumber);
314
315   Double_t pedMean[11];
316   Double_t pedSigma[11];
317   for ( Int_t i=0 ; i<11 ; i++) {pedMean[i]=0.;pedSigma[i]=1.;};
318   Double_t injCharge[11];
319   Double_t injChargeErr[11];
320   for ( Int_t i=0 ; i<11 ; i++) {injCharge[i]=0.;injChargeErr[i]=1.;};
321
322   //  print out in .log file
323   detail=Form("\n%s : ------  MUONTRKGAINda for Gain computing (Last Run = %d) ------\n",fPrefixLDC.Data(),fRunNumber); printf("%s",detail);
324   (*fFilcout)<<"\n\n//=================================================" << endl;
325   (*fFilcout)<<"//    MUONTRKGAINda: Gain Computing  Run = " << fRunNumber << endl;
326   (*fFilcout)<<"//    RootDataFile  = "<< fRootDataFileName.Data() << endl;
327   (*fFilcout)<<"//=================================================" << endl;
328   (*fFilcout)<<"//* Date          : " << fDate->AsString("l") << "\n" << endl;
329
330   (*fFilcout) << " Entries = " << nEntries << " DAC values \n" << endl; 
331   for (Int_t i = 0; i < nEntries; ++i) {
332     (*fFilcout) << " Run = " << run[i]->GetFirst() << "    DAC = " << run[i]->GetSecond() << endl;
333     numrun[i] = run[i]->GetFirst();
334     injCharge[i] = run[i]->GetSecond();
335     injChargeErr[i] = 0.01*injCharge[i];
336     if(injChargeErr[i] <= 1.) injChargeErr[i]=1.;
337   }
338   detail=Form("%s : .... Fitting .... \n",fPrefixLDC.Data()); printf("%s",detail);
339
340
341   // why 2 files ? (Ch. F.)  => second file contains detailed results
342     FILE *pfilen = 0;
343     if(fPrintLevel>1)
344       {
345         filename=Form("%s.param",fPrefixDA.Data());
346         detail=Form("%s : Second fit parameter file        = %s\n",fPrefixLDC.Data(),filename.Data()); printf("%s",detail);
347         //       cout << " Second fit parameter file        = " << filename.Data() << "\n";
348         pfilen = fopen (filename.Data(),"w");
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
357
358   // file outputs for gain
359
360   ofstream pfilew;
361   pfilew.open(shuttleFile.Data());
362   // Write Header Data of the .par file
363   pfilew << WriteGainHeader(fnInit,nEntries,nbpf2,numrun,injCharge);
364
365   // print mean and sigma values in file
366   FILE *pfilep = 0;
367   if(fPrintLevel>1)
368     {
369       filename=Form("%s.peak",fPrefixDA.Data());
370       detail=Form("%s : File containing Peak mean values = %s\n",fPrefixLDC.Data(),filename.Data()); printf("%s",detail);
371       //      cout << " File containing Peak mean values = " << filename << "\n";
372       pfilep = fopen (filename,"w");
373
374       fprintf(pfilep,"//==============================================================================================================================\n");
375       fprintf(pfilep,"//   * Run           : %d \n",fRunNumber); 
376       fprintf(pfilep,"//==============================================================================================================================\n");
377       fprintf(pfilep,"// BP  MANU  CH.    Ped.     <0>      <1>      <2>      <3>      <4>      <5>      <6>      <7>      <8>      <9>     <10> \n"); 
378       fprintf(pfilep,"//==============================================================================================================================\n");
379       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]);
380       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]);
381       fprintf(pfilep,"//==============================================================================================================================\n");
382     }
383
384   Double_t chi2    = 0.;
385   Double_t chi2P2  = 0.;
386   Double_t prChi2  = 0; 
387   Double_t prChi2P2 =0;
388   Double_t a0=0.,a1=1.,a2=0.;
389   Int_t busPatchId ;
390   Int_t manuId     ;
391   Int_t channelId ;
392   Int_t threshold = 0;
393   Int_t q = 0;
394   Int_t p1 =0;
395   Int_t p2 =0;
396   Double_t gain=10.; // max value (= bad value)
397   Double_t capa=0.2; // internal capacitor (pF)
398
399   //  plot out 
400
401   TTree* tg = 0x0;
402   if(fPlotLevel>0)
403     {
404       fHistoFileName=Form("%s.root",fPrefixDA.Data());
405       new TFile(fHistoFileName.Data(),"RECREATE","MUON Tracking gains");
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 ;
444   Int_t unparabolicfit=0;
445
446   while ( ( p = dynamic_cast<AliMUONVCalibParam*>(next() ) ) )
447     {
448       ped[0]  = p;
449       unparabolicfit=0;
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
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           // Fit Method:  Linear fit over gAlinbpf1 points + parabolic fit  over nbpf2  points) 
487           // nInit=1 : 1st pt DAC=0 excluded
488
489           // 1. - linear fit over gAlinbpf1 points
490
491           Double_t par[4] = {0.,0.5,0.,ADCMax()};
492           Int_t nbs   = nEntries - fnInit;
493           if(nbs < fnbpf1)fnbpf1=nbs;
494
495           Int_t fitproceed=1;
496           Int_t nbpf2Dynamic=nbpf2;
497           Int_t adcLimit=4090; // when RMS < 0.5 (in other cases mean values forced to 4095, see DA_PED)
498           for (Int_t j = 0; j < nbs; ++j)
499             {
500               Int_t k = j + fnInit;
501               x[j]    = pedMean[k];
502               if(x[j]<=0.){fitproceed=0; break;}
503               //              if(x[j]>= ADCMax())
504               if(x[j]>= adcLimit)
505                 {
506                   if(j < nbs-1){fitproceed=0; break;}
507                   else { nbpf2Dynamic=nbpf2-1; break;}
508                 }
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                       
521               TF1 *f1 = new TF1("f1",funcLin,0.,ADCMax(),2);
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
544               //checking:         if(busPatchId ==1841 && manuId==4)nbpf2Dynamic=2;
545               if(nbpf2Dynamic > 2)
546                 {
547                   for (Int_t j = 0; j < nbpf2Dynamic; j++)
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
557                   TF1 *f2 = new TF1("f2",funcParabolic,0.,ADCMax(),1);
558                   graphErr = new TGraphErrors(nbpf2Dynamic, xp, yp, xpErr, ypErr);
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
568                   prChi2P2 = TMath::Prob(chi2P2, nbpf2Dynamic-1);
569                   a2 = par[0];
570                 }
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                 }
578
579               par[0] = a0;
580               par[1] = a1;
581               par[2] = a2;
582               par[3] = xLim;
583
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];
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;
620               threshold=ADCMax();       
621
622               // bad calibration counter
623               char bpmanuname[256];
624               AliMUONErrorCounter* uncalcounter;
625
626               snprintf(bpmanuname,256,"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             }
642           gain=1./(par[1]*capa); // mv/fC
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)
649                   //                        if(p1>10 && p2>10  and  nplot < 100)
650                   //    if(p1>=1 and p1<=2  and  nplot < 100)
651 //                if((p1==1 || p2==1) and  nplot < 100)
652                             if(nbpf2Dynamic<nbpf2  and  nplot < 100)
653                     {
654                       nplot++;
655                       //              cout << " nplot = " << nplot << endl;
656                       TF1 *f2Calib = new TF1("f2Calib",funcCalib,0.,ADCMax(),NFITPARAMS);
657
658                       graphErr = new TGraphErrors(nEntries,pedMean,injCharge,pedSigma,injChargeErr);
659
660                       snprintf(graphName,256,"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                       snprintf(graphName,256,"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                 }
678               tg->Fill();
679             }
680           pfilew << WriteGainData(busPatchId,manuId,channelId,par[1],par[2],threshold,q);
681         }
682       nmanu++;
683       Int_t step=500;
684       if(nmanu % step == 0)printf("%s : Nb manu = %d\n",fPrefixLDC.Data(),nmanu);
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   if(nmanu && nGoodChannel) 
708     {
709       Double_t ratio_limit=0.25;
710       Double_t ratio=float (nBadChannel)/float (nmanu*64);
711
712       detail=Form("\n%s : Nb of channels in raw data = %d (%d Manu)",fPrefixLDC.Data(),nmanu*64,nmanu); 
713       (*fFilcout) << detail ; printf("%s",detail);
714       detail=Form("\n%s : Nb of calibrated channel   = %d  (%4.2f <a1< %4.2f and %4.2f <a2< %4.2f)",fPrefixLDC.Data(),nGoodChannel,goodA1Min,goodA1Max, goodA2Min,goodA2Max);
715       (*fFilcout) << detail ; printf("%s",detail);
716       detail=Form("\n%s : Nb of uncalibrated channel = %d  [%6.4f] (%d unfitted channels) ",fPrefixLDC.Data(),nBadChannel,ratio,noFitChannel);
717       (*fFilcout) << detail ; printf("%s",detail);
718
719       if(ratio > ratio_limit) { status=-1;
720         detail=Form("\n%s : !!!!! WARNING : Nb of uncalibrated channels very large : %6.4f > %6.4f (status= %d) ",fPrefixLDC.Data(),ratio,ratio_limit,status);
721         (*fFilcout) << detail ; printf("%s",detail); 
722       }
723
724       Double_t meanA1         = sumA1/(nGoodChannel);
725       Double_t meanProbChi2   = sumProbChi2/(nGoodChannel);
726       Double_t meanA2         = sumA2/(nGoodChannel);
727       Double_t meanProbChi2P2 = sumProbChi2P2/(nGoodChannel);
728       Double_t capaManu = 0.2; // pF
729       (*fFilcout) << "\n linear fit   : <a1> = " << meanA1 << "     Prob(chi2)>  = " <<  meanProbChi2 <<  "    <gain>  = " <<  1./(meanA1*capaManu) 
730                   << " mV/fC (capa= " << capaManu << " pF)" << endl;
731       (*fFilcout)   <<" parabolic fit: <a2> = " << meanA2  << "    Prob(chi2)>  = " <<  meanProbChi2P2 << "\n" <<  endl;
732  
733       detail=Form("\n%s : <gain>  = %7.5f  mV/fC (capa= %3.1f pF)  Prob(chi2) = %5.3f\n" ,fPrefixLDC.Data(),1./(meanA1*capaManu),capaManu,meanProbChi2);
734       printf("%s",detail);
735     }
736   else 
737     { status=-1;
738       detail=Form("\n%s : !!!!! ERROR :  Nb of Manu = %d or Nb calibrated channel = %d !!!!! (status= %d)\n",fPrefixLDC.Data(),nmanu,nGoodChannel,status);
739       (*fFilcout) << detail ; printf("%s",detail);
740     }
741   pfilew.close();
742
743   if(fPlotLevel>0){tg->Write();histoFile->Close();}
744   if(fPrintLevel>1){fclose(pfilep); fclose(pfilen);}
745   SetStatusDA(status);
746 }