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