MUON/MUONTRKda.cxx

   1 /*
   2
   3 Version 1 for MUONTRKda MUON tracking
   4 Working version for pedestal
   5 Framework for gain computing
   6 (Ch. Finck)
   7
   8
   9 Rem:  AliMUON2DMap stores all channels, even those which are not connected
  10       if pedMean == -1, channel not connected to a pad
  11
  12
  13 */
  14 extern "C" {
  15 #include <daqDA.h>
  16 }
  17
  18 #include "event.h"
  19 #include "monitor.h"
  20
  21 #include <Riostream.h>
  22 #include <stdio.h>
  23 #include <stdlib.h>
  24
  25 //AliRoot
  26 #include "AliMUONRawStreamTracker.h"
  27 #include "AliMUONDspHeader.h"
  28 #include "AliMUONBlockHeader.h"
  29 #include "AliMUONBusStruct.h"
  30 #include "AliMUONDDLTracker.h"
  31 #include "AliMUONV2DStore.h"
  32 #include "AliMUON2DMap.h"
  33 #include "AliMUONCalibParamNF.h"
  34 #include "AliMUONObjectPair.h"
  35 #include "AliMUONVDataIterator.h"
  36 #include "AliMpDDLStore.h"
  37 #include "AliMpIntPair.h"
  38
  39 #include "AliRawReaderDate.h"
  40
  41
  42 //ROOT
  43 #include "TFile.h"
  44 #include "TTree.h"
  45 #include "TH1F.h"
  46 #include "TString.h"
  47 #include "TStopwatch.h"
  48 #include "TMath.h"
  49 #include "TTimeStamp.h"
  50 #include "TGraphErrors.h"
  51 #include "TF1.h"
  52 #include "TROOT.h"
  53 #include "TPluginManager.h"
  54
  55
  56 // global variables
  57 AliMUONV2DStore* pedestalStore =  new AliMUON2DMap(kFALSE);
  58 const Int_t kNchannels = 64;
  59 Int_t nManu = 0;
  60 Int_t nChannel = 0;
  61 UInt_t runNumber = 0;
  62 Int_t nEvents = 0;
  63 TH1F* pedMeanHisto = 0x0;
  64 TH1F* pedSigmaHisto = 0x0;
  65 Char_t histoFileName[256];
  66 TString command("ped");
  67
  68 // funtions
  69
  70 //________________
  71 Double_t fitFunc(Double_t *x, Double_t *par)
  72 {
  73     //fit function for gains
  74
  75     Float_t xx = x[0];
  76     Double_t f;
  77
  78     if (xx < par[3])
  79         f = par[0] + par[1]*xx;
  80     else
  81         f= par[0] + par[1]*xx + par[2]*xx*xx;
  82
  83     return f;
  84 }
  85
  86
  87
  88 //__________
  89 void MakePed(Int_t busPatchId, Int_t manuId, Int_t channelId, Int_t charge)
  90 {
  91
  92     AliMUONVCalibParam* ped =
  93         static_cast<AliMUONVCalibParam*>(pedestalStore->Get(busPatchId, manuId));
  94
  95     if (!ped) {
  96       nManu++;
  97       ped = new AliMUONCalibParamNF(2, kNchannels, -1.); // put default wise -1, not connected channel
  98       pedestalStore->Set(busPatchId, manuId, ped, kFALSE);
  99     }
 100
 101     if (nEvents == 1) {
 102       ped->SetValueAsFloat(channelId, 0, 0.);
 103       ped->SetValueAsFloat(channelId, 1, 0.);
 104     }
 105
 106     Float_t pedMean  = ped->ValueAsFloat(channelId, 0) + charge;
 107     Float_t pedSigma = ped->ValueAsFloat(channelId, 1) + charge*charge;
 108
 109     ped->SetValueAsFloat(channelId, 0, pedMean);
 110     ped->SetValueAsFloat(channelId, 1, pedSigma);
 111
 112 }
 113
 114 //________________
 115 void MakePedStore(TString flatOutputFile = "")
 116 {
 117   TTimeStamp date;
 118   Float_t pedMean;
 119   Float_t pedSigma;
 120   ofstream fileout;
 121   Int_t busPatchId;
 122   Int_t manuId;
 123   Int_t channelId;
 124
 125  // histo
 126   sprintf(histoFileName,"mutrkped-%d.root",runNumber);
 127   TFile*  histoFile = new TFile(histoFileName,"RECREATE","MUON Tracking pedestals");
 128
 129   Char_t name[255];
 130   Char_t title[255];
 131   sprintf(name,"pedmean_allch");
 132   sprintf(title,"Pedestal mean all channels");
 133   Int_t nx = 1000;
 134   Int_t xmin = 0;
 135   Int_t xmax = 1000;
 136   pedMeanHisto = new TH1F(name,title,nx,xmin,xmax);
 137   pedMeanHisto->SetDirectory(histoFile);
 138
 139   sprintf(name,"pedsigma_allch");
 140   sprintf(title,"Pedestal sigma all channels");
 141   nx = 200;
 142   xmin = 0;
 143   xmax = 50;
 144   pedSigmaHisto = new TH1F(name,title,nx,xmin,xmax);
 145   pedSigmaHisto->SetDirectory(histoFile);
 146
 147   TTree* tree = new TTree("t","Pedestal tree");
 148   tree->Branch("bp",&busPatchId,"bp/I");
 149   tree->Branch("manu",&manuId,",manu/I");
 150   tree->Branch("channel",&channelId,",channel/I");
 151
 152   if (!flatOutputFile.IsNull()) {
 153     fileout.open(flatOutputFile.Data());
 154     fileout<<"//===========================================================================" << endl;
 155     fileout<<"//                       Pedestal file calculated by MUONTRKda"<<endl;
 156     fileout<<"//===========================================================================" << endl;
 157     fileout<<"//       * Run           : " << runNumber << endl;
 158     fileout<<"//       * Date          : " << date.AsString("l") <<endl;
 159     fileout<<"//       * Statictics    : " << nEvents << endl;
 160     fileout<<"//       * # of MANUS    : " << nManu << endl;
 161     fileout<<"//       * # of channels : " << nChannel << endl;
 162     fileout<<"//"<<endl;
 163     fileout<<"//---------------------------------------------------------------------------" << endl;
 164     fileout<<"//---------------------------------------------------------------------------" << endl;
 165     fileout<<"//format : BUS_PATCH MANU_ID CHANNEL MEAN SIGMA"<<endl;
 166     fileout<<"//---------------------------------------------------------------------------" << endl;
 167
 168   }
 169
 170   // iterator over pedestal
 171   AliMUONVDataIterator* it = pedestalStore->Iterator();
 172
 173   AliMUONObjectPair* p;
 174
 175   while ( ( p = dynamic_cast<AliMUONObjectPair*>(it->Next() ) ) )
 176   {
 177     AliMUONVCalibParam* ped = dynamic_cast<AliMUONVCalibParam*>(p->Value());
 178     AliMpIntPair* pair      = dynamic_cast<AliMpIntPair*>(p->Key());
 179     busPatchId              = pair->GetFirst();
 180     manuId                  = pair->GetSecond();
 181
 182     for (channelId = 0; channelId < ped->Size() ; ++channelId) {
 183
 184       pedMean  = ped->ValueAsFloat(channelId, 0);
 185
 186       if (pedMean > 0) { // connected channels
 187
 188         ped->SetValueAsFloat(channelId, 0, pedMean/(Float_t)nEvents);
 189
 190         pedMean  = ped->ValueAsFloat(channelId, 0);
 191         pedSigma = ped->ValueAsFloat(channelId, 1);
 192
 193         ped->SetValueAsFloat(channelId, 1, TMath::Sqrt(TMath::Abs(pedSigma/(Float_t)nEvents - pedMean*pedMean)));
 194
 195         pedMean  = ped->ValueAsFloat(channelId, 0) + 0.5 ;
 196         pedSigma = ped->ValueAsFloat(channelId, 1);
 197
 198
 199         if (!flatOutputFile.IsNull()) {
 200           fileout << "\t" << busPatchId << "\t" << manuId <<"\t"<< channelId << "\t"
 201                   << pedMean <<"\t"<< pedSigma << endl;
 202         }
 203
 204         pedMeanHisto->Fill(pedMean);
 205         pedSigmaHisto->Fill(pedSigma);
 206
 207         tree->Fill();
 208       }
 209     }
 210
 211     delete p; // cos create a new object for each iteration
 212   }
 213
 214   // file outputs
 215   if (!flatOutputFile.IsNull())
 216       fileout.close();
 217
 218   histoFile->Write();
 219   histoFile->Close();
 220
 221 //  delete tree;
 222
 223
 224 // not usable need an update of DATE ->5.28, but it compiles
 225 // uncomment when DATE version ok.
 226 // setenv DATE_FES_PATH
 227 // setenv DATE_RUN_NUMBER
 228 // setenv DATE_ROLE_NAME
 229 // setenv DATE_DETECTOR_CODE
 230
 231 //   if (!flatOutputFile.IsNull()) {
 232
 233 //     flatOutputFile.Prepend("./");
 234
 235 //     status = daqDA_FES_storeFile(flatOutputFile.Data(),"MUONTRKda_results");
 236 //     if (status) {
 237 //       printf("Failed to export file : %d\n",status);
 238 //       return -1;
 239 //     }
 240 //  }
 241
 242 }
 243
 244 //________________
 245 void MakePedStoreForGain()
 246 {
 247     // store pedestal map in root file
 248
 249     Int_t injCharge = 200;
 250
 251     TTree* tree = 0x0;
 252
 253     // compute and store pedestals
 254     MakePedStore();
 255
 256     // store in root file
 257     sprintf(histoFileName,"mutrkgain.root");
 258
 259     TString mode("UPDATE");
 260
 261     if (command.Contains("cre")) {
 262         mode = "RECREATE";
 263     }
 264     TFile* histoFile = new TFile(histoFileName, mode.Data(), "MUON Tracking Gains");
 265
 266     // second argument should be the injected charge, taken from config crocus file
 267     // put also info about run number could be usefull
 268     AliMpIntPair* pair   = new AliMpIntPair(runNumber, injCharge);
 269
 270     if (mode.CompareTo("UPDATE") == 0) {
 271       tree = (TTree*)histoFile->Get("t");
 272       tree->SetBranchAddress("run",&pair);
 273       tree->SetBranchAddress("ped",&pedestalStore);
 274
 275     } else {
 276       tree = new TTree("t","Pedestal tree");
 277       tree->Branch("run", "AliMpIntPair",&pair);
 278       tree->Branch("ped", "AliMUON2DMap",&pedestalStore);
 279       tree->SetBranchAddress("run",&pair);
 280       tree->SetBranchAddress("ped",&pedestalStore);
 281
 282     }
 283
 284     tree->Fill();
 285     tree->Write("t", TObject::kOverwrite); // overwrite the tree
 286     histoFile->Close();
 287
 288     delete pair;
 289 }
 290
 291 //________________
 292 void MakeGainStore(TString flatOutputFile)
 293 {
 294
 295     // open file mutrkgain.root
 296     // read again the pedestal for the calibration runs (9 runs ?)
 297     // need the injection charge from config file (to be done)
 298     // For each channel make a TGraphErrors (mean, sigma) vs injected charge
 299     // Fit with a polynomial fct
 300     // store the result in a flat file.
 301
 302     Float_t pedMean[10];
 303     Float_t pedSigma[10];
 304     Float_t injCharge[10];
 305     Float_t injChargeErr[10];
 306
 307     ofstream fileout;
 308     TTimeStamp date;
 309
 310     if (!flatOutputFile.IsNull()) {
 311       fileout.open(flatOutputFile.Data());
 312       fileout<<"//===========================================================================" << endl;
 313       fileout<<"//                       Pedestal file calculated by MUONTRKda"<<endl;
 314       fileout<<"//===========================================================================" << endl;
 315       fileout<<"//       * Run           : " << runNumber << endl;
 316       fileout<<"//       * Date          : " << date.AsString("l") <<endl;
 317       fileout<<"//       * Statictics    : " << nEvents << endl;
 318       fileout<<"//       * # of MANUS    : " << nManu << endl;
 319       fileout<<"//       * # of channels : " << nChannel << endl;
 320       fileout<<"//"<<endl;
 321       fileout<<"//---------------------------------------------------------------------------" << endl;
 322       fileout<<"//---------------------------------------------------------------------------" << endl;
 323       fileout<<"//format : BUS_PATCH MANU_ID CHANNEL GAIN0 GAIN1 GAIN2 XLIM CHI2"<<endl;
 324       fileout<<"//---------------------------------------------------------------------------" << endl;
 325
 326     }
 327
 328
 329     sprintf(histoFileName,"mutrkgain.root");
 330     TFile*  histoFile = new TFile(histoFileName);
 331
 332     AliMUON2DMap* map[10];
 333     AliMUONVCalibParam* ped[10];
 334     AliMpIntPair* pair;
 335     AliMpIntPair* run[10];
 336
 337     //read back from root file
 338     TTree* tree = (TTree*)histoFile->Get("t");
 339     Int_t nEntries = tree->GetEntries();
 340
 341     // read back info
 342     for (Int_t i = 0; i < nEntries; ++i) {
 343       map[i] = 0x0;
 344       run[i] = 0x0;
 345       tree->SetBranchAddress("ped",&map[i]);
 346       tree->SetBranchAddress("run",&run[i]);
 347       tree->GetEvent(i);
 348     }
 349
 350     // Q = f(ADC)
 351     TF1* func = new TF1("func",fitFunc, 0., 2500., 4);
 352     // TF1* func = new TF1("func","pol1");
 353
 354     // iterates over the first pedestal run
 355     AliMUONVDataIterator* it = map[0]->Iterator();
 356
 357     AliMUONObjectPair* p;
 358
 359     while ( ( p = dynamic_cast<AliMUONObjectPair*>(it->Next() ) ) )
 360     {
 361       ped[0]  = dynamic_cast<AliMUONVCalibParam*>(p->Value());
 362       pair    = dynamic_cast<AliMpIntPair*>(p->Key());
 363
 364       Int_t busPatchId = pair->GetFirst();
 365       Int_t manuId     = pair->GetSecond();
 366
 367       // read back pedestal from the other runs for the given (bupatch, manu)
 368       for (Int_t i = 1; i < nEntries; ++i) {
 369         ped[i] = static_cast<AliMUONVCalibParam*>(map[i]->Get(busPatchId, manuId));
 370       }
 371
 372       // compute for each channel the gain parameters
 373       for ( Int_t channelId = 0; channelId < ped[0]->Size() ; ++channelId ) {
 374
 375         Int_t n = 0;
 376         for (Int_t i = 0; i < nEntries; ++i) {
 377
 378           if (!ped[i]) continue; //shouldn't happen.
 379           pedMean[i]      = ped[i]->ValueAsFloat(channelId, 0);
 380           pedSigma[i]     = ped[i]->ValueAsFloat(channelId, 1);
 381           injCharge[i]    = (Float_t)run[i]->GetSecond();
 382           injChargeErr[i] = 1.;
 383
 384           if (pedMean[i] < 0) continue; // not connected
 385
 386           if (pedSigma[i] <= 0) pedSigma[i] = 1.; // should not happen.
 387           n++;
 388         }
 389
 390         if (n > 4) {
 391           // if (n > 1) {
 392           //fit
 393           TGraph *gain = new TGraphErrors(n, pedMean, pedSigma, injCharge, injChargeErr);
 394           //should set some initial parameters
 395           func->SetParameter(0,-300);  // a0
 396           func->SetParameter(1, 1.);   // a1
 397           func->SetParameter(2, 0.00001);// a2
 398           func->SetParameter(3, 1100.); // xlim in ADC
 399
 400           gain->Fit("func","q");
 401           cout  << setw(8) << func->GetParameter(0)  <<"\t"<< setw(8) << func->GetParameter(1) << endl;
 402
 403           if (!flatOutputFile.IsNull()) {
 404             fileout << "\t" << busPatchId << "\t" << manuId <<"\t"<< channelId << "\t"
 405                     << setw(8) << func->GetParameter(0)  <<"\t"<< setw(8) << func->GetParameter(1)
 406                     << setw(8) << func->GetParameter(2)  <<"\t"<< setw(5) << func->GetParameter(3)
 407                     << "\t" << func->GetChisquare() << endl;
 408           }
 409           delete gain;
 410         }
 411
 412       }
 413       delete p;
 414     }
 415
 416     // file outputs for gain
 417     if (!flatOutputFile.IsNull())
 418         fileout.close();
 419
 420     delete func;
 421
 422 }
 423
 424 //*************************************************************//
 425
 426 // main routine
 427 int main(Int_t argc, Char_t **argv)
 428 {
 429
 430     // needed for streamer application
 431     gROOT->GetPluginManager()->AddHandler("TVirtualStreamerInfo",
 432                                           "*",
 433                                           "TStreamerInfo",
 434                                           "RIO",
 435                                           "TStreamerInfo()");
 436
 437
 438     Int_t printLevel = 0;  // Global variable defined as extern in the others .cxx files
 439     Int_t skipEvents = 0;
 440     Int_t maxEvents  = 1000000;
 441     Float_t nSigma = 3;
 442     Int_t threshold = -1;
 443     Char_t inputFile[256];
 444     TString flatOutputFile;
 445     TString crocusOutputFile;
 446     TString crocusConfigFile;
 447
 448     Int_t totalParityErr = 0;
 449     Int_t totalGlitchErr = 0;
 450
 451 // option handler
 452
 453    // decode the input line
 454   for (Int_t i = 1; i < argc; i++) // argument 0 is the executable name
 455   {
 456       Char_t* arg;
 457
 458       arg = argv[i];
 459       if (arg[0] != '-') continue;
 460       switch (arg[1])
 461         {
 462         case 'f' :
 463           i++;
 464           sprintf(inputFile,argv[i]);
 465           break;
 466         case 'a' :
 467           i++;
 468           flatOutputFile = argv[i];
 469           break;
 470         case 'o' :
 471           i++;
 472           crocusOutputFile = argv[i];
 473           break;
 474         case 'c' :
 475           i++;
 476           crocusConfigFile = argv[i];
 477           break;
 478         case 'e' :
 479           i++;
 480           command = argv[i];
 481           break;
 482         case 'd' :
 483           i++;
 484           printLevel=atoi(argv[i]);
 485           break;
 486         case 's' :
 487           i++;
 488           skipEvents=atoi(argv[i]);
 489           break;
 490         case 'n' :
 491           i++;
 492           sscanf(argv[i],"%d",&maxEvents);
 493           break;
 494         case 'p' :
 495           i++;
 496           sscanf(argv[i],"%f",&nSigma);
 497           break;
 498         case 't' :
 499           i++;
 500           sscanf(argv[i],"%d",&threshold);
 501           break;
 502         case 'h' :
 503           i++;
 504           printf("\n******************* %s usage **********************",argv[0]);
 505           printf("\n%s -options, the available options are :",argv[0]);
 506           printf("\n-h help                   (this screen)");
 507           printf("\n");
 508           printf("\n Input");
 509           printf("\n-f <raw data file>        (default = %s)",inputFile);
 510           printf("\n-c <Crocus config. file>  (default = %s)",crocusConfigFile.Data());
 511           printf("\n");
 512           printf("\n Output");
 513           printf("\n-a <Flat ASCII file>      (default = %s)",flatOutputFile.Data());
 514           printf("\n-o <CROUCUS cmd file>     (default = %s)",crocusOutputFile.Data());
 515           printf("\n");
 516           printf("\n Options");
 517           printf("\n-d <print level>          (default = %d)",printLevel);
 518           printf("\n-s <skip events>          (default = %d)",skipEvents);
 519           printf("\n-n <max events>           (default = %d)",maxEvents);
 520           printf("\n-p <n sigmas>             (default = %f)",nSigma);
 521           printf("\n-t <threshold (-1 = no)>  (default = %d)",threshold);
 522           printf("\n-e <execute ped/gain>     (default = %s)",command.Data());
 523           printf("\n-e <gain create>           make gain & create a new root file");
 524           printf("\n-e <gain>                  make gain & update root file");
 525           printf("\n-e <gain compute>          make gain & compute gains");
 526
 527           printf("\n\n");
 528           exit(-1);
 529         default :
 530           printf("%s : bad argument %s (please check %s -h)\n",argv[0],argv[i],argv[0]);
 531           argc = 2; exit(-1); // exit if error
 532         } // end of switch
 533     } // end of for i
 534
 535   // set command to lower case
 536   command.ToLower();
 537
 538   // decoding the events
 539
 540   Int_t status;
 541   Int_t nDateEvents = 0;
 542
 543   void* event;
 544
 545    // containers
 546   AliMUONDDLTracker*       ddlTracker = 0x0;
 547   AliMUONBlockHeader*      blkHeader  = 0x0;
 548   AliMUONDspHeader*        dspHeader  = 0x0;
 549   AliMUONBusStruct*        busStruct  = 0x0;
 550
 551   pedMeanHisto = 0x0;
 552   pedSigmaHisto = 0x0;
 553
 554   TStopwatch timers;
 555
 556   timers.Start(kTRUE);
 557
 558   // once we have a configuration file in db
 559   // copy locally a file from daq detector config db
 560   // The current detector is identified by detector code in variable
 561   // DATE_DETECTOR_CODE. It must be defined.
 562   // If environment variable DAQDA_TEST_DIR is defined, files are copied from DAQDA_TEST_DIR
 563   // instead of the database. The usual environment variables are not needed.
 564   if (!crocusConfigFile.IsNull()) {
 565     status = daqDA_DB_getFile("myconfig", crocusConfigFile.Data());
 566     if (status) {
 567       printf("Failed to get config file : %d\n",status);
 568       return -1;
 569     }
 570   }
 571
 572
 573   status = monitorSetDataSource(inputFile);
 574   if (status) {
 575     cerr << "ERROR : monitorSetDataSource status (hex) = " << hex << status
 576               << " " << monitorDecodeError(status) << endl;
 577     return -1;
 578   }
 579   status = monitorDeclareMp("MUON Tracking monitoring");
 580   if (status) {
 581     cerr << "ERROR : monitorDeclareMp status (hex) = " << hex << status
 582               << " " << monitorDecodeError(status) << endl;
 583     return -1;
 584   }
 585
 586   cout << "MUONTRKda : Reading data from file " << inputFile <<endl;
 587
 588   while(1)
 589   {
 590     if (nEvents >= maxEvents) break;
 591     if (nEvents && nEvents % 100 == 0)
 592         cout<<"Cumulated events " << nEvents << endl;
 593
 594     // check shutdown condition
 595     if (daqDA_checkShutdown())
 596         break;
 597
 598     // Skip Events if needed
 599     while (skipEvents) {
 600       status = monitorGetEventDynamic(&event);
 601       skipEvents--;
 602     }
 603
 604     // starts reading
 605     status = monitorGetEventDynamic(&event);
 606     if (status < 0)  {
 607       cout<<"EOF found"<<endl;
 608       break;
 609     }
 610
 611     nDateEvents++;
 612
 613     // decoding rawdata headers
 614     AliRawReader *rawReader = new AliRawReaderDate(event);
 615
 616     Int_t eventType = rawReader->GetType();
 617     runNumber = rawReader->GetRunNumber();
 618
 619     if (eventType != PHYSICS_EVENT)
 620         continue; // for the moment
 621
 622     nEvents++;
 623
 624     // decoding MUON payload
 625     AliMUONRawStreamTracker* rawStream  = new AliMUONRawStreamTracker(rawReader);
 626
 627     // loops over DDL
 628     while((status = rawStream->NextDDL())) {
 629
 630       if (printLevel) printf("\niDDL %d\n", rawStream->GetDDL());
 631
 632       ddlTracker =  rawStream->GetDDLTracker();
 633
 634       // loop over block (CRT) structure
 635       Int_t nBlock = ddlTracker->GetBlkHeaderEntries();
 636       for(Int_t iBlock = 0; iBlock < nBlock ;iBlock++){
 637
 638         blkHeader = ddlTracker->GetBlkHeaderEntry(iBlock);
 639          if (printLevel) printf("Block Total length %d\n",blkHeader->GetTotalLength());
 640
 641         // loop over DSP (FRT) structure
 642         Int_t nDsp = blkHeader->GetDspHeaderEntries();
 643         for(Int_t iDsp = 0; iDsp < nDsp ;iDsp++){   //DSP loop
 644
 645           dspHeader =  blkHeader->GetDspHeaderEntry(iDsp);
 646            if (printLevel) printf("Dsp length %d\n",dspHeader->GetTotalLength());
 647
 648           // loop over BusPatch structure
 649           Int_t nBusPatch = dspHeader->GetBusPatchEntries();
 650           for(Int_t iBusPatch = 0; iBusPatch < nBusPatch; iBusPatch++) {
 651
 652             busStruct = dspHeader->GetBusPatchEntry(iBusPatch);
 653
 654              if (printLevel) printf("busPatchId %d", busStruct->GetBusPatchId());
 655              if (printLevel) printf(" BlockId %d", busStruct->GetBlockId());
 656              if (printLevel) printf(" DspId %d\n", busStruct->GetDspId());
 657
 658              Int_t busPatchId = busStruct->GetBusPatchId();
 659
 660             // loop over data
 661             Int_t dataSize = busStruct->GetLength();
 662             for (Int_t iData = 0; iData < dataSize; iData++) {
 663
 664               if (nEvents == 1)
 665                   nChannel++;
 666               Int_t manuId     = busStruct->GetManuId(iData);
 667               Int_t channelId  = busStruct->GetChannelId(iData);
 668               Int_t charge     = busStruct->GetCharge(iData);
 669
 670               if (printLevel) printf("manuId: %d, channelId: %d charge: %d\n", manuId,
 671                                       channelId, charge);
 672
 673               MakePed(busPatchId, manuId, channelId, charge);
 674
 675             } // iData
 676           } // iBusPatch
 677         } // iDsp
 678       } // iBlock
 679
 680       totalParityErr += rawStream->GetPayLoad()->GetParityErrors();
 681       totalGlitchErr += rawStream->GetPayLoad()->GetGlitchErrors();
 682     } // NextDDL
 683
 684     delete rawReader;
 685     delete rawStream;
 686
 687   } // while (1)
 688
 689   nEvents -= totalGlitchErr; // skipped event cos of Glitch errors
 690
 691
 692   cout << " Total number of parity errors in the Run: " << totalParityErr << endl;
 693   cout << " Total number of glitch errors in the Run: " << totalGlitchErr << endl;
 694
 695   if (command.CompareTo("ped") == 0)
 696       MakePedStore(flatOutputFile);
 697
 698   // option gain -> update root file with pedestal results
 699   // gain + create -> recreate root file
 700   // gain + comp -> update root file and compute gain parameters
 701
 702   if (command.Contains("gain"))
 703       MakePedStoreForGain();
 704
 705   if (command.Contains("comp"))
 706       MakeGainStore(flatOutputFile);
 707
 708
 709   delete pedestalStore;
 710
 711   timers.Stop();
 712
 713   if (!(crocusConfigFile.IsNull()))
 714       cout << "MUONTRKda : CROCUS command file generated : " << crocusOutputFile.Data() << endl;
 715   else
 716       cout << "MUONTRKda : WARNING no CROCUS command file generated" << endl;
 717
 718   cout << "MUONTRKda : Flat ASCII file generated     : " << flatOutputFile << endl;
 719   cout << "MUONTRKda : Histo file generated          : " << histoFileName  << endl;
 720   cout << "MUONTRKda : Nb of DATE events     = "         << nDateEvents    << endl;
 721   cout << "MUONTRKda : Nb of events used     = "         << nEvents        << endl;
 722
 723   printf("Execution time : R:%7.2fs C:%7.2fs\n", timers.RealTime(), timers.CpuTime());
 724
 725   return status;
 726 }