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-10-08 New version: MUONTRKda.cxx,v 1.14
  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   FILE *pfilef = 0;
 559   if(gPrintLevel==2)
 560     {
 561       sprintf(filename,"%s/%s_%d.param",gOutFolder,filenam,gRunNumber);
 562       cout << " fit parameter file               = " << filename << "\n";
 563       pfilen = fopen (filename,"w");
 564
 565       fprintf(pfilen,"//===================================================================\n");
 566       fprintf(pfilen,"//  BP MANU CH. a0     a1       a2      xlim  P(chi2) p1 P(chi2)2  p2\n");
 567       fprintf(pfilen,"//===================================================================\n");
 568       fprintf(pfilen,"//   * Run           : %d \n",gRunNumber);
 569       fprintf(pfilen,"//===================================================================\n");
 570
 571       sprintf(filename,"%s/%s_%d.bad",gOutFolder,filenam,gRunNumber);
 572       cout << " Bad channel file                 = " << filename << "\n";
 573       pfilef = fopen (filename,"w");
 574
 575       fprintf(pfilef,"//=================================================\n");
 576       fprintf(pfilef,"//  Bad Channel file calculated by MUONTRKda \n");
 577       fprintf(pfilef,"//=================================================\n");
 578       fprintf(pfilef,"//   * Run           : %d \n",gRunNumber);
 579       fprintf(pfilef,"//   * Date          : %s \n",date.AsString("l"));
 580       fprintf(pfilef,"//=======================================\n");
 581       fprintf(pfilef,"// BP MANU CH.   a1      a2     thres. Q\n");
 582       fprintf(pfilef,"//=======================================\n");
 583     }
 584
 585   FILE *pfilew=0;
 586   if(flatOutputFile.IsNull())
 587     {
 588       sprintf(filename,"%s_%d.par",filenam,gRunNumber);
 589       flatOutputFile=filename;
 590     }
 591   if(!flatOutputFile.IsNull())
 592     {
 593       pfilew = fopen (flatOutputFile.Data(),"w");
 594
 595       fprintf(pfilew,"//================================================\n");
 596       fprintf(pfilew,"//  Calibration file calculated by MUONTRKda \n");
 597       fprintf(pfilew,"//=================================================\n");
 598       fprintf(pfilew,"//   * Run           : %d \n",gRunNumber);
 599       fprintf(pfilew,"//   * Date          : %s \n",date.AsString("l"));
 600       fprintf(pfilew,"//   * Statictics    : %d \n",gNEvents);
 601       fprintf(pfilew,"//   * # of MANUS    : %d \n",gNManu);
 602       fprintf(pfilew,"//   * # of channels : %d \n",gNChannel);
 603       fprintf(pfilew,"//-------------------------------------------------\n");
 604       if(nInit==0)
 605         fprintf(pfilew,"//   %d DAC values  fit:  %d pts (1st order) %d pts (2nd order) \n",nEntries,nbpf1,nbpf2);
 606       if(nInit==1)
 607         fprintf(pfilew,"//   %d DAC values  fit: %d pts (1st order) %d pts (2nd order) DAC=0 excluded\n",nEntries,nbpf1,nbpf2);
 608       fprintf(pfilew,"//   RUN     DAC   \n");
 609       fprintf(pfilew,"//-----------------\n");
 610       for (Int_t i = 0; i < nEntries; ++i) {
 611         tree->SetBranchAddress("run",&run[i]);
 612         fprintf(pfilew,"//   %d    %5.0f \n",num_RUN[i],injCharge[i]);
 613       }
 614       fprintf(pfilew,"//=======================================\n");
 615       fprintf(pfilew,"// BP MANU CH.   a1      a2     thres. Q\n");
 616       fprintf(pfilew,"//=======================================\n");
 617     }
 618
 619   FILE *pfilep = 0;
 620   if(gPrintLevel==2)
 621     {
 622       sprintf(filename,"%s/%s_%d.peak",gOutFolder,filenam,gRunNumber);
 623       cout << " File containing Peak mean values = " << filename << "\n";
 624       pfilep = fopen (filename,"w");
 625
 626       fprintf(pfilep,"//==============================================================================================================================\n");
 627       fprintf(pfilep,"//   * Run           : %d \n",gRunNumber);
 628       fprintf(pfilep,"//==============================================================================================================================\n");
 629       fprintf(pfilep,"// BP  MANU  CH.    Ped.     <0>      <1>      <2>      <3>      <4>      <5>      <6>      <7>      <8>      <9>     <10> \n");
 630       fprintf(pfilep,"//==============================================================================================================================\n");
 631       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]);
 632       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]);
 633       fprintf(pfilep,"//==============================================================================================================================\n");
 634     }
 635
 636
 637
 638   //  plot out
 639
 640   TFile* gainFile = 0x0;
 641   sprintf(gRootFileName,"%s/%s_%d.root",gOutFolder,filenam,gRunNumber);
 642   gainFile = new TFile(gRootFileName,"RECREATE");
 643
 644   Double_t chi2    = 0.;
 645   Double_t chi2P2  = 0.;
 646   Double_t prChi2  = 0;
 647   Double_t prChi2P2 =0;
 648   Double_t a0=0.,a1=1.,a2=0.;
 649   Int_t busPatchId ;
 650   Int_t manuId     ;
 651   Int_t channelId ;
 652   Int_t threshold = 0;
 653   Int_t Q = 0;
 654   Int_t p1 =0;
 655   Int_t p2 =0;
 656   Double_t gain=0;
 657   Double_t capa=0.2; // internal capacitor (pF)
 658
 659   TTree *tg = new TTree("tg","TTree avec class Manu_DiMu");
 660
 661   tg->Branch("bp",&busPatchId, "busPatchId/I");
 662   tg->Branch("manu",&manuId, "manuId/I");
 663   tg->Branch("channel",&channelId, "channelId/I");
 664
 665   tg->Branch("a0",&a0, "a0/D");
 666   tg->Branch("a1",&a1, "a1/D");
 667   tg->Branch("a2",&a2, "a2/D");
 668   tg->Branch("Pchi2",&prChi2, "prChi2/D");
 669   tg->Branch("Pchi2_2",&prChi2P2, "prChi2P2/D");
 670   tg->Branch("Threshold",&threshold, "threshold/I");
 671   tg->Branch("Q",&Q, "Q/I");
 672   tg->Branch("p1",&p1, "p1/I");
 673   tg->Branch("p2",&p2, "p2/I");
 674   tg->Branch("gain",&gain, "gain/D");
 675
 676   char graphName[256];
 677
 678   // iterates over the first pedestal run
 679   TIter next(map[0]->CreateIterator());
 680   AliMUONVCalibParam* p;
 681
 682   Int_t    nmanu         = 0;
 683   Int_t    nGoodChannel   = 0;
 684   Int_t    nBadChannel   = 0;
 685   Int_t    nBadChannel_a1   = 0;
 686   Int_t    nBadChannel_a2   = 0;
 687   Int_t    noFitChannel   = 0;
 688   Int_t    nplot=0;
 689   Double_t sumProbChi2   = 0.;
 690   Double_t sumA1         = 0.;
 691   Double_t sumProbChi2P2 = 0.;
 692   Double_t sumA2         = 0.;
 693
 694   Double_t x[11], xErr[11], y[11], yErr[11];
 695   Double_t xp[11], xpErr[11], yp[11], ypErr[11];
 696
 697   Int_t uncalcountertotal=0 ;
 698
 699   while ( ( p = dynamic_cast<AliMUONVCalibParam*>(next() ) ) )
 700     {
 701       ped[0]  = p;
 702
 703       busPatchId = p->ID0();
 704       manuId     = p->ID1();
 705
 706       // read back pedestal from the other runs for the given (bupatch, manu)
 707       for (Int_t i = 1; i < nEntries; ++i) {
 708         ped[i] = static_cast<AliMUONVCalibParam*>(map[i]->FindObject(busPatchId, manuId));
 709       }
 710
 711       // compute for each channel the gain parameters
 712       for ( channelId = 0; channelId < ped[0]->Size() ; ++channelId )
 713         {
 714
 715           Int_t n = 0;
 716           for (Int_t i = 0; i < nEntries; ++i) {
 717
 718             if (!ped[i]) continue; //shouldn't happen.
 719             pedMean[i]      = ped[i]->ValueAsDouble(channelId, 0);
 720             pedSigma[i]     = ped[i]->ValueAsDouble(channelId, 1);
 721
 722             if (pedMean[i] < 0) continue; // not connected
 723
 724             if (pedSigma[i] <= 0) pedSigma[i] = 1.; // should not happen.
 725             n++;
 726           }
 727
 728
 729           // print_peak_mean_values
 730           if(gPrintLevel==2)
 731             {
 732
 733               fprintf(pfilep,"%4i%5i%5i%10.3f",busPatchId,manuId,channelId,0.);
 734               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]);
 735               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]);
 736             }
 737
 738           // makegain
 739
 740
 741           // Fit Method:  Linear fit over nbpf1 points + parabolic fit  over nbpf2  points)
 742           // nInit=1 : 1st pt DAC=0 excluded
 743
 744           // 1. - linear fit over nbpf1 points
 745
 746           Double_t par[4] = {0.,0.5,0.,gkADCMax};
 747           Int_t nbs   = nEntries - nInit;
 748           if(nbs < nbpf1)nbpf1=nbs;
 749
 750           Int_t FitProceed=1;
 751           for (Int_t j = 0; j < nbs; ++j)
 752             {
 753               Int_t k = j + nInit;
 754               x[j]    = pedMean[k];
 755               if(x[j]==0.)FitProceed=0;
 756               xErr[j] = pedSigma[k];
 757               y[j]    = injCharge[k];
 758               yErr[j] = injChargeErr[k];
 759
 760             }
 761
 762           TGraphErrors *graphErr;
 763           if(!FitProceed) { p1=0; p2=0; noFitChannel++;}
 764
 765           if(FitProceed)
 766             {
 767
 768               TF1 *f1 = new TF1("f1",funcLin,0.,gkADCMax,2);
 769               graphErr = new TGraphErrors(nbpf1, x, y, xErr, yErr);
 770
 771               f1->SetParameters(0,0);
 772
 773               graphErr->Fit("f1","RQ");
 774
 775               chi2 = f1->GetChisquare();
 776               f1->GetParameters(par);
 777
 778               delete graphErr;
 779               graphErr=0;
 780               delete f1;
 781
 782               prChi2 = TMath::Prob(chi2, nbpf1 - 2);
 783
 784               Double_t xLim = pedMean[nInit + nbpf1 - 1];
 785               Double_t yLim = par[0]+par[1] * xLim;
 786
 787               a0 = par[0];
 788               a1 = par[1];
 789
 790               // 2. - Translation : new origin (xLim, yLim) + parabolic fit over nbf2 points
 791
 792               if(nbpf2 > 1)
 793                 {
 794                   for (Int_t j = 0; j < nbpf2; j++)
 795                     {
 796                       Int_t k  = j + (nInit + nbpf1) - 1;
 797                       xp[j]    = pedMean[k] - xLim;
 798                       xpErr[j] = pedSigma[k];
 799
 800                       yp[j]    = injCharge[k] - yLim - par[1]*xp[j];
 801                       ypErr[j] = injChargeErr[k];
 802
 803                     }
 804
 805                   TF1 *f2 = new TF1("f2",funcParabolic,0.,gkADCMax,1);
 806                   graphErr = new TGraphErrors(nbpf2, xp, yp, xpErr, ypErr);
 807
 808                   graphErr->Fit(f2,"RQ");
 809                   chi2P2 = f2->GetChisquare();
 810                   f2->GetParameters(par);
 811
 812                   delete graphErr;
 813                   graphErr=0;
 814                   delete f2;
 815
 816                   prChi2P2 = TMath::Prob(chi2P2, nbpf2-1);
 817                   a2 = par[0];
 818
 819                   //      delete graphErr;
 820
 821                 }
 822
 823               par[0] = a0;
 824               par[1] = a1;
 825               par[2] = a2;
 826               par[3] = xLim;
 827
 828
 829               // Prints
 830
 831               p1 = TMath::Nint(ceil(prChi2*14))+1;
 832               p2 = TMath::Nint(ceil(prChi2P2*14))+1;
 833               Q  = p1*16 + p2;  // fit quality
 834
 835               Double_t x0 = -par[0]/par[1]; // value of x corresponding to à 0 fC
 836               threshold = TMath::Nint(ceil(par[3]-x0)); // linear if x < threshold
 837
 838               if(gPrintLevel==2)
 839                 {
 840                   fprintf(pfilen,"%4i %4i %2i",busPatchId,manuId,channelId);
 841                   fprintf(pfilen," %6.2f %6.4f %10.3e %4.2f  %5.3f   %x  %5.3f  %x\n",
 842                           par[0], par[1], par[2], par[3], prChi2, p1, prChi2P2, p2);
 843                 }
 844
 845               // some tests
 846
 847               if(par[1]< goodA1Min ||  par[1]> goodA1Max)
 848                 {
 849                   p1=0;
 850                   nBadChannel_a1++;
 851                   if (gPrintLevel && nBadChannel_a1 < 1)
 852                     {
 853                       cout << " !!!!! " << nBadChannel_a1 << " !!!!!!!! Bad Calib.: BP= " << busPatchId << " Manu_Id= " << manuId <<
 854                         " Ch.= " << channelId << ":";
 855                       cout << "  a1 = " << par[1] << "    out of limit : [" <<  goodA1Min << "," << goodA1Max <<
 856                         "]" << endl;
 857                     }
 858                 }
 859
 860
 861               if(par[2]< goodA2Min ||  par[2]> goodA2Max)
 862                 {
 863                   p2=0;
 864                   nBadChannel_a2++;
 865                   if (gPrintLevel && nBadChannel_a2 < 1)
 866                     {
 867                       cout << " !!!!! " << nBadChannel_a2 << " !!!!!!!! Bad Calib.: BP= " << busPatchId << " Manu_Id= " << manuId
 868                            << " Ch.= " << channelId << ":";
 869                       cout << "  a2 = " << par[2] << "    out of limit : [" <<  goodA2Min << "," << goodA2Max
 870                            << "]" << endl;
 871
 872                       for (Int_t j = 0; j < nbpf2; ++j)
 873                         {cout << j << " " << x[j] << " " << xErr[j] << " " << y[j] << " " << yErr[j] << endl;}
 874                     }
 875                 }
 876             } // FitProceed
 877
 878           if(FitProceed && p1>0 && p2>0)
 879             {
 880               nGoodChannel++;
 881               sumProbChi2   += prChi2;
 882               sumA1         += par[1];
 883               gain=1./(par[1]*capa);
 884               sumProbChi2P2   += prChi2P2;
 885               sumA2         += par[2];
 886             }
 887           else // bad calibration
 888             {
 889               nBadChannel++;
 890               Q=0;
 891               if(gPrintLevel==2)fprintf(pfilef,"%4i %5i %2i %7.4f %10.3e %4i %2x\n",busPatchId,manuId,channelId,par[1],par[2],threshold,Q);
 892               par[1]=0.5; a1=0.5; p1=0;
 893               par[2]=0.;        a2=0.; p2=0;
 894               threshold=gkADCMax;
 895
 896               char bpmanuname[256];
 897               ErrorCounter* uncalcounter;
 898
 899               sprintf(bpmanuname,"bp%dmanu%d",busPatchId,manuId);
 900               if (!(uncalcounter = (ErrorCounter*)gUncalBuspatchManuTable->FindObject(bpmanuname)))
 901                 {
 902                   // New buspatch_manu name
 903                   uncalcounter= new ErrorCounter (busPatchId,manuId);
 904                   uncalcounter->SetName(bpmanuname);
 905                   gUncalBuspatchManuTable->Add(uncalcounter);
 906                 }
 907               else
 908                 {
 909                   // Existing buspatch_manu name
 910                   uncalcounter->Increment();
 911                 }
 912               //                            uncalcounter->Print_uncal()
 913               uncalcountertotal ++;
 914             }
 915
 916           if(gPlotLevel){
 917             //                if(Q==0  and  nplot < 100)
 918             //    if(p1>1 && p2==0  and  nplot < 100)
 919             if(p1>1 && p2>1  and  nplot < 100)
 920               //        if(p1>=1 and p1<=2  and  nplot < 100)
 921               {
 922                 nplot++;
 923                 //            cout << " nplot = " << nplot << endl;
 924                 TF1 *f2Calib = new TF1("f2Calib",funcCalib,0.,gkADCMax,NFITPARAMS);
 925
 926                 graphErr = new TGraphErrors(nEntries,pedMean,injCharge,pedSigma,injChargeErr);
 927
 928                 sprintf(graphName,"BusPatch_%d_Manu_%d_Ch_%d",busPatchId, manuId,channelId);
 929
 930                 graphErr->SetTitle(graphName);
 931                 graphErr->SetMarkerColor(3);
 932                 graphErr->SetMarkerStyle(12);
 933                 graphErr->Write(graphName);
 934
 935                 sprintf(graphName,"f2_BusPatch_%d_Manu_%d_Ch_%d",busPatchId, manuId,channelId);
 936                 f2Calib->SetTitle(graphName);
 937                 f2Calib->SetLineColor(4);
 938                 f2Calib->SetParameters(par);
 939                 f2Calib->Write(graphName);
 940
 941                 delete graphErr;
 942                 graphErr=0;
 943                 delete f2Calib;
 944               }
 945           }
 946
 947
 948           tg->Fill();
 949
 950           if (!flatOutputFile.IsNull())
 951             {
 952               fprintf(pfilew,"%4i %5i %2i %7.4f %10.3e %4i %2x\n",busPatchId,manuId,channelId,par[1],par[2],threshold,Q);
 953             }
 954
 955         }
 956       nmanu++;
 957       if(fmod(nmanu,100)==0)std::cout << " Nb manu = " << nmanu << std::endl;
 958     }
 959
 960   // file outputs for gain
 961   if (!flatOutputFile.IsNull())  fclose(pfilew);
 962   if(gPrintLevel==2){ fclose(pfilen); fclose(pfilef); fclose(pfilep); }
 963
 964   tg->Write();
 965   histoFile->Close();
 966
 967   //OutPut
 968   if (gPrintLevel)
 969     {
 970       // print in logfile
 971       if (gUncalBuspatchManuTable->GetSize())
 972         {
 973           gUncalBuspatchManuTable->Sort();  // use compare
 974           TIterator* iter = gUncalBuspatchManuTable->MakeIterator();
 975           ErrorCounter* uncalcounter;
 976           filcouc << "\n List of problematic BusPatch and Manu " << endl;
 977           filcouc << " ========================================" << endl;
 978           filcouc << "        BP       Manu        Nb Channel  " << endl ;
 979           filcouc << " ========================================" << endl;
 980           while((uncalcounter = (ErrorCounter*) iter->Next()))
 981             {
 982               filcouc << "\t" << uncalcounter->BusPatch() << "\t " << uncalcounter->ManuId() << "\t\t"   << uncalcounter->Events() << endl;
 983             }
 984           filcouc << " ========================================" << endl;
 985
 986           filcouc << " Number of bad calibrated Manu    = " << gUncalBuspatchManuTable->GetSize() << endl ;
 987           filcouc << " Number of bad calibrated channel = " << uncalcountertotal << endl;
 988
 989         }
 990
 991
 992       filcouc << "\n Nb of channels in raw data = " << nmanu*64 << " (" << nmanu << " Manu)" <<  endl;
 993       filcouc << " Nb of calibrated channel   = " << nGoodChannel << " (" << goodA1Min << "<a1<" << goodA1Max
 994               << " and " << goodA2Min << "<a2<" << goodA2Max << ") " << endl;
 995       filcouc << " Nb of uncalibrated channel = " << nBadChannel << " (" << noFitChannel << " unfitted)" << endl;
 996
 997       cout << "\n Nb of channels in raw data = " << nmanu*64 << " (" << nmanu << " Manu)" <<  endl;
 998       cout << " Nb of calibrated channel   = " << nGoodChannel << " (" << goodA1Min << "<a1<" << goodA1Max
 999            << " and " << goodA2Min << "<a2<" << goodA2Max << ") " << endl;
1000       cout << " Nb of uncalibrated channel = " << nBadChannel << " (" << noFitChannel << " unfitted)" << endl;
1001
1002       Double_t meanA1         = sumA1/(nGoodChannel);
1003       Double_t meanProbChi2   = sumProbChi2/(nGoodChannel);
1004       Double_t meanA2         = sumA2/(nGoodChannel);
1005       Double_t meanProbChi2P2 = sumProbChi2P2/(nGoodChannel);
1006
1007       Double_t capaManu = 0.2; // pF
1008       filcouc << "\n linear fit   : <a1> = " << meanA1 << "\t  <gain>  = " <<  1./(meanA1*capaManu)
1009               << " mV/fC (capa= " << capaManu << " pF)" << endl;
1010       filcouc <<   "        Prob(chi2)>  = " <<  meanProbChi2 << endl;
1011       filcouc << "\n parabolic fit: <a2> = " << meanA2  << endl;
1012       filcouc <<   "        Prob(chi2)>  = " <<  meanProbChi2P2 << "\n" << endl;
1013
1014       cout << "\n  <gain>  = " <<  1./(meanA1*capaManu)
1015            << " mV/fC (capa= " << capaManu << " pF)"
1016            <<  "  Prob(chi2)>  = " <<  meanProbChi2 << endl;
1017     }
1018
1019   filcouc.close();
1020
1021   return  ;
1022
1023 }
1024
1025 //*************************************************************//
1026
1027 // main routine
1028 int main(Int_t argc, Char_t **argv)
1029 {
1030
1031   // needed for streamer application
1032   gROOT->GetPluginManager()->AddHandler("TVirtualStreamerInfo",
1033                                         "*",
1034                                         "TStreamerInfo",
1035                                         "RIO",
1036                                         "TStreamerInfo()");
1037
1038   TFitter *minuitFit = new TFitter(NFITPARAMS);
1039   TVirtualFitter::SetFitter(minuitFit);
1040
1041   //    ofstream filcout;
1042
1043   Int_t skipEvents = 0;
1044   Int_t maxEvents  = 1000000;
1045   Int_t MaxDateEvents  = 1000000;
1046   Int_t injCharge = 0;
1047   Char_t inputFile[256]="";
1048   Char_t inputFile1[256]="";
1049   Char_t inputFile2[256]="";
1050
1051   Int_t gGlitchErrors= 0;
1052   Int_t gParityErrors= 0;
1053   Int_t gPaddingErrors= 0;
1054   Int_t recoverParityErrors = 1;
1055
1056   TString fesOutputFile;
1057
1058   // option handler
1059
1060   // decode the input line
1061   for (Int_t i = 1; i < argc; i++) // argument 0 is the executable name
1062     {
1063       Char_t* arg;
1064
1065       arg = argv[i];
1066       if (arg[0] != '-') continue;
1067       switch (arg[1])
1068         {
1069         case 'f' :
1070           i++;
1071           sprintf(inputFile1,argv[i]);
1072           break;
1073         case 'z' :
1074           i++;
1075           sprintf(inputFile2,argv[i]);
1076           break;
1077         case 'a' :
1078           i++;
1079           flatOutputFile = argv[i];
1080           break;
1081         case 'b' :
1082           i++;
1083           sprintf(gOutFolder,argv[i]);
1084           break;
1085         case 'c' :
1086           i++;
1087           gFES=atoi(argv[i]);
1088           break;
1089         case 'd' :
1090           i++;
1091           gPrintLevel=atoi(argv[i]);
1092           break;
1093         case 'e' :
1094           i++;
1095           gCommand = argv[i];
1096           break;
1097         case 'g' :
1098           i++;
1099           gPlotLevel=atoi(argv[i]);
1100           break;
1101         case 'i' :
1102           i++;
1103           nbpf1=atoi(argv[i]);
1104           break;
1105         case 's' :
1106           i++;
1107           skipEvents=atoi(argv[i]);
1108           break;
1109         case 'l' :
1110           i++;
1111           injCharge=atoi(argv[i]);
1112           break;
1113         case 'm' :
1114           i++;
1115           sscanf(argv[i],"%d",&MaxDateEvents);
1116           break;
1117         case 'n' :
1118           i++;
1119           sscanf(argv[i],"%d",&maxEvents);
1120           break;
1121         case 'r' :
1122           i++;
1123           sprintf(gHistoFileName_gain,argv[i]);
1124           break;
1125         case 'p' :
1126           i++;
1127           sscanf(argv[i],"%d",&recoverParityErrors);
1128           break;
1129         case 'h' :
1130           i++;
1131           printf("\n******************* %s usage **********************",argv[0]);
1132           printf("\n%s -options, the available options are :",argv[0]);
1133           printf("\n-h help                    (this screen)");
1134           printf("\n");
1135           printf("\n Input");
1136           printf("\n-f <raw data file>         (default = %s)",inputFile1);
1137           printf("\n-z <raw data file2>        (default = %s)",inputFile2);
1138           printf("\n");
1139           printf("\n Output");
1140           printf("\n-a <Flat ASCII file>       (default = %s)",flatOutputFile.Data());
1141           printf("\n");
1142           printf("\n Options");
1143           printf("\n-b <output directory>      (default = %s)",gOutFolder);
1144           printf("\n-c <FES switch>            (default = %d)",gFES);
1145           printf("\n-d <print level>           (default = %d)",gPrintLevel);
1146           printf("\n-g <plot level>            (default = %d)",gPlotLevel);
1147           printf("\n-i <nb linear points>      (default = %d)",nbpf1);
1148           printf("\n-l <DAC level>             (default = %d)",injCharge);
1149           printf("\n-m <max date events>       (default = %d)",MaxDateEvents);
1150           printf("\n-s <skip events>           (default = %d)",skipEvents);
1151           printf("\n-n <max events>            (default = %d)",maxEvents);
1152           printf("\n-r root file data for gain (default = %s)",gHistoFileName_gain);
1153           printf("\n-e <execute ped/gain>      (default = %s)",gCommand.Data());
1154           printf("\n-e <gain create>           make gain & create a new root file");
1155           printf("\n-e <gain>                  make gain & update root file");
1156           printf("\n-e <gain compute>          make gain & compute gains");
1157           printf("\n-p <Recover parity errors> (default = %d)",recoverParityErrors);
1158
1159           printf("\n\n");
1160           exit(-1);
1161         default :
1162           printf("%s : bad argument %s (please check %s -h)\n",argv[0],argv[i],argv[0]);
1163           argc = 2; exit(-1); // exit if error
1164         } // end of switch
1165     } // end of for i
1166
1167   // set gCommand to lower case
1168   gCommand.ToLower();
1169
1170
1171   // decoding the events
1172
1173   Int_t status=0;
1174   void* event;
1175
1176   gPedMeanHisto = 0x0;
1177   gPedSigmaHisto = 0x0;
1178
1179   TStopwatch timers;
1180
1181   timers.Start(kTRUE);
1182
1183   UShort_t manuId;
1184   UChar_t channelId;
1185   UShort_t charge;
1186   TString key("MUONTRKda :");
1187
1188   AliMUONRawStreamTrackerHP* rawStream  = 0;
1189
1190   int nrawdata=1;
1191   if (strlen(inputFile2) > 0) nrawdata=2;  //    2 raw data files to read (offline case)
1192
1193   if (gCommand.CompareTo("comp") != 0)
1194     {
1195
1196       for(int iraw = 0; iraw < nrawdata; ++iraw)
1197         {
1198           if(iraw==0) strcpy(inputFile,inputFile1);
1199           if(iraw==1) strcpy(inputFile,inputFile2);
1200
1201           status = monitorSetDataSource(inputFile);
1202
1203           if (status) {
1204             cerr << "ERROR : monitorSetDataSource status (hex) = " << hex << status
1205                  << " " << monitorDecodeError(status) << endl;
1206             return -1;
1207           }
1208           status = monitorDeclareMp("MUON Tracking monitoring");
1209           if (status) {
1210             cerr << "ERROR : monitorDeclareMp status (hex) = " << hex << status
1211                  << " " << monitorDecodeError(status) << endl;
1212             return -1;
1213           }
1214
1215
1216           cout << "\nMUONTRKda : Reading data from file " << inputFile << " gNEvents = " << gNEvents << endl;
1217
1218           while(1)
1219             {
1220               if (gNDateEvents >= MaxDateEvents) break;
1221               if (gNEvents >= maxEvents) break;
1222               if (gNDateEvents>0 &&  gNDateEvents % 100 == 0)
1223                 cout<<"Cumulated:  DATE events = " << gNDateEvents << "   Used events = " << gNEvents << endl;
1224
1225               // check shutdown condition
1226               if (daqDA_checkShutdown())
1227                 break;
1228
1229               // Skip Events if needed
1230               while (skipEvents) {
1231                 status = monitorGetEventDynamic(&event);
1232                 skipEvents--;
1233               }
1234
1235               // starts reading
1236               status = monitorGetEventDynamic(&event);
1237               if (status < 0)  {
1238                 cout<<"EOF found"<<endl;
1239                 break;
1240               }
1241
1242               // decoding rawdata headers
1243               AliRawReader *rawReader = new AliRawReaderDate(event);
1244
1245               Int_t eventType = rawReader->GetType();
1246               gRunNumber = rawReader->GetRunNumber();
1247
1248               // Output log file initialisations
1249
1250               if(gNDateEvents==0)
1251                 {
1252                   if (gCommand.CompareTo("ped") == 0){
1253                     sprintf(flatFile,"%s/MUONTRKda_ped_%d.log",gOutFolder,gRunNumber);
1254                     logOutputFile=flatFile;
1255
1256                     filcout.open(logOutputFile.Data());
1257                     filcout<<"//=================================================" << endl;
1258                     filcout<<"//        MUONTRKda for Pedestal run = "   << gRunNumber << endl;
1259                     cout<<"\n ********  MUONTRKda for Pedestal run = " << gRunNumber << "\n" << endl;
1260                   }
1261
1262                   if (gCommand.Contains("gain")){
1263                     sprintf(flatFile,"%s/%s_%d_DAC_%d.log",gOutFolder,filenam,gRunNumber,injCharge);
1264                     logOutputFile=flatFile;
1265
1266                     filcout.open(logOutputFile.Data());
1267                     filcout<<"//=================================================" << endl;
1268                     filcout<<"//        MUONTRKda for Gain run = " << gRunNumber << "  (DAC=" << injCharge << ")" << endl;
1269                     cout<<"\n ********  MUONTRKda for Gain run = " << gRunNumber << "  (DAC=" << injCharge << ")\n" << endl;
1270                   }
1271
1272                   filcout<<"//=================================================" << endl;
1273                   filcout<<"//   * Date          : " << date.AsString("l") << "\n" << endl;
1274                   cout<<" * Date          : " << date.AsString("l") << "\n" << endl;
1275
1276                 }
1277
1278               gNDateEvents++;
1279
1280               if (eventType != PHYSICS_EVENT)
1281                 continue; // for the moment
1282
1283               // decoding MUON payload
1284               rawStream  = new AliMUONRawStreamTrackerHP(rawReader);
1285               rawStream->DisableWarnings();
1286               rawStream->EnabbleErrorLogger();
1287
1288               // First lopp over DDL's to find good events
1289               // Error counters per event (counters in the decoding lib are for each DDL)
1290               Bool_t eventIsErrorMessage = kFALSE;
1291               int eventGlitchErrors = 0;
1292               int eventParityErrors = 0;
1293               int eventPaddingErrors = 0;
1294               rawStream->First();
1295               do
1296                 {
1297                   if (rawStream->IsErrorMessage()) eventIsErrorMessage = kTRUE;
1298                   eventGlitchErrors += rawStream->GetGlitchErrors();
1299                   eventParityErrors += rawStream->GetParityErrors();
1300                   eventPaddingErrors += rawStream->GetPaddingErrors();
1301                 } while(rawStream->NextDDL());
1302
1303               AliMUONRawStreamTrackerHP::AliBusPatch* busPatch;
1304               if (!eventIsErrorMessage)
1305                 {
1306                   // Good events (no error) -> compute pedestal for all channels
1307                   rawStream->First();
1308                   while( (busPatch = (AliMUONRawStreamTrackerHP::AliBusPatch*) rawStream->Next()))
1309                     {
1310                       for(int i = 0; i < busPatch->GetLength(); ++i)
1311                         {
1312                           if (gNEvents == 0) gNChannel++;
1313                           busPatch->GetData(i, manuId, channelId, charge);
1314                           MakePed(busPatch->GetBusPatchId(), (Int_t)manuId, (Int_t)channelId, (Int_t)charge);
1315                         }
1316                     }
1317                   gNEvents++;
1318                 }
1319               else
1320                 {
1321                   // Events with errors
1322                   if (recoverParityErrors && eventParityErrors && !eventGlitchErrors&& !eventPaddingErrors)
1323                     {
1324                       // Recover parity errors -> compute pedestal for all good buspatches
1325                       if ( TEST_SYSTEM_ATTRIBUTE( rawReader->GetAttributes(),
1326                                                   ATTR_ORBIT_BC ))
1327                         {
1328                           filcout <<"Event recovered -> Period:"<<EVENT_ID_GET_PERIOD( rawReader->GetEventId() )
1329                                   <<" Orbit:"<<EVENT_ID_GET_ORBIT( rawReader->GetEventId() )
1330                                   <<" BunchCrossing:"<<EVENT_ID_GET_BUNCH_CROSSING( rawReader->GetEventId() )<<endl;
1331                         }
1332                       else
1333                         {
1334                           filcout <<"Event recovered -> nbInRun:"<<EVENT_ID_GET_NB_IN_RUN( rawReader->GetEventId() )
1335                                   <<" burstNb:"<<EVENT_ID_GET_BURST_NB( rawReader->GetEventId() )
1336                                   <<" nbInBurst:"<<EVENT_ID_GET_NB_IN_BURST( rawReader->GetEventId() )<<endl;
1337                         }
1338                       rawStream->First();
1339                       while( (busPatch = (AliMUONRawStreamTrackerHP::AliBusPatch*) rawStream->Next()))
1340                         {
1341                           // Check the buspatch -> if error not use it in the pedestal calculation
1342                           int errorCount = 0;
1343                           for(int i = 0; i < busPatch->GetLength(); ++i)
1344                             {
1345                               if (!busPatch->IsParityOk(i)) errorCount++;
1346                             }
1347                           if (!errorCount)
1348                             {
1349                               // Good buspatch
1350                               for(int i = 0; i < busPatch->GetLength(); ++i)
1351                                 {
1352                                   if (gNEvents == 0) gNChannel++;
1353                                   busPatch->GetData(i, manuId, channelId, charge);
1354                                   // if (busPatch->GetBusPatchId()==1719 && manuId == 1 && channelId == 0) cout <<"Recovered charge "<<charge<<endl;
1355                                   MakePed(busPatch->GetBusPatchId(), (Int_t)manuId, (Int_t)channelId, (Int_t)charge);
1356                                 }
1357                             }
1358                           else
1359                             {
1360                               char bpname[256];
1361                               ErrorCounter* errorCounter;
1362                               // Bad buspatch -> not used (just print)
1363                               filcout<<"bpId "<<busPatch->GetBusPatchId()<<" words "<<busPatch->GetLength()
1364                                      <<" parity errors "<<errorCount<<endl;
1365                               // Number of events where this buspatch is missing
1366                               sprintf(bpname,"bp%d",busPatch->GetBusPatchId());
1367                               if (!(errorCounter = (ErrorCounter*)gErrorBuspatchTable->FindObject(bpname)))
1368                                 {
1369                                   // New buspatch
1370                                   errorCounter = new ErrorCounter(busPatch->GetBusPatchId());
1371                                   errorCounter->SetName(bpname);
1372                                   gErrorBuspatchTable->Add(errorCounter);
1373                                 }
1374                               else
1375                                 {
1376                                   // Existing buspatch
1377                                   errorCounter->Increment();
1378                                 }
1379                               // errorCounter->Print();
1380                             } // end of if (!errorCount)
1381                         } // end of while( (busPatch = (AliMUONRawStreamTrackerHP ...
1382                       gNEvents++;
1383                       gNEventsRecovered++;
1384                     } //end of if (recoverParityErrors && eventParityErrors && !eventGlitchErrors&& !eventPaddingErrors)
1385                   else
1386                     {
1387                       // Fatal errors reject the event
1388                       if ( TEST_SYSTEM_ATTRIBUTE( rawReader->GetAttributes(),
1389                                                   ATTR_ORBIT_BC ))
1390                         {
1391                           filcout <<"Event rejected -> Period:"<<EVENT_ID_GET_PERIOD( rawReader->GetEventId() )
1392                                   <<" Orbit:"<<EVENT_ID_GET_ORBIT( rawReader->GetEventId() )
1393                                   <<" BunchCrossing:"<<EVENT_ID_GET_BUNCH_CROSSING( rawReader->GetEventId() )<<endl;
1394                         }
1395                       else
1396                         {
1397                           filcout <<"Event rejected -> nbInRun:"<<EVENT_ID_GET_NB_IN_RUN( rawReader->GetEventId() )
1398                                   <<" burstNb:"<<EVENT_ID_GET_BURST_NB( rawReader->GetEventId() )
1399                                   <<" nbInBurst:"<<EVENT_ID_GET_NB_IN_BURST( rawReader->GetEventId() )<<endl;
1400
1401                         }
1402                     } // end of if (!rawStream->GetGlitchErrors() && !rawStream->GetPaddingErrors() ...
1403                   filcout<<"Number of errors : Glitch "<<eventGlitchErrors
1404                          <<" Parity "<<eventParityErrors
1405                          <<" Padding "<<eventPaddingErrors<<endl;
1406                   filcout<<endl;
1407                 } // end of if (!rawStream->IsErrorMessage())
1408
1409               if (eventGlitchErrors)  gGlitchErrors++;
1410               if (eventParityErrors)  gParityErrors++;
1411               if (eventPaddingErrors) gPaddingErrors++;
1412
1413               delete rawReader;
1414               delete rawStream;
1415
1416             } // while (1)
1417         }
1418
1419
1420
1421       if (gCommand.CompareTo("ped") == 0)
1422         {
1423           sprintf(flatFile,"MUONTRKda_ped_%d.ped",gRunNumber);
1424           if(flatOutputFile.IsNull())flatOutputFile=flatFile;
1425           MakePedStore(flatOutputFile);
1426         }
1427
1428       // option gain -> update root file with pedestal results
1429       // gain + create -> recreate root file
1430       // gain + comp -> update root file and compute gain parameters
1431
1432       if (gCommand.Contains("gain"))
1433         {
1434           MakePedStoreForGain(injCharge);
1435         }
1436
1437
1438       delete gPedestalStore;
1439
1440       delete minuitFit;
1441       TVirtualFitter::SetFitter(0);
1442
1443       timers.Stop();
1444
1445       cout << "\nMUONTRKda : Nb of DATE events           = " << gNDateEvents    << endl;
1446       cout << "MUONTRKda : Nb of Glitch errors         = "   << gGlitchErrors  << endl;
1447       cout << "MUONTRKda : Nb of Parity errors         = "   << gParityErrors  << endl;
1448       cout << "MUONTRKda : Nb of Padding errors        = "   << gPaddingErrors << endl;
1449       cout << "MUONTRKda : Nb of events recovered      = "   << gNEventsRecovered<< endl;
1450       cout << "MUONTRKda : Nb of events without errors = "   << gNEvents-gNEventsRecovered<< endl;
1451       cout << "MUONTRKda : Nb of events used           = "   << gNEvents        << endl;
1452
1453       filcout << "\nMUONTRKda : Nb of DATE events           = " << gNDateEvents    << endl;
1454       filcout << "MUONTRKda : Nb of Glitch errors         = "   << gGlitchErrors << endl;
1455       filcout << "MUONTRKda : Nb of Parity errors         = "   << gParityErrors << endl;
1456       filcout << "MUONTRKda : Nb of Padding errors        = "   << gPaddingErrors << endl;
1457       filcout << "MUONTRKda : Nb of events recovered      = "   << gNEventsRecovered<< endl;
1458       filcout << "MUONTRKda : Nb of events without errors = "   << gNEvents-gNEventsRecovered<< endl;
1459       filcout << "MUONTRKda : Nb of events used           = "   << gNEvents        << endl;
1460
1461       if (gCommand.CompareTo("ped") == 0)
1462         {
1463           cout << "\nMUONTRKda : Output logfile             : " << logOutputFile  << endl;
1464           cout << "MUONTRKda : Pedestal Histo file        : " << gHistoFileName  << endl;
1465           cout << "MUONTRKda : Pedestal file (to SHUTTLE) : " << flatOutputFile << endl;
1466         }
1467       else
1468         {
1469           cout << "\nMUONTRKda : Output logfile          : " << logOutputFile  << endl;
1470           cout << "MUONTRKda : DAC data (root file)    : " << gHistoFileName_gain  << endl;
1471           cout << "MUONTRKda : Dummy file (to SHUTTLE) : " << flatOutputFile << endl;
1472         }
1473
1474     }
1475
1476   // Compute gain parameters
1477
1478
1479   if (gCommand.Contains("comp"))
1480     {
1481       flatOutputFile="";
1482
1483       MakeGainStore();
1484
1485       cout << "\nMUONTRKda : Output logfile          : " << logOutputFile_comp  << endl;
1486       cout << "MUONTRKda : Root Histo. file        : " << gRootFileName  << endl;
1487       cout << "MUONTRKda : Gain file (to SHUTTLE)  : " << flatOutputFile << endl;
1488     }
1489
1490
1491   if(gFES) // Store IN FES
1492     {
1493       printf("\n *****  STORE FILE in FES ****** \n");
1494
1495       // be sure that env variable DAQDALIB_PATH is set in script file
1496       //       gSystem->Setenv("DAQDALIB_PATH", "$DATE_SITE/infoLogger");
1497
1498       if (!flatOutputFile.IsNull())
1499         {
1500           if (gCommand.CompareTo("ped") == 0)
1501             status = daqDA_FES_storeFile(flatOutputFile.Data(),"PEDESTALS");
1502           else
1503             status = daqDA_FES_storeFile(flatOutputFile.Data(),"GAINS");
1504
1505           if (status)
1506             {
1507               printf(" Failed to export file : %d\n",status);
1508             }
1509           else if(gPrintLevel) printf(" %s successfully exported to FES  \n",flatOutputFile.Data());
1510         }
1511     }
1512
1513   filcout.close();
1514
1515   printf("\nExecution time : R:%7.2fs C:%7.2fs\n", timers.RealTime(), timers.CpuTime());
1516
1517   return status;
1518 }