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