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