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