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