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