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