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 #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 //______________________________________________________________________________
  58 Double_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 //______________________________________________________________________________
  65 Double_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 //______________________________________________________________________________
  72 Double_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
  85 ClassImp(AliMUONGain)
  86 /// \endcond
  87
  88 //______________________________________________________________________________
  89 AliMUONGain::AliMUONGain()
  90 : AliMUONPedestal(),
  91 fInjCharge(0),
  92 fnInit(1),
  93 fnEntries(11),
  94 fnbpf1(6),
  95 fPrintLevel(0),
  96 fPlotLevel(0)
  97 {
  98 /// Default constructor
  99
 100   sprintf(fRootDataFileName," "); //Gain
 101 }
 102 //  AliMUONPedestal& operator=(const AliMUONPedestal& other); Copy ctor
 103
 104 //______________________________________________________________________________
 105 AliMUONGain::~AliMUONGain()
 106 {
 107 /// Destructor
 108 }
 109
 110 //______________________________________________________________________________
 111 TString 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 //______________________________________________________________________________
 129 void 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 //______________________________________________________________________________
 208 TString 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 //______________________________________________________________________________
 245 TString 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 //__________
 257 void 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 }