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