MUON/MUONTRKda.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 /* 21/11/07
  19
  20 Version for MUONTRKda MUON tracking
  21 (A. Baldisseri, J.-L. Charvet & Ch. Finck)
  22
  23
  24 Rem:  AliMUON2DMap stores all channels, even those which are not connected
  25       if pedMean == -1, channel not connected to a pad
  26
  27
  28 */
  29 extern "C" {
  30 #include <daqDA.h>
  31 }
  32
  33 #include "event.h"
  34 #include "monitor.h"
  35
  36 #include <Riostream.h>
  37 #include <stdio.h>
  38 #include <stdlib.h>
  39
  40 //AliRoot
  41 #include "AliMUONLogger.h"
  42 #include "AliMUONRawStreamTracker.h"
  43 #include "AliMUONDspHeader.h"
  44 #include "AliMUONBlockHeader.h"
  45 #include "AliMUONBusStruct.h"
  46 #include "AliMUONDDLTracker.h"
  47 #include "AliMUONVStore.h"
  48 #include "AliMUON2DMap.h"
  49 #include "AliMUONCalibParamND.h"
  50 #include "AliMpIntPair.h"
  51 #include "AliMpConstants.h"
  52 #include "AliRawReaderDate.h"
  53
  54 //ROOT
  55 #include "TFile.h"
  56 #include "TSystem.h"
  57 #include "TTree.h"
  58 #include "TH1F.h"
  59 #include "TString.h"
  60 #include "TStopwatch.h"
  61 #include "TMath.h"
  62 #include "TTimeStamp.h"
  63 #include "TGraphErrors.h"
  64 #include "TF1.h"
  65 #include "TROOT.h"
  66 #include "TPluginManager.h"
  67 #include "TFitter.h"
  68
  69 #define  NFITPARAMS 4
  70
  71 // global variables
  72 const Int_t gkNChannels = AliMpConstants::ManuNofChannels();
  73 const Int_t gkADCMax    = 4095;
  74
  75 AliMUONVStore* gPedestalStore =  new AliMUON2DMap(kFALSE);
  76
  77 Int_t  gNManu       = 0;
  78 Int_t  gNChannel    = 0;
  79 UInt_t gRunNumber   = 0;
  80 Int_t  gNEvents     = 0;
  81 Int_t  gGlitchErrors= 0;
  82 Int_t  gParityErrors= 0;
  83 Int_t  gPaddingErrors= 0;
  84 Int_t  gNDateEvents = 0;
  85 Int_t  gPrintLevel  = 1;  // global printout variable
  86 Int_t  gPlotLevel  = 0;  // global plot variable
  87
  88 TH1F*  gPedMeanHisto  = 0x0;
  89 TH1F*  gPedSigmaHisto = 0x0;
  90 Char_t gHistoFileName[256];
  91
  92 // used by makegain
  93 Char_t gHistoFileName_gain[256]="MUONTRKda_gain_data.root";
  94 Char_t gRootFileName[256];
  95 Char_t gOutFolder[256]=".";
  96 Char_t filename[256];
  97 Char_t filenam[256]="MUONTRKda_gain";
  98 // Char_t filenam[256]="MUONTRKda_gain.param"; // if gPrintLevel  = 2
  99 // Char_t filenam_bad[256]="MUONTRKda_gain.bad"; // if gPrintLevel  = 2
 100 // Char_t filenam_peak[256]="MUONTRKda_gain.peak"; // if gPrintLevel  = 3
 101 // Char_t filenam_out[256]="MUONTRKda_gain.out";
 102 Char_t flatFile[256];
 103
 104 ofstream filcout;
 105
 106 TString flatOutputFile;
 107 TString logOutputFile;
 108 TString gCommand("ped");
 109 TTimeStamp date;
 110
 111 // funtions
 112
 113
 114 //________________
 115 Double_t funcLin(Double_t *x, Double_t *par)
 116 {
 117   return par[0] + par[1]*x[0];
 118 }
 119
 120 //________________
 121 Double_t funcParabolic(Double_t *x, Double_t *par)
 122 {
 123   return par[0]*x[0]*x[0];
 124 }
 125
 126 //________________
 127 Double_t funcCalib(Double_t *x, Double_t *par)
 128 {
 129   Double_t xLim= par[3];
 130
 131   if(x[0] <= xLim) return par[0] + par[1]*x[0];
 132
 133   Double_t yLim = par[0]+ par[1]*xLim;
 134   return yLim + par[1]*(x[0] - xLim) + par[2]*(x[0] - xLim)*(x[0] - xLim);
 135 }
 136
 137
 138 //__________
 139 void MakePed(Int_t busPatchId, Int_t manuId, Int_t channelId, Int_t charge)
 140 {
 141
 142     AliMUONVCalibParam* ped =
 143         static_cast<AliMUONVCalibParam*>(gPedestalStore->FindObject(busPatchId, manuId));
 144
 145     if (!ped) {
 146       gNManu++;
 147       ped = new AliMUONCalibParamND(2, gkNChannels,busPatchId, manuId, -1.); // put default wise -1, not connected channel
 148       gPedestalStore->Add(ped);
 149     }
 150
 151     if (gNEvents == 0) {
 152       ped->SetValueAsDouble(channelId, 0, 0.);
 153       ped->SetValueAsDouble(channelId, 1, 0.);
 154     }
 155
 156     Double_t pedMean  = ped->ValueAsDouble(channelId, 0) + charge;
 157     Double_t pedSigma = ped->ValueAsDouble(channelId, 1) + charge*charge;
 158
 159     ped->SetValueAsDouble(channelId, 0, pedMean);
 160     ped->SetValueAsDouble(channelId, 1, pedSigma);
 161
 162 }
 163
 164 //________________
 165 void MakePedStore(TString flatOutputFile = "")
 166 {
 167 //   TTimeStamp date;
 168   Double_t pedMean;
 169   Double_t pedSigma;
 170   ofstream fileout;
 171   Int_t busPatchId;
 172   Int_t manuId;
 173   Int_t channelId;
 174
 175  // histo
 176   TFile*  histoFile = 0;
 177   TTree* tree = 0;
 178   if (gCommand.CompareTo("ped") == 0)
 179     {
 180       sprintf(gHistoFileName,"%s/MUONTRKda_ped_%d.root",gOutFolder,gRunNumber);
 181       histoFile = new TFile(gHistoFileName,"RECREATE","MUON Tracking pedestals");
 182 //     }
 183
 184   Char_t name[255];
 185   Char_t title[255];
 186   sprintf(name,"pedmean_allch");
 187   sprintf(title,"Pedestal mean all channels");
 188   Int_t nx = 4096;
 189   Int_t xmin = 0;
 190   Int_t xmax = 4095;
 191   gPedMeanHisto = new TH1F(name,title,nx,xmin,xmax);
 192   gPedMeanHisto->SetDirectory(histoFile);
 193
 194   sprintf(name,"pedsigma_allch");
 195   sprintf(title,"Pedestal sigma all channels");
 196   nx = 201;
 197   xmin = 0;
 198   xmax = 200;
 199   gPedSigmaHisto = new TH1F(name,title,nx,xmin,xmax);
 200   gPedSigmaHisto->SetDirectory(histoFile);
 201
 202 //   TTree* tree = new TTree("t","Pedestal tree");
 203   tree = new TTree("t","Pedestal tree");
 204   tree->Branch("bp",&busPatchId,"bp/I");
 205   tree->Branch("manu",&manuId,",manu/I");
 206   tree->Branch("channel",&channelId,",channel/I");
 207   tree->Branch("pedMean",&pedMean,",pedMean/D");
 208   tree->Branch("pedSigma",&pedSigma,",pedSigma/D");
 209
 210     }
 211
 212   if (!flatOutputFile.IsNull()) {
 213     fileout.open(flatOutputFile.Data());
 214     fileout<<"//===========================================================================" << endl;
 215     fileout<<"//                       Pedestal file calculated by MUONTRKda"<<endl;
 216     fileout<<"//===========================================================================" << endl;
 217     fileout<<"//       * Run           : " << gRunNumber << endl;
 218     fileout<<"//       * Date          : " << date.AsString("l") <<endl;
 219     fileout<<"//       * Statictics    : " << gNEvents << endl;
 220     fileout<<"//       * # of MANUS    : " << gNManu << endl;
 221     fileout<<"//       * # of channels : " << gNChannel << endl;
 222     fileout<<"//"<<endl;
 223     fileout<<"//---------------------------------------------------------------------------" << endl;
 224     fileout<<"//---------------------------------------------------------------------------" << endl;
 225     fileout<<"//      BP     MANU     CH.      MEAN    SIGMA"<<endl;
 226     fileout<<"//---------------------------------------------------------------------------" << endl;
 227
 228   }
 229
 230   // iterator over pedestal
 231   TIter next(gPedestalStore->CreateIterator());
 232   AliMUONVCalibParam* ped;
 233
 234   while ( ( ped = dynamic_cast<AliMUONVCalibParam*>(next() ) ) )
 235   {
 236     busPatchId              = ped->ID0();
 237     manuId                  = ped->ID1();
 238
 239     for (channelId = 0; channelId < ped->Size() ; ++channelId) {
 240
 241       pedMean  = ped->ValueAsDouble(channelId, 0);
 242
 243       if (pedMean > 0) { // connected channels
 244
 245         ped->SetValueAsDouble(channelId, 0, pedMean/(Double_t)gNEvents);
 246
 247         pedMean  = ped->ValueAsDouble(channelId, 0);
 248         pedSigma = ped->ValueAsDouble(channelId, 1);
 249
 250         ped->SetValueAsDouble(channelId, 1, TMath::Sqrt(TMath::Abs(pedSigma/(Double_t)gNEvents - pedMean*pedMean)));
 251
 252         pedMean  = ped->ValueAsDouble(channelId, 0) + 0.5 ;
 253 //      pedMean  = ped->ValueAsDouble(channelId, 0) ;
 254         pedSigma = ped->ValueAsDouble(channelId, 1);
 255
 256
 257         if (!flatOutputFile.IsNull()) {
 258           fileout << "\t" << busPatchId << "\t" << manuId <<"\t"<< channelId << "\t"
 259                   << pedMean <<"\t"<< pedSigma << endl;
 260         }
 261
 262   if (gCommand.CompareTo("ped") == 0)
 263     {
 264         gPedMeanHisto->Fill(pedMean);
 265         gPedSigmaHisto->Fill(pedSigma);
 266
 267         tree->Fill();
 268       }
 269       }
 270     }
 271   }
 272
 273   // file outputs
 274   if (!flatOutputFile.IsNull())
 275       fileout.close();
 276
 277   if (gCommand.CompareTo("ped") == 0)
 278     {
 279       histoFile->Write();
 280       histoFile->Close();
 281     }
 282
 283 //  delete tree;
 284
 285 }
 286
 287 //________________
 288 void MakePedStoreForGain(Int_t injCharge)
 289 {
 290     // store pedestal map in root file
 291
 292 //     Int_t injCharge = 200;
 293
 294     TTree* tree = 0x0;
 295
 296     // compute and store pedestals
 297     sprintf(flatFile,"%s/%s_%d_DAC_%d.ped",gOutFolder,filenam,gRunNumber,injCharge);
 298     cout << "MUONTRKda : Flat file  generated   = " << flatFile << "\n";
 299     MakePedStore(flatFile);
 300
 301     TString mode("UPDATE");
 302
 303     if (gCommand.Contains("cre")) {
 304         mode = "RECREATE";
 305     }
 306     TFile* histoFile = new TFile(gHistoFileName_gain, mode.Data(), "MUON Tracking Gains");
 307
 308     // second argument should be the injected charge, taken from config crocus file
 309     // put also info about run number could be usefull
 310     AliMpIntPair* pair   = new AliMpIntPair(gRunNumber, injCharge);
 311
 312     if (mode.CompareTo("UPDATE") == 0) {
 313       tree = (TTree*)histoFile->Get("t");
 314       tree->SetBranchAddress("run",&pair);
 315       tree->SetBranchAddress("ped",&gPedestalStore);
 316
 317     } else {
 318       tree = new TTree("t","Pedestal tree");
 319       tree->Branch("run", "AliMpIntPair",&pair);
 320       tree->Branch("ped", "AliMUON2DMap",&gPedestalStore);
 321       tree->SetBranchAddress("run",&pair);
 322       tree->SetBranchAddress("ped",&gPedestalStore);
 323
 324     }
 325
 326     tree->Fill();
 327     tree->Write("t", TObject::kOverwrite); // overwrite the tree
 328     histoFile->Close();
 329
 330     delete pair;
 331 }
 332
 333 //________________
 334 // void MakeGainStore(TString flatOutputFile)
 335 void MakeGainStore()
 336 {
 337     Double_t goodA1Min =  0.5;
 338     Double_t goodA1Max =  2.;
 339     Double_t goodA2Min = -0.5E-03;
 340     Double_t goodA2Max =  1.E-03;
 341
 342     // open file mutrkgain.root
 343     // read again the pedestal for the calibration runs (9 runs ?)
 344     // need the injection charge from config file (to be done)
 345     // For each channel make a TGraphErrors (mean, sigma) vs injected charge
 346     // Fit with a polynomial fct
 347     // store the result in a flat file.
 348
 349
 350     TFile*  histoFile = new TFile(gHistoFileName_gain);
 351
 352     AliMUON2DMap* map[11];
 353     AliMUONVCalibParam* ped[11];
 354     AliMpIntPair* run[11];
 355
 356     //read back from root file
 357     TTree* tree = (TTree*)histoFile->Get("t");
 358     Int_t nEntries = tree->GetEntries();
 359
 360     // read back info
 361     for (Int_t i = 0; i < nEntries; ++i) {
 362       map[i] = 0x0;
 363       run[i] = 0x0;
 364       tree->SetBranchAddress("ped",&map[i]);
 365       tree->SetBranchAddress("run",&run[i]);
 366       tree->GetEvent(i);
 367 //       std::cout << map[i] << " " << run[i] << std::endl;
 368     }
 369     gRunNumber=(UInt_t)run[0]->GetFirst();
 370
 371     // some print
 372     cout<<"\n ********  MUONTRKda for Gain computing (Run = " << gRunNumber << ")\n" << endl;
 373     cout<<" * Date          : " << date.AsString("l") << "\n" << endl;
 374     cout << " Entries = " << nEntries << " DAC values \n" << endl;
 375     for (Int_t i = 0; i < nEntries; ++i) {
 376       cout<< " Run = " << (Double_t)run[i]->GetFirst() << "    DAC = " << (Double_t)run[i]->GetSecond() << endl;
 377     }
 378     cout << "" << endl;
 379
 380
 381     Double_t pedMean[11];
 382     Double_t pedSigma[11];
 383     Double_t injCharge[11];
 384     Double_t injChargeErr[11];
 385
 386 // full print out
 387
 388     sprintf(filename,"%s/%s_%d.log",gOutFolder,filenam,gRunNumber);
 389     logOutputFile=filename;
 390
 391     filcout.open(logOutputFile.Data());
 392     filcout<<"//====================================================" << endl;
 393     filcout<<"//        MUONTRKda for Gain computing (Run = " << gRunNumber << ")" << endl;
 394     filcout<<"//====================================================" << endl;
 395     filcout<<"//   * Date          : " << date.AsString("l") << "\n" << endl;
 396
 397 //     cout<<"\n ********  MUONTRKda for Gain computing (Run = " << gRunNumber << ")\n" << endl;
 398 //     cout<<" * Date          : " << date.AsString("l") << "\n" << endl;
 399
 400
 401
 402     // why 2 files ? (Ch. F.)
 403     FILE *pfilen = 0;
 404     FILE *pfilef = 0;
 405     if(gPrintLevel>=2)
 406     {
 407       sprintf(filename,"%s/%s_%d.param",gOutFolder,filenam,gRunNumber);
 408       cout << " fit parameter file               = " << filename << "\n";
 409       pfilen = fopen (filename,"w");
 410
 411       fprintf(pfilen,"//===================================================================\n");
 412       fprintf(pfilen,"//  BP MANU CH. a0     a1       a2      xlim  P(chi2) P(chi2)2    Q\n");
 413       fprintf(pfilen,"//===================================================================\n");
 414       fprintf(pfilen,"//   * Run           : %d \n",gRunNumber);
 415       fprintf(pfilen,"//===================================================================\n");
 416
 417       sprintf(filename,"%s/%s_%d.bad",gOutFolder,filenam,gRunNumber);
 418       cout << " Bad channel file                 = " << filename << "\n";
 419       pfilef = fopen (filename,"w");
 420
 421       fprintf(pfilef,"//=================================================\n");
 422       fprintf(pfilef,"//  Bad Channel file calculated by MUONTRKda \n");
 423       fprintf(pfilef,"//=================================================\n");
 424       fprintf(pfilef,"//   * Run           : %d \n",gRunNumber);
 425       fprintf(pfilef,"//   * Date          : %s \n",date.AsString("l"));
 426       fprintf(pfilef,"//=======================================\n");
 427       fprintf(pfilef,"// BP MANU CH.   a1      a2     thres. Q\n");
 428       fprintf(pfilef,"//=======================================\n");
 429     }
 430
 431     FILE *pfilew=0;
 432     if(flatOutputFile.IsNull())
 433       {
 434         sprintf(filename,"%s_%d.par",filenam,gRunNumber);
 435         flatOutputFile=filename;
 436       }
 437     if(!flatOutputFile.IsNull())
 438     {
 439       pfilew = fopen (flatOutputFile.Data(),"w");
 440
 441       fprintf(pfilew,"//=================================================\n");
 442       fprintf(pfilew,"//  Calibration file calculated by MUONTRKda \n");
 443       fprintf(pfilew,"//=================================================\n");
 444       fprintf(pfilew,"//   * Run           : %d \n",gRunNumber);
 445       fprintf(pfilew,"//   * Date          : %s \n",date.AsString("l"));
 446       fprintf(pfilew,"//   * Statictics    : %d \n",gNEvents);
 447       fprintf(pfilew,"//   * # of MANUS    : %d \n",gNManu);
 448       fprintf(pfilew,"//   * # of channels : %d \n",gNChannel);
 449       fprintf(pfilew,"//-------------------------------------------------\n");
 450       fprintf(pfilew,"//=======================================\n");
 451       fprintf(pfilew,"// BP MANU CH.   a1      a2     thres. Q\n");
 452       fprintf(pfilew,"//=======================================\n");
 453     }
 454
 455     FILE *pfilep = 0;
 456     if(gPrintLevel==3)
 457     {
 458       sprintf(filename,"%s/%s_%d.peak",gOutFolder,filenam,gRunNumber);
 459       cout << " File containing Peak mean values = " << filename << "\n";
 460       pfilep = fopen (filename,"w");
 461
 462       fprintf(pfilep,"//===============================================================================================================================\n");
 463       fprintf(pfilep,"//   * Run           : %d \n",gRunNumber);
 464       fprintf(pfilep,"//===============================================================================================================================\n");
 465       fprintf(pfilep,"// BP  MANU  CH.    Ped.     <0>      <1>      <2>      <3>      <4>      <5>      <6>      <7>      <8>      <9>     <10> \n");
 466 //       fprintf(pfilep,"//                      %9.1f%9.1f%9.1f%9.1f%9.1f%9.1f%9.1f%9.1f%9.1f%9.1f%9.1f  fC\n",level_fC[0],level_fC[1],level_fC[2],level_fC[3],level_fC[4],level_fC[5],level_fC[6],level_fC[7],level_fC[8],level_fC[9],level_fC[10]);
 467 //       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",level_err[0],level_err[1],level_err[2],level_err[3],level_err[4],level_err[5],level_err[6],level_err[7],level_err[8],level_err[9],level_err[10]);
 468       fprintf(pfilep,"//===============================================================================================================================\n");
 469     }
 470
 471
 472
 473 //  plot out
 474
 475     TFile* gainFile = 0x0;
 476     sprintf(gRootFileName,"%s/%s_%d.root",gOutFolder,filenam,gRunNumber);
 477     gainFile = new TFile(gRootFileName,"RECREATE");
 478
 479     Double_t chi2    = 0.;
 480     Double_t chi2P2  = 0.;
 481     Double_t prChi2  = 0;
 482     Double_t prChi2P2 =0;
 483     Double_t a0,a1,a2;
 484     Int_t busPatchId ;
 485     Int_t manuId     ;
 486     Int_t channelId ;
 487     Int_t threshold = 0;
 488     Int_t Q = 0;
 489     Int_t p1 ;
 490     Int_t p2 ;
 491     Double_t gain;
 492     Double_t capa=0.2; // internal capacitor (pF)
 493
 494     TTree *tg = new TTree("tg","TTree avec class Manu_DiMu");
 495
 496     tg->Branch("bp",&busPatchId, "busPatchId/I");
 497     tg->Branch("manu",&manuId, "manuId/I");
 498     tg->Branch("channel",&channelId, "channelId/I");
 499
 500     tg->Branch("a0",&a0, "a0/D");
 501     tg->Branch("a1",&a1, "a1/D");
 502     tg->Branch("a2",&a2, "a2/D");
 503     tg->Branch("Pchi2",&prChi2, "prChi2/D");
 504     tg->Branch("Pchi2_2",&prChi2P2, "prChi2P2/D");
 505     tg->Branch("Threshold",&threshold, "threshold/I");
 506     tg->Branch("Q",&Q, "Q/I");
 507     tg->Branch("p1",&p1, "p1/I");
 508     tg->Branch("p2",&p2, "p2/I");
 509     tg->Branch("gain",&gain, "gain/D");
 510
 511 // bad BusPatch and manu
 512     Int_t num_tot_BP=800;
 513     Int_t num_tot_Manu=1500;
 514 //     Int_t bad_channel[num_tot_BP][num_tot_Manu];
 515     Int_t bad_channel[800][1500];
 516     for ( Int_t i = 0; i < num_tot_BP ; i++ )
 517       { for ( Int_t j = 0; j < num_tot_Manu ; j++ )   bad_channel[i][j]=0;}
 518
 519
 520     char graphName[256];
 521
 522     // iterates over the first pedestal run
 523     TIter next(map[0]->CreateIterator());
 524     AliMUONVCalibParam* p;
 525
 526     Int_t    nmanu         = 0;
 527     Int_t    nGoodChannel   = 0;
 528     Int_t    nGoodChannel_a1   = 0;
 529     Int_t    nBadChannel   = 0;
 530     Int_t    nBadChannel_a1   = 0;
 531     Int_t    nBadChannel_a2   = 0;
 532     Int_t    nplot=0;
 533     Double_t sumProbChi2   = 0.;
 534     Double_t sumA1         = 0.;
 535     Double_t sumProbChi2P2 = 0.;
 536     Double_t sumA2         = 0.;
 537
 538     Double_t x[11], xErr[11], y[11], yErr[11];
 539     Double_t xp[11], xpErr[11], yp[11], ypErr[11];
 540
 541     while ( ( p = dynamic_cast<AliMUONVCalibParam*>(next() ) ) )
 542     {
 543       ped[0]  = p;
 544
 545       busPatchId = p->ID0();
 546       manuId     = p->ID1();
 547
 548       // read back pedestal from the other runs for the given (bupatch, manu)
 549       for (Int_t i = 1; i < nEntries; ++i) {
 550         ped[i] = static_cast<AliMUONVCalibParam*>(map[i]->FindObject(busPatchId, manuId));
 551       }
 552
 553       // compute for each channel the gain parameters
 554       for ( channelId = 0; channelId < ped[0]->Size() ; ++channelId ) {
 555
 556         gain=0.4;
 557
 558         Int_t n = 0;
 559         for (Int_t i = 0; i < nEntries; ++i) {
 560
 561           if (!ped[i]) continue; //shouldn't happen.
 562           pedMean[i]      = ped[i]->ValueAsDouble(channelId, 0);
 563           pedSigma[i]     = ped[i]->ValueAsDouble(channelId, 1);
 564           injCharge[i]    = (Double_t)run[i]->GetSecond();
 565           injChargeErr[i] = 0.01*injCharge[i];
 566           if(injChargeErr[i] <= 1.) injChargeErr[i]=1.;
 567
 568 //        if(n<2)cout << nEntries << " " << i << " " << injCharge[i] << endl;
 569
 570 //        cout << busPatchId << "\t" << manuId <<"\t"<< channelId << "\t" << n << " " << pedMean[i] << " " << pedSigma[i] << " " << injCharge[i] << " " << injChargeErr[i] << endl;
 571
 572           if (pedMean[i] < 0) continue; // not connected
 573
 574           if (pedSigma[i] <= 0) pedSigma[i] = 1.; // should not happen.
 575           n++;
 576         }
 577
 578
 579         // print_peak_mean_values
 580         if(gPrintLevel==3)
 581           {
 582
 583             fprintf(pfilep,"%4i%5i%5i%10.3f",busPatchId,manuId,channelId,0.);
 584             fprintf(pfilep,"%9.1f%9.1f%9.1f%9.1f%9.1f%9.1f%9.1f%9.1f%9.1f%9.1f%9.1f  mV\n",pedMean[0],pedMean[1],pedMean[2],pedMean[3],pedMean[4],pedMean[5],pedMean[6],pedMean[7],pedMean[8],pedMean[9],pedMean[10]);
 585             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",pedSigma[0],pedSigma[1],pedSigma[2],pedSigma[3],pedSigma[4],pedSigma[5],pedSigma[6],pedSigma[7],pedSigma[8],pedSigma[9],pedSigma[10]);
 586           }
 587
 588         // makegain
 589
 590
 591         // Fit Method:  Linear fit over 6 points + fit parabolic function  over 3  points)
 592
 593         // 1. - linear fit over 6 points
 594
 595         Double_t par[4] = {0.,0.,0.,gkADCMax};
 596
 597         Int_t nInit = 1;
 598         Int_t nbs   = nEntries - nInit;
 599         Int_t nbpf1 = 6; // linear fit over nbf1 points
 600
 601         for (Int_t j = 0; j < nbs; ++j)
 602         {
 603           Int_t k = j + nInit;
 604           x[j]    = pedMean[k];
 605           xErr[j] = pedSigma[k];
 606           y[j]    = injCharge[k];
 607           yErr[j] = injChargeErr[k];
 608
 609         }
 610
 611         TF1 *f1 = new TF1("f1",funcLin,0.,gkADCMax,2);
 612         TGraphErrors *graphErr = new TGraphErrors(nbpf1, x, y, xErr, yErr);
 613
 614         f1->SetParameters(0,0);
 615
 616         graphErr->Fit("f1","RQ");
 617
 618         chi2 = f1->GetChisquare();
 619         f1->GetParameters(par);
 620
 621         delete graphErr;
 622         graphErr=0;
 623         delete f1;
 624
 625         prChi2 = TMath::Prob(chi2, nbpf1 - 2);
 626
 627         Double_t xLim = pedMean[nInit + nbpf1 - 1];
 628         Double_t yLim = par[0]+par[1] * xLim;
 629
 630         a0 = par[0];
 631         a1 = par[1];
 632
 633         // 2. - Translation : new origin (xLim, yLim) + parabolic fit over nbf2 points
 634
 635 //      par[0] = 0;
 636         Int_t nbpf2 = nEntries - (nInit + nbpf1) + 1;
 637
 638         if(nbpf2 > 1)
 639         {
 640           for (Int_t j = 0; j < nbpf2; j++)
 641           {
 642             Int_t k  = j + (nInit + nbpf1) - 1;
 643             xp[j]    = pedMean[k] - xLim;
 644             xpErr[j] = pedSigma[k];
 645
 646             yp[j]    = injCharge[k] - yLim - par[1]*xp[j];
 647             ypErr[j] = injChargeErr[k];
 648
 649           }
 650
 651           TF1 *f2 = new TF1("f2",funcParabolic,0.,gkADCMax,1);
 652           TGraphErrors *graphErr = new TGraphErrors(nbpf2, xp, yp, xpErr, ypErr);
 653
 654           graphErr->Fit(f2,"RQ");
 655           chi2P2 = f2->GetChisquare();
 656           f2->GetParameters(par);
 657
 658           delete graphErr;
 659           graphErr=0;
 660           delete f2;
 661
 662           prChi2P2 = TMath::Prob(chi2P2, nbpf2-1);
 663
 664
 665               // ------------- print out in log file
 666 //        if (busPatchId == 6 && manuId == 116 && ( channelId >= 17 && channelId <= 20) )
 667 //          {
 668 //            filcout << " \n ********! Print_out.: BP= " << busPatchId << " Manu_Id= " << manuId
 669 //                      << " Ch.= " << channelId << ":" << endl;
 670
 671 //            for (Int_t j = 0; j < nbpf1; ++j)
 672 //              {filcout << j << " " << x[j] << " " << xErr[j] << " " << y[j] << " " << yErr[j] << endl;}
 673 //            filcout << "  a0,a1 = " << a0 << " , " << a1 << " pr_chi2 = " <<  prChi2 << endl ;
 674
 675 //            for (Int_t j = 0; j < nbpf2; ++j)
 676 //              {filcout << j << " " << xp[j] << " " << xpErr[j] << " " << yp[j] << " " << ypErr[j] << endl;}
 677 //            filcout << "  a2 = " << par[0] << " pr_chi2_2 = " <<  prChi2P2 << endl;
 678
 679 //          }
 680         // ------------------------------------------
 681
 682
 683
 684           a2 = par[0];
 685
 686           par[0] = a0;
 687           par[1] = a1;
 688           par[2] = a2;
 689           par[3] = xLim;
 690
 691 //        delete graphErr;
 692
 693         }
 694
 695         // Prints
 696
 697         p1 = TMath::Nint(ceil(prChi2*14))+1;
 698         p2 = TMath::Nint(ceil(prChi2P2*14))+1;
 699
 700         Double_t x0 = -par[0]/par[1]; // value of x corresponding to à 0 fC
 701         threshold = TMath::Nint(ceil(par[3]-x0)); // linear if x < threshold
 702
 703         if(gPrintLevel>=2)
 704         {
 705           fprintf(pfilen,"%4i %4i %2i",busPatchId,manuId,channelId);
 706           fprintf(pfilen," %6.2f %6.4f %10.3e %4.2f  %5.3f   %x  %5.3f  %x\n",
 707                   par[0], par[1], par[2], par[3], prChi2, p1, prChi2P2, p2);
 708         }
 709
 710         // some tests
 711
 712         if(par[1]< goodA1Min ||  par[1]> goodA1Max)
 713         {
 714           p1=0;
 715           nBadChannel_a1++;
 716           if (gPrintLevel && nBadChannel_a1 < 1)
 717           {
 718             cout << " !!!!! " << nBadChannel_a1 << " !!!!!!!! Bad Calib.: BP= " << busPatchId << " Manu_Id= " << manuId <<
 719                 " Ch.= " << channelId << ":";
 720             cout << "  a1 = " << par[1] << "    out of limit : [" <<  goodA1Min << "," << goodA1Max <<
 721                 "]" << endl;
 722           }
 723         }
 724
 725         if(par[2]< goodA2Min ||  par[2]> goodA2Max)
 726         {
 727           p2=0;
 728           nBadChannel_a2++;
 729           if (gPrintLevel && nBadChannel_a2 < 1)
 730           {
 731             cout << " !!!!! " << nBadChannel_a2 << " !!!!!!!! Bad Calib.: BP= " << busPatchId << " Manu_Id= " << manuId
 732                  << " Ch.= " << channelId << ":";
 733             cout << "  a2 = " << par[2] << "    out of limit : [" <<  goodA2Min << "," << goodA2Max
 734                  << "]" << endl;
 735
 736             for (Int_t j = 0; j < nbpf2; ++j)
 737               {cout << j << " " << x[j] << " " << xErr[j] << " " << y[j] << " " << yErr[j] << endl;}
 738           }
 739         }
 740
 741         Q  = p1*16 + p2;  // fit quality
 742         if(p1==0)Q=0;  // bad linear fit <=> bad calibration
 743
 744         if(p1>0 && p2>0)
 745           {
 746             nGoodChannel++;
 747             sumProbChi2P2   += prChi2P2;
 748             sumA2         += par[2];
 749           }
 750         else
 751           {
 752             nBadChannel++;
 753             if(busPatchId < num_tot_BP  && manuId < num_tot_Manu)  bad_channel[busPatchId][manuId]++;
 754             else{cout << " Warning : busPatch = " << busPatchId << " Manu = " << manuId << endl;}
 755             if(gPrintLevel>=2)fprintf(pfilef,"%4i %5i %2i %7.4f %10.3e %4i %2x\n",busPatchId,manuId,channelId,par[1],par[2],threshold,Q);
 756           }
 757
 758
 759         if(p1>0)
 760           {
 761             nGoodChannel_a1++;
 762             sumProbChi2   += prChi2;
 763             sumA1         += par[1];
 764             gain=1./(par[1]*capa);
 765           }
 766
 767
 768         tg->Fill();
 769
 770         if (!flatOutputFile.IsNull())
 771           {
 772             fprintf(pfilew,"%4i %5i %2i %7.4f %10.3e %4i %2x\n",busPatchId,manuId,channelId,par[1],par[2],threshold,Q);
 773           }
 774
 775         // Plots
 776
 777         if(gPlotLevel){
 778           if(Q==0  and  nplot < 100)
 779 //        if(p1>1 && p2==0  and  nplot < 100)
 780 //        if(p1>1 && p2>1  and  nplot < 100)
 781             {
 782               nplot++;
 783               TF1 *f2Calib = new TF1("f2Calib",funcCalib,0.,gkADCMax,NFITPARAMS);
 784
 785               TGraphErrors *graphErr = new TGraphErrors(nEntries,pedMean,injCharge,pedSigma,injChargeErr);
 786
 787               sprintf(graphName,"BusPatch_%d_Manu_%d_Ch_%d",busPatchId, manuId,channelId);
 788
 789               graphErr->SetTitle(graphName);
 790               graphErr->SetMarkerColor(3);
 791               graphErr->SetMarkerStyle(12);
 792               graphErr->Write(graphName);
 793
 794               sprintf(graphName,"f2_BusPatch_%d_Manu_%d_Ch_%d",busPatchId, manuId,channelId);
 795               f2Calib->SetTitle(graphName);
 796               f2Calib->SetLineColor(4);
 797               f2Calib->SetParameters(par);
 798               f2Calib->Write(graphName);
 799
 800               delete graphErr;
 801               delete f2Calib;
 802             }
 803         }
 804       }
 805       nmanu++;
 806       if(fmod(nmanu,100)==0)std::cout << " Nb manu = " << nmanu << std::endl;
 807     }
 808
 809     // file outputs for gain
 810 //     if (!flatOutputFile.IsNull())  fileout.close();
 811     if (!flatOutputFile.IsNull())  fclose(pfilew);
 812     if(gPrintLevel==2){ fclose(pfilen); fclose(pfilef);}
 813     if(gPrintLevel==3)  fclose(pfilep);
 814
 815     tg->Write();
 816     histoFile->Close();
 817
 818     //OutPut
 819     if (gPrintLevel)
 820     {
 821       filcout << "\n List of problematic BusPatch and Manu " << endl;
 822       filcout << " ========================================" << endl;
 823       filcout << "        BP       Manu        Nb Channel  " << endl ;
 824       filcout << " ========================================" << endl;
 825       for ( Int_t i = 0 ; i < num_tot_BP ; i++ )
 826         { for ( Int_t j = 0 ; j < num_tot_Manu ; j++ )
 827             if (bad_channel[i][j] != 0 ) filcout << "\t" << i << "\t " << j << "\t\t" << bad_channel[i][j] << endl;}
 828       filcout << " ========================================" << endl;
 829
 830
 831       filcout << "\n Nb of channels in raw data     = " << nmanu*64 << " (" << nmanu << " Manu)" <<  endl;
 832       filcout << "\n Nb of fully calibrated channel = " << nGoodChannel << " (" << goodA1Min << "<a1<" << goodA1Max
 833            << " and " << goodA2Min << "<a2<" << goodA2Max << ") " << endl;
 834       filcout << "\n Nb of Bad channel              = " << nBadChannel << endl;
 835
 836       filcout << "\n Nb of Good a1 channels  = " << nGoodChannel_a1 << " (" << goodA1Min << "<a1<" << goodA1Max <<  ") " << endl;
 837
 838       Double_t meanA1         = sumA1/(nGoodChannel_a1);
 839       Double_t meanProbChi2   = sumProbChi2/(nGoodChannel_a1);
 840       Double_t meanA2         = sumA2/(nGoodChannel);
 841       Double_t meanProbChi2P2 = sumProbChi2P2/(nGoodChannel);
 842
 843       Double_t capaManu = 0.2; // pF
 844       filcout << "\n linear fit   : <a1> = " << meanA1 << "\t  <gain>  = " <<  1./(meanA1*capaManu)
 845            << " mV/fC (capa= " << capaManu << " pF)" << endl;
 846       filcout <<   "        Prob(chi2)>  = " <<  meanProbChi2 << endl;
 847       filcout << "\n parabolic fit: <a2> = " << meanA2  << endl;
 848       filcout <<   "        Prob(chi2)>  = " <<  meanProbChi2P2 << "\n" << endl;
 849
 850     }
 851
 852
 853     return  ;
 854
 855 }
 856
 857 //*************************************************************//
 858
 859 // main routine
 860 int main(Int_t argc, Char_t **argv)
 861 {
 862
 863     // needed for streamer application
 864     gROOT->GetPluginManager()->AddHandler("TVirtualStreamerInfo",
 865                                           "*",
 866                                           "TStreamerInfo",
 867                                           "RIO",
 868                                           "TStreamerInfo()");
 869
 870     TFitter *minuitFit = new TFitter(NFITPARAMS);
 871     TVirtualFitter::SetFitter(minuitFit);
 872
 873     Int_t skipEvents = 0;
 874     Int_t maxEvents  = 1000000;
 875     Int_t MaxDateEvents  = 1000000;
 876     Int_t injCharge = 0;
 877     Double_t nSigma = 3;
 878     Int_t threshold = -1;
 879     Char_t inputFile[256];
 880
 881 //     TString flatOutputFile;
 882     TString fesOutputFile;
 883     TString crocusOutputFile;
 884     TString crocusConfigFile;
 885
 886 // option handler
 887
 888    // decode the input line
 889   for (Int_t i = 1; i < argc; i++) // argument 0 is the executable name
 890   {
 891       Char_t* arg;
 892
 893       arg = argv[i];
 894       if (arg[0] != '-') continue;
 895       switch (arg[1])
 896         {
 897         case 'f' :
 898           i++;
 899           sprintf(inputFile,argv[i]);
 900           break;
 901         case 'a' :
 902           i++;
 903           flatOutputFile = argv[i];
 904           break;
 905         case 'b' :
 906           i++;
 907           sprintf(gOutFolder,argv[i]);
 908           break;
 909         case 'o' :
 910           i++;
 911           crocusOutputFile = argv[i];
 912           break;
 913         case 'c' :
 914           i++;
 915           crocusConfigFile = argv[i];
 916           break;
 917 //      case 'v' :
 918 //        i++;
 919 //        logOutputFile = argv[i];
 920 //        break;
 921         case 'e' :
 922           i++;
 923           gCommand = argv[i];
 924           break;
 925         case 'd' :
 926           i++;
 927           gPrintLevel=atoi(argv[i]);
 928           break;
 929         case 'g' :
 930           i++;
 931           gPlotLevel=atoi(argv[i]);
 932           break;
 933         case 's' :
 934           i++;
 935           skipEvents=atoi(argv[i]);
 936           break;
 937         case 'l' :
 938           i++;
 939           injCharge=atoi(argv[i]);
 940           break;
 941         case 'm' :
 942           i++;
 943           sscanf(argv[i],"%d",&MaxDateEvents);
 944           break;
 945         case 'n' :
 946           i++;
 947           sscanf(argv[i],"%d",&maxEvents);
 948           break;
 949         case 'p' :
 950           i++;
 951           sscanf(argv[i],"%lf",&nSigma);
 952           break;
 953         case 'r' :
 954           i++;
 955           sprintf(gHistoFileName_gain,argv[i]);
 956           break;
 957         case 't' :
 958           i++;
 959           sscanf(argv[i],"%d",&threshold);
 960           break;
 961         case 'h' :
 962           i++;
 963           printf("\n******************* %s usage **********************",argv[0]);
 964           printf("\n%s -options, the available options are :",argv[0]);
 965           printf("\n-h help                   (this screen)");
 966           printf("\n");
 967           printf("\n Input");
 968           printf("\n-f <raw data file>        (default = %s)",inputFile);
 969           printf("\n-c <Crocus config. file>  (default = %s)",crocusConfigFile.Data());
 970           printf("\n");
 971           printf("\n Output");
 972           printf("\n-a <Flat ASCII file>      (default = %s)",flatOutputFile.Data());
 973           printf("\n-o <CROCUS cmd file>      (default = %s)",crocusOutputFile.Data());
 974 //        printf("\n-v <log output file>      (default = %s)",logOutputFile.Data());
 975           printf("\n");
 976           printf("\n Options");
 977           printf("\n-b <output directory>     (default = %s)",gOutFolder);
 978           printf("\n-d <print level>          (default = %d)",gPrintLevel);
 979           printf("\n-g <plot level>           (default = %d)",gPlotLevel);
 980           printf("\n-l <DAC level>            (default = %d)",injCharge);
 981           printf("\n-m <max date events>      (default = %d)",MaxDateEvents);
 982           printf("\n-s <skip events>          (default = %d)",skipEvents);
 983           printf("\n-n <max events>           (default = %d)",maxEvents);
 984           printf("\n-p <n sigmas>             (default = %f)",nSigma);
 985           printf("\n-r root file data for gain(default = %s)",gHistoFileName_gain);
 986           printf("\n-t <threshold (-1 = no)>  (default = %d)",threshold);
 987           printf("\n-e <execute ped/gain>     (default = %s)",gCommand.Data());
 988           printf("\n-e <gain create>           make gain & create a new root file");
 989           printf("\n-e <gain>                  make gain & update root file");
 990           printf("\n-e <gain compute>          make gain & compute gains");
 991
 992           printf("\n\n");
 993           exit(-1);
 994         default :
 995           printf("%s : bad argument %s (please check %s -h)\n",argv[0],argv[i],argv[0]);
 996           argc = 2; exit(-1); // exit if error
 997         } // end of switch
 998     } // end of for i
 999
1000   // set gCommand to lower case
1001   gCommand.ToLower();
1002
1003
1004   // decoding the events
1005
1006   Int_t status;
1007   void* event;
1008
1009   gPedMeanHisto = 0x0;
1010   gPedSigmaHisto = 0x0;
1011
1012   TStopwatch timers;
1013
1014   timers.Start(kTRUE);
1015
1016   // once we have a configuration file in db
1017   // copy locally a file from daq detector config db
1018   // The current detector is identified by detector code in variable
1019   // DATE_DETECTOR_CODE. It must be defined.
1020   // If environment variable DAQDA_TEST_DIR is defined, files are copied from DAQDA_TEST_DIR
1021   // instead of the database. The usual environment variables are not needed.
1022   if (!crocusConfigFile.IsNull()) {
1023     status = daqDA_DB_getFile("myconfig", crocusConfigFile.Data());
1024     if (status) {
1025       printf("Failed to get config file : %d\n",status);
1026       return -1;
1027     }
1028   }
1029
1030
1031   status = monitorSetDataSource(inputFile);
1032   if (status) {
1033     cerr << "ERROR : monitorSetDataSource status (hex) = " << hex << status
1034               << " " << monitorDecodeError(status) << endl;
1035     return -1;
1036   }
1037   status = monitorDeclareMp("MUON Tracking monitoring");
1038   if (status) {
1039     cerr << "ERROR : monitorDeclareMp status (hex) = " << hex << status
1040               << " " << monitorDecodeError(status) << endl;
1041     return -1;
1042   }
1043
1044   Int_t busPatchId;
1045   UShort_t manuId;
1046   UChar_t channelId;
1047   UShort_t charge;
1048   TString key("MUONTRKda :");
1049
1050   if (gCommand.CompareTo("comp") != 0)
1051     {
1052       cout << "\nMUONTRKda : Reading data from file " << inputFile <<endl;
1053
1054       while(1)
1055         {
1056           if (gNDateEvents >= MaxDateEvents) break;
1057           if (gNEvents >= maxEvents) break;
1058           if (gNEvents && gNEvents % 100 == 0)
1059             cout<<"Cumulated events " << gNEvents << endl;
1060
1061           // check shutdown condition
1062           if (daqDA_checkShutdown())
1063             break;
1064
1065           // Skip Events if needed
1066           while (skipEvents) {
1067             status = monitorGetEventDynamic(&event);
1068             skipEvents--;
1069           }
1070
1071           // starts reading
1072           status = monitorGetEventDynamic(&event);
1073           if (status < 0)  {
1074             cout<<"EOF found"<<endl;
1075             break;
1076           }
1077
1078 //        gNDateEvents++;
1079
1080           // decoding rawdata headers
1081           AliRawReader *rawReader = new AliRawReaderDate(event);
1082
1083           Int_t eventType = rawReader->GetType();
1084           gRunNumber = rawReader->GetRunNumber();
1085
1086           // Output log file initialisations
1087
1088           if(gNDateEvents==0)
1089             {
1090               if (gCommand.CompareTo("ped") == 0){
1091                 sprintf(flatFile,"%s/MUONTRKda_ped_%d.log",gOutFolder,gRunNumber);
1092                 logOutputFile=flatFile;
1093
1094                 filcout.open(logOutputFile.Data());
1095                 filcout<<"//=================================================" << endl;
1096                 filcout<<"//        MUONTRKda for Pedestal run = "   << gRunNumber << endl;
1097                 cout<<"\n ********  MUONTRKda for Pedestal run = " << gRunNumber << "\n" << endl;
1098               }
1099
1100               if (gCommand.Contains("gain")){
1101                 sprintf(flatFile,"%s/%s_%d_DAC_%d.log",gOutFolder,filenam,gRunNumber,injCharge);
1102                 logOutputFile=flatFile;
1103
1104                 filcout.open(logOutputFile.Data());
1105                 filcout<<"//=================================================" << endl;
1106                 filcout<<"//        MUONTRKda for Gain run = " << gRunNumber << "  (DAC=" << injCharge << ")" << endl;
1107                 cout<<"\n ********  MUONTRKda for Gain run = " << gRunNumber << "  (DAC=" << injCharge << ")\n" << endl;
1108               }
1109
1110               filcout<<"//=================================================" << endl;
1111               filcout<<"//   * Date          : " << date.AsString("l") << "\n" << endl;
1112               cout<<" * Date          : " << date.AsString("l") << "\n" << endl;
1113
1114             }
1115
1116           gNDateEvents++;
1117
1118
1119
1120           if (eventType != PHYSICS_EVENT)
1121             continue; // for the moment
1122
1123           // decoding MUON payload
1124           AliMUONRawStreamTracker* rawStream  = new AliMUONRawStreamTracker(rawReader);
1125           rawStream->DisableWarnings();
1126
1127           // loops over DDL
1128           rawStream->First();  // if GlitchError ? what we are doing ?
1129           while( (status = rawStream->Next(busPatchId, manuId, channelId, charge)) ) {
1130
1131             if (gNEvents == 0)
1132               gNChannel++;
1133
1134             //       if (gPrintLevel) printf("manuId: %d, channelId: %d charge: %d\n", manuId,
1135             //                       channelId, charge);
1136
1137             MakePed(busPatchId, (Int_t)manuId, (Int_t)channelId, (Int_t)charge);
1138
1139           } // Next digit
1140
1141           if (!rawStream->IsErrorMessage()) {
1142             gNEvents++;
1143           }
1144
1145           if (rawStream->GetPayLoad()->GetGlitchErrors())  gGlitchErrors++;
1146           if (rawStream->GetPayLoad()->GetParityErrors())  gParityErrors++;
1147           if (rawStream->GetPayLoad()->GetPaddingErrors()) gPaddingErrors++;
1148
1149           AliMUONLogger* log = rawStream->GetPayLoad()->GetErrorLogger();
1150           log->Print(key, filcout);
1151
1152           delete rawReader;
1153           delete rawStream;
1154
1155         } // while (1)
1156     }
1157
1158
1159
1160     if (gCommand.CompareTo("ped") == 0)
1161       {
1162         sprintf(flatFile,"MUONTRKda_ped_%d.ped",gRunNumber);
1163         if(flatOutputFile.IsNull())flatOutputFile=flatFile;
1164         MakePedStore(flatOutputFile);
1165       }
1166
1167   // option gain -> update root file with pedestal results
1168   // gain + create -> recreate root file
1169   // gain + comp -> update root file and compute gain parameters
1170
1171     if (gCommand.Contains("gain"))
1172       {
1173         MakePedStoreForGain(injCharge);
1174       }
1175
1176     if (gCommand.Contains("comp"))
1177       {
1178
1179 //      if(flatOutputFile.IsNull())flatOutputFile="MUONTRKda_gain.par";
1180 //      MakeGainStore(flatOutputFile);
1181         MakeGainStore();
1182       }
1183
1184
1185   delete gPedestalStore;
1186
1187   delete minuitFit;
1188   TVirtualFitter::SetFitter(0);
1189
1190   timers.Stop();
1191
1192   if (gCommand.CompareTo("comp") != 0)
1193     {
1194       cout << "MUONTRKda : Nb of DATE events     = "         << gNDateEvents    << endl;
1195       cout << "MUONTRKda : Nb of Glitch errors   = "         << gGlitchErrors   << endl;
1196       cout << "MUONTRKda : Nb of Parity errors   = "         << gParityErrors   << endl;
1197       cout << "MUONTRKda : Nb of Padding errors  = "         << gPaddingErrors  << endl;
1198       cout << "MUONTRKda : Nb of events used     = "         << gNEvents        << endl;
1199
1200       filcout << "MUONTRKda : Nb of DATE events     = "         << gNDateEvents    << endl;
1201       filcout << "MUONTRKda : Nb of Glitch errors   = "         << gGlitchErrors   << endl;
1202       filcout << "MUONTRKda : Nb of Parity errors   = "         << gParityErrors   << endl;
1203       filcout << "MUONTRKda : Nb of Padding errors  = "         << gPaddingErrors  << endl;
1204       filcout << "MUONTRKda : Nb of events used     = "         << gNEvents        << endl;
1205
1206     }
1207
1208
1209   cout << "\nMUONTRKda : Output file generated   : " << logOutputFile  << endl;
1210
1211   if (gCommand.CompareTo("ped") == 0)
1212     {
1213       if (!(crocusConfigFile.IsNull()))
1214         cout << "MUONTRKda : CROCUS command file generated : " << crocusOutputFile.Data() << endl;
1215       else
1216         cout << "MUONTRKda : WARNING no CROCUS command file generated" << endl;
1217       cout << "\nMUONTRKda : Histo file generated for pedestal : " << gHistoFileName  << endl;
1218     }
1219   else
1220     {
1221       cout << "\nMUONTRKda : Root data file for gain calculation    : " << gHistoFileName_gain  << endl;
1222     }
1223
1224   if (gCommand.CompareTo("comp") == 0)
1225     {
1226       cout << "MUONTRKda : Root Histo. file generated             : " << gRootFileName  << endl;
1227       cout << "MUONTRKda : Flat ASCII file generated (to SHUTTLE) : " << flatOutputFile << endl;
1228     }
1229
1230
1231
1232   // Store IN FES
1233
1234   if (gCommand.CompareTo("comp") == 0 || gCommand.CompareTo("ped") == 0)
1235     {
1236       printf("\n *****  STORE FILE in FES ****** \n");
1237
1238       // to be sure that env variable is set
1239       gSystem->Setenv("DAQDALIB_PATH", "$DATE_SITE/infoLogger");
1240
1241       if (!flatOutputFile.IsNull())
1242         {
1243
1244           //       flatOutputFile.Prepend("./");
1245
1246           status = daqDA_FES_storeFile(flatOutputFile.Data(),flatOutputFile.Data());
1247           if (status)
1248             {
1249               printf(" Failed to export file : %d\n",status);
1250             }
1251           else if(gPrintLevel) printf("Export file: %s\n",flatOutputFile.Data());
1252         }
1253     }
1254
1255   filcout.close();
1256   printf("\nExecution time : R:%7.2fs C:%7.2fs\n", timers.RealTime(), timers.CpuTime());
1257
1258   return status;
1259 }