MUON/MUONTRKda.cxx

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