TRD/AliTRDqaBlackEvents.cxx

   1 /**************************************************************************
   2  * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
   3  *                                                                        *
   4  * Author: The ALICE Off-line Project.                                    *
   5  * Contributors are mentioned in the code where appropriate.              *
   6  *                                                                        *
   7  * Permission to use, copy, modify and distribute this software and its   *
   8  * documentation strictly for non-commercial purposes is hereby granted   *
   9  * withount fee, provided thats the abov copyright notice appears in all   *
  10  * copies and that both the copyright notice and this permission notice   *
  11  * appear in the supporting documentation. The authors make no claims     *
  12  * about the suitability of this software for any purpose. It is          *
  13  * provided "as is without express or implied warranty.                  *
  14  **************************************************************************/
  15
  16 /* $Id: AliTRDqaBlackEvents.cxx 23387 2008-01-17 17:25:16Z cblume $ */
  17
  18 ////////////////////////////////////////////////////////////////////////////
  19 //                                                                        //
  20 //  QA of black events                                                    //
  21 //                                                                        //
  22 //  Author:                                                               //
  23 //    Sylwester Radomski (radomski@physi.uni-heidelberg.de)               //
  24 //                                                                        //
  25 ////////////////////////////////////////////////////////////////////////////
  26
  27 #include "TH1D.h"
  28 #include "TH2D.h"
  29 #include "TH2S.h"
  30 #include "TF1.h"
  31 #include "TFile.h"
  32 #include "TCanvas.h"
  33 #include "TPad.h"
  34 #include "TLatex.h"
  35 #include "TStyle.h"
  36 #include "TGraph.h"
  37
  38 #include "AliRawReader.h"
  39
  40 #include "AliTRDrawStreamOld.h"
  41 #include "AliTRDqaBlackEvents.h"
  42
  43 ClassImp(AliTRDqaBlackEvents)
  44
  45 ///////////////////////////////////////////////////////////////////////////////////////////////////
  46
  47 AliTRDqaBlackEvents::AliTRDqaBlackEvents()
  48   :TObject()
  49   ,fnEvents(0)
  50   ,fCreateFull(0)
  51   ,fThresh(0)
  52   ,fCount(0)
  53   ,fRefEv(0)
  54   ,fRefFileName(0x0)
  55   ,fOccupancy(0)
  56   ,fDetRob(0)
  57   ,fTBEvent(0)
  58   ,fRefHistPed(0)
  59   ,fRefHistNoise(0)
  60   ,fErrorHC(0)
  61   ,fErrorMCM(0)
  62   ,fErrorADC(0)
  63   ,fErrorSMHC(0)
  64   ,fErrorSMMCM(0)
  65   ,fErrorSMADC(0)
  66   ,fErrorGraphHC(0)
  67   ,fErrorGraphMCM(0)
  68   ,fErrorGraphADC(0)
  69   ,fGraphMCM(0)
  70   ,fMcmTracks(0)
  71   ,fMapMCM(0)
  72   ,fFracMCM(0)
  73   ,fSMHCped(0)
  74   ,fSMHCerr(0)
  75   ,fNoiseTotal(0)
  76   ,fPP(0)
  77   ,fMinNoise(0.5)
  78   ,fMaxNoise(2)
  79   ,fFitType(0)
  80 {
  81   //
  82   // Constructor
  83   // to create the histograms call Init()
  84   //
  85
  86   for (Int_t i = 0; i < kDET; i++) {
  87     fPed[i]            = 0x0;
  88     fNoise[i]          = 0x0;
  89     fChPP[i]           = 0x0;
  90     fNPointDist[i]     = 0x0;
  91     fChPed[i]          = 0x0;
  92     fChNoise[i]        = 0x0;
  93     fNPoint[i]         = 0x0;
  94     fData[i]           = 0x0;
  95     fSignal[i]         = 0x0;
  96     fnEntriesRM[i]     = 0x0;
  97     fnEntriesRMDist[i] = 0x0;
  98     fChPedRes[i]       = 0x0;
  99     fChNoiseRes[i]     = 0x0;
 100     fErrorLocHC[i]     = 0x0;
 101     fErrorLocMCM[i]    = 0x0;
 102     fErrorLocADC[i]    = 0x0;
 103   }
 104   for (Int_t i = 0; i < 3; i++) {
 105     fGraphPP[i]        = 0x0;
 106     fSMLink[i]         = 0x0;
 107     fGrLink[i]         = 0x0;
 108     fppThresh[i]       = 0;
 109     fnPP[i]            = 0;
 110     fnLink[i]          = 0;
 111   }
 112   for (Int_t i = 0; i < 2; i++) {
 113     fnErrorHC[i]       = 0;
 114     fnErrorMCM[i]      = 0;
 115     fnErrorADC[i]      = 0;
 116   }
 117   for (Int_t i = 0; i < kSM; i++) {
 118     fSmNoiseRms[i]     = 0x0;
 119     fSmNoiseFit[i]     = 0x0;
 120     fSmPP[i]           = 0x0;
 121   }
 122   for (Int_t i = 0; i < kSM+1; i++) {
 123     fNumberADC[i]      = 0x0;
 124     fnADCinSM[i]       = 0;
 125   }
 126   for (Int_t i = 0; i < 1000; i++) {
 127     fEvNoDist[i]       = 0;
 128   }
 129   for (Int_t i = 0; i < kDET*kROB*kMCM; i++) {
 130     fFullSignal[i]     = 0x0;
 131     fFullCounter[i]    = 0;
 132   }
 133   //strncpy(fRefFileName,"",256);
 134
 135 }
 136
 137 ///////////////////////////////////////////////////////////////////////////////////////////////////
 138
 139 AliTRDqaBlackEvents::AliTRDqaBlackEvents(const AliTRDqaBlackEvents &qa)
 140   :TObject(qa)
 141   ,fnEvents(0)
 142   ,fCreateFull(0)
 143   ,fThresh(0)
 144   ,fCount(0)
 145   ,fRefEv(0)
 146   ,fRefFileName(0x0)
 147   ,fOccupancy(0)
 148   ,fDetRob(0)
 149   ,fTBEvent(0)
 150   ,fRefHistPed(0)
 151   ,fRefHistNoise(0)
 152   ,fErrorHC(0)
 153   ,fErrorMCM(0)
 154   ,fErrorADC(0)
 155   ,fErrorSMHC(0)
 156   ,fErrorSMMCM(0)
 157   ,fErrorSMADC(0)
 158   ,fErrorGraphHC(0)
 159   ,fErrorGraphMCM(0)
 160   ,fErrorGraphADC(0)
 161   ,fGraphMCM(0)
 162   ,fMcmTracks(0)
 163   ,fMapMCM(0)
 164   ,fFracMCM(0)
 165   ,fSMHCped(0)
 166   ,fSMHCerr(0)
 167   ,fNoiseTotal(0)
 168   ,fPP(0)
 169   ,fMinNoise(0.5)
 170   ,fMaxNoise(2)
 171   ,fFitType(0)
 172 {
 173   //
 174   // Copy constructor
 175   // to create the histograms call Init()
 176   //
 177
 178   for (Int_t i = 0; i < kDET; i++) {
 179     fPed[i]            = 0x0;
 180     fNoise[i]          = 0x0;
 181     fChPP[i]           = 0x0;
 182     fNPointDist[i]     = 0x0;
 183     fChPed[i]          = 0x0;
 184     fChNoise[i]        = 0x0;
 185     fNPoint[i]         = 0x0;
 186     fData[i]           = 0x0;
 187     fSignal[i]         = 0x0;
 188     fnEntriesRM[i]     = 0x0;
 189     fnEntriesRMDist[i] = 0x0;
 190     fChPedRes[i]       = 0x0;
 191     fChNoiseRes[i]     = 0x0;
 192     fErrorLocHC[i]     = 0x0;
 193     fErrorLocMCM[i]    = 0x0;
 194     fErrorLocADC[i]    = 0x0;
 195   }
 196   for (Int_t i = 0; i < 3; i++) {
 197     fGraphPP[i]        = 0x0;
 198     fSMLink[i]         = 0x0;
 199     fGrLink[i]         = 0x0;
 200     fppThresh[i]       = 0;
 201     fnPP[i]            = 0;
 202     fnLink[i]          = 0;
 203   }
 204   for (Int_t i = 0; i < 2; i++) {
 205     fnErrorHC[i]       = 0;
 206     fnErrorMCM[i]      = 0;
 207     fnErrorADC[i]      = 0;
 208   }
 209   for (Int_t i = 0; i < kSM; i++) {
 210     fSmNoiseRms[i]     = 0x0;
 211     fSmNoiseFit[i]     = 0x0;
 212     fSmPP[i]           = 0x0;
 213   }
 214   for (Int_t i = 0; i < kSM+1; i++) {
 215     fNumberADC[i]      = 0x0;
 216     fnADCinSM[i]       = 0;
 217   }
 218   for (Int_t i = 0; i < 1000; i++) {
 219     fEvNoDist[i]       = 0;
 220   }
 221   for (Int_t i = 0; i < kDET*kROB*kMCM; i++) {
 222     fFullSignal[i]     = 0x0;
 223     fFullCounter[i]    = 0;
 224   }
 225   //strncpy(fRefFileName,"",256);
 226
 227 }
 228
 229 ///////////////////////////////////////////////////////////////////////////////////////////////////
 230
 231 void AliTRDqaBlackEvents::Init()
 232 {
 233   //
 234   // creates histograms
 235   //
 236
 237   //TFile *file = new
 238   //Info("Init", "Statring");
 239
 240   fnEvents = 0;
 241
 242   // histograms for chambers
 243   for(Int_t det=0; det<kDET; det++) {
 244
 245     fNPoint[det]  = new TH2D(Form("entries_%d", det), "",  16, -0.5, 15.5, 144, -0.5, 143.5);
 246     //fData[det]    = new TH3F(Form("data_%d", det), "", 16, -0.5, 15.5, 144, -0.5, 143.5, 50, -0.5, 49.5);
 247
 248     // pedestal noise maps using RMS and Fit
 249     fChPed[det]   = new TH2D(Form("ped_%d", det), "", 16, -0.5, 15.5, 144, -0.5, 143.5);
 250     fChNoise[det] = new TH2D(Form("noise_%d", det), "", 16, -0.5, 15.5, 144, -0.5, 143.5);
 251
 252     //fChPed[det]   = new TH2D(Form("ped_%d", det), "", 16, -0.5, 15.5, 144, -0.5, 143.5);
 253     //fChNoise[det] = new TH2D(Form("noise_%d", det), "", 16, -0.5, 15.5, 144, -0.5, 143.5);
 254
 255     // distribution per detector
 256     fPed[det]     = new TH1D(Form("pedDist_%d", det), ";pedestals (ADC counts)", 100, 5, 15);
 257     fNoise[det]   = new TH1D(Form("noiseDist_%d", det), ";noise (ADC counts)", 100, 0, 5);
 258     fSignal[det]  = new TH1D(Form("signal_%d", det), ";signal (ADC counts)", 100, -0.5, 99.5);
 259     fChPP[det]    = new TH1D(Form("pp_%d", det), ";pp (ADC)", 200, -0.5, 199.5);
 260
 261     fnEntriesRM[det] = new TH2D(Form("entriesRM_%d", det), ";ROB,MCM", 8, -0.5, 7.5, 16, -0.5, 15.5);
 262
 263     // histograms after reference subtraction
 264     fChPedRes[det]    = new TH2D(Form("pedRef_%d", det), "", 16, -0.5, 15.5, 144, -0.5, 143.5);
 265     fChNoiseRes[det]  = new TH2D(Form("noiseRef_%d", det), "", 16, -0.5, 15.5, 144, -0.5, 143.5);
 266
 267
 268     // error codes
 269     fErrorLocMCM[det] = new TH2D(Form("errorLocMCM_%d", det), "", 16, -0.5, 15.5, 144, -0.5, 143.5);
 270     fErrorLocADC[det] = new TH2D(Form("errorLocADC_%d", det), "", 16, -0.5, 15.5, 144, -0.5, 143.5);
 271     fErrorLocHC[det]  = new TH2D(Form("errorLocHC_%d", det), "", 16, -0.5, 15.5, 144, -0.5, 143.5);
 272   }
 273
 274   // histogram for each MCM
 275   for(Int_t i=0; i < kDET * kROB * kMCM; i++)
 276     fFullCounter[i] = 0;
 277
 278   // histograms from the whole detector
 279   fOccupancy = new TH1D("occupancy", "", 20, -0.5, 19.5);
 280   fDetRob    = new TH2D("DetRob", ";detector;ROB", kDET, -0.5, 539.5, 8, -0.5, 7.5);
 281   fTBEvent   = new TH2D("tbEvent", ";event ID;time bin", 100, -0.5, 99.5, 30, -0.5, 29.5);
 282
 283   // errors statistics and location
 284   fErrorHC  = new TH1D("errorHC", ";error ID;", 18, -3.5, 14.5);
 285   fErrorMCM = new TH1D("errorMCM", ";error ID;", 18, -3.5, 14.5);
 286   fErrorADC = new TH1D("errorADC", ";error ID;", 18, -3.5, 14.5);
 287
 288   fErrorSMHC  = new TH1D("errorSM_HC", ";SM id", 18, -0.5, 17.5);
 289   fErrorSMMCM = new TH1D("errorSM_MCM", ";SM id", 18, -0.5, 17.5);
 290   fErrorSMADC = new TH1D("errorSM_ADC", ";SM id", 18, -0.5, 17.5);
 291
 292
 293   fErrorGraphHC  = new TGraph();
 294   fErrorGraphMCM = new TGraph();
 295   fErrorGraphADC = new TGraph();
 296
 297   fGraphMCM = new TGraph();
 298
 299   for(Int_t i=0; i<3; i++) {
 300     fGraphPP[i] = new TGraph();
 301   }
 302
 303
 304   fMapMCM = new TH2D("mapMCM", ";det;mcm", 540, -0.5, 539.5, kROB*kMCM, -0.5, kROB*kMCM-0.5);
 305   fFracMCM = new TH1D("fracMCM", ";frequency", 100, 0, 1);
 306
 307
 308   fErrorGraphHC->GetHistogram()->SetTitle("fraction of events with HC error;event number");
 309   fErrorGraphMCM->GetHistogram()->SetTitle("fraction of events with MCM error;event number;");
 310   fErrorGraphADC->GetHistogram()->SetTitle("fraction of events with ADC error;event number;");
 311
 312
 313   fSMHCped = new TH2D("smHcPed", ";super module;half chamber", 18, -0.5, 17.5, 60, -0.5, 59.5);
 314
 315   // link monitor
 316   const char *linkName[3] = {"smLink", "smBeaf", "smData"};
 317   const char *linkGrName[3] = {"grSmLink", "grSmBeaf", "grSmData"};
 318   for(Int_t i=0; i<3; i++) {
 319     fSMLink[i] = new TH2D(linkName[i], ";super module;link", 18, -0.5, 17.5, 60, -0.5, 59.5);
 320     fGrLink[i] = new TGraph();
 321     fGrLink[i]->SetName(linkGrName[i]);
 322   }
 323
 324   //fZSsize = new TH1D("zssizeSingle", ";threshold;nADC", 40, -0.5, 39.5);
 325
 326
 327   //Info("Init", "Done");
 328
 329   // number of ADC channels fired per SM and in total
 330   for(Int_t sm=0; sm<kSM+1; sm++)
 331     fNumberADC[sm] = new TGraph();
 332
 333   //
 334   fNoiseTotal = new TH1D("noiseTotal", "noise (ADC)", 250, 0, 10);
 335   fPP = new TH1D("peakPeak", "p-p (ADC)", 200, -0.5, 199.5);
 336
 337   for(Int_t sm=0; sm<kSM; sm++) {
 338     fSmNoiseRms[sm] = new TH1D(Form("noiseRms_sm%d", sm), ";noise from RMS (ADC)", 100, 0, 10);
 339     fSmNoiseFit[sm] = new TH1D(Form("noiseFit_sm%d", sm), ";noise frim Fit (ADC)", 100, 0, 10);
 340     fSmPP[sm] = new TH1D(Form("peakPeak_sm%d", sm), ";peak-peak (ADC)", 200, -0.5, 199.5);
 341   }
 342
 343   // event number consistancy
 344   for(Int_t i=0; i<1000; i++) {
 345     fEvNoDist[i] = new TH1D(Form("mcmEvDist_%d", i), ";#Delta Events", 201, -100.5, 100.5);
 346   }
 347   fRefEv = -1;
 348
 349   fMcmTracks = new TObjArray();
 350
 351   // clean data direct
 352   for(Int_t det=0; det<kDET; det++)
 353     for(Int_t row=0; row<kROW; row++)
 354       for(Int_t pad=0; pad<kPAD; pad++)
 355         for(Int_t ch=0; ch<kCH; ch++) {
 356           fDataDirect[det][row][pad][ch] = 0;
 357           fSignalDirect[det][ch] = 0;  // overdone
 358         }
 359 }
 360
 361
 362 ///////////////////////////////////////////////////////////////////////////////////////////////////
 363
 364 void AliTRDqaBlackEvents::Reset()
 365 {
 366   //
 367   // Resets the histograms
 368   //
 369
 370   for(Int_t i=0; i<kDET; i++) {
 371     //fData[i]->Reset();
 372     fChPed[i]->Reset();
 373     fChNoise[i]->Reset();
 374   }
 375 }
 376
 377
 378 ///////////////////////////////////////////////////////////////////////////////////////////////////
 379
 380 void AliTRDqaBlackEvents::SetRefFile(const char *filename) {
 381
 382   //strncpy(fRefFileName,filename,256);
 383   fRefFileName = filename;
 384
 385 }
 386
 387 ///////////////////////////////////////////////////////////////////////////////////////////////////
 388
 389 void AliTRDqaBlackEvents::ReadRefHists(Int_t det) {
 390   //
 391   // Read the reference histograms
 392   //
 393
 394   fRefHistPed = 0;
 395   fRefHistNoise = 0;
 396
 397   TFile *file = 0x0;
 398   if (fRefFileName) file = TFile::Open(fRefFileName);
 399   if (!file) return;
 400
 401   fRefHistPed   = (TH2D*)file->Get(Form("ped_%d",det));
 402   fRefHistNoise = (TH2D*)file->Get(Form("noise_%d", det));
 403
 404   if (file) file->Close();
 405 }
 406
 407 ///////////////////////////////////////////////////////////////////////////////////////////////////
 408
 409 void AliTRDqaBlackEvents::StartEvent()
 410 {
 411   //
 412   // start an event
 413   //
 414
 415   // clear the mcm data
 416   for(Int_t i=0; i < kDET * kROB * kMCM; i++) {
 417     if (fFullSignal[i]) fFullSignal[i]->Reset();
 418     fFullCounter[i] = 0;
 419   }
 420
 421   for(Int_t i=0; i<2; i++) {
 422     fnErrorHC[i] = 0;
 423     fnErrorMCM[i] = 0;
 424     fnErrorADC[i] = 0;
 425   }
 426
 427
 428   Int_t ppThresh[3] = {10, 20, 40};
 429   for(Int_t i=0; i<3; i++) {
 430     fppThresh[i] = ppThresh[i];
 431     fnPP[i] = 0;
 432     fnLink[i] = 0;
 433   }
 434
 435   for(Int_t sm=0; sm<kSM+1; sm++) fnADCinSM[sm] = 0;
 436
 437   if (fRefEv > 0) fRefEv++;
 438   fEvNoDist[999]->Reset();  // keep only the last event
 439
 440 }
 441
 442 ///////////////////////////////////////////////////////////////////////////////////////////////////
 443
 444 void AliTRDqaBlackEvents::AddBuffer(AliTRDrawStreamOld *data, AliRawReader * const reader)
 445 {
 446
 447   //printf ("try to read data\n");
 448   Int_t nextBuff  = data->NextBuffer();
 449   //printf("done ...\n");
 450
 451   if (nextBuff == 0) return;
 452
 453   Int_t sm = reader->GetEquipmentId() - 1024;
 454   //printf("reading SM %d\n", sm);
 455   AliInfo(Form("reading SM %d", sm));
 456
 457   if (sm < 0 || sm > 17) return;
 458
 459   // lopp over stacks, links ...
 460
 461   for (Int_t istack = 0; istack < 5; istack++) {
 462     for (Int_t ilink = 0; ilink < 12; ilink++) {
 463
 464       //printf("HC = %d %d\n", istack, ilink);
 465
 466       Int_t det = sm * 30 + istack * 6 + ilink/2;
 467
 468       // check if data delivered
 469       if (!(data->IsLinkActiveInStack(istack, ilink))) continue;
 470       fSMLink[0]->Fill(sm, istack * 12 + ilink);
 471       fnLink[0]++;
 472
 473       // check if beaf-beaf
 474       if (data->GetLinkMonitorError(istack, ilink)) {
 475         fSMLink[1]->Fill(sm, istack * 12 + ilink);
 476         fnLink[1]++;
 477         continue;
 478       }
 479
 480       // fill histogram with HC header errors
 481       Int_t nErrHc = 0;
 482       Int_t nErrHcTot = 0;
 483
 484       nErrHc = FillBits(fErrorHC, data->GetH0ErrorCode(istack, ilink), 0);
 485       if (!nErrHc) fErrorHC->Fill(-3);
 486       nErrHcTot += nErrHc;
 487
 488       nErrHc = FillBits(fErrorHC, data->GetH1ErrorCode(istack, ilink), 2);
 489       if (!nErrHc) fErrorHC->Fill(-2);
 490       nErrHcTot += nErrHc;
 491
 492       nErrHc = FillBits(fErrorHC, data->GetHCErrorCode(istack, ilink), 4);
 493       if (!nErrHc) fErrorHC->Fill(-1);
 494       nErrHcTot += nErrHc;
 495
 496       // trending
 497       fnErrorHC[0]++;
 498       if (nErrHcTot > 0) {
 499         fnErrorHC[1]++;
 500         fErrorSMHC->Fill(sm);
 501       }
 502
 503       // data integrity protection
 504
 505       //if (data->GetHCErrorCode(istack, ilink) > 0) continue;
 506       if (data->GetH0ErrorCode(istack, ilink) > 0) continue;
 507       if (data->GetH1ErrorCode(istack, ilink) > 0) continue;
 508
 509       fSMLink[2]->Fill(sm, istack * 12 + ilink);
 510       fnLink[2]++;
 511
 512
 513       for (Int_t imcm = 0; imcm < data->GetHCMCMmax(istack, ilink); imcm++ ){
 514
 515         //printf("mcm = %d %d %d\n", istack, ilink, imcm);
 516
 517         // fill MCM error code
 518
 519         Int_t nErrMcm = 0;
 520         Int_t nErrMcmTot = 0;
 521
 522         nErrMcm = FillBits(fErrorMCM, data->GetMCMhdErrorCode(istack, ilink, imcm), 0);
 523         if (!nErrMcm) fErrorMCM->Fill(-3);
 524         nErrMcmTot += nErrMcm;
 525
 526         nErrMcm = FillBits(fErrorMCM, data->GetMCMADCMaskErrorCode(istack, ilink, imcm), 5);
 527         if (!nErrMcm) fErrorMCM->Fill(-2);
 528         nErrMcmTot += nErrMcm;
 529
 530         nErrMcm = FillBits(fErrorMCM, data->GetMCMErrorCode(istack, ilink, imcm), 10);
 531         if (!nErrMcm) fErrorMCM->Fill(-1);
 532         nErrMcmTot += nErrMcm;
 533
 534         // trending
 535         fnErrorMCM[0]++;
 536         if (nErrMcmTot > 0) {
 537           fnErrorMCM[1]++;
 538           fErrorSMMCM->Fill(sm);
 539         }
 540
 541         // MCM protection
 542         if ( (data->GetMCMhdErrorCode(istack,ilink,imcm)) & 2 ) continue;
 543         //if ((data->GetMCMADCMaskErrorCode(istack,ilink,imcm))) continue;
 544         //if ((data->GetMCMErrorCode(istack,ilink,imcm))) continue;
 545
 546         Int_t mcmEvent = data->GetEventNumber(istack, ilink, imcm);
 547
 548         // set the reference event number
 549         if (fRefEv < 0) {
 550           fRefEv = mcmEvent;
 551           printf("Reference Event Number = %d (%d %d %d)\n", fRefEv, istack, ilink, imcm);
 552         }
 553
 554         // fill event distribution
 555         if (!(fnEvents%10)) {
 556           fEvNoDist[fnEvents/10]->Fill(mcmEvent - fRefEv);
 557         }
 558
 559         fEvNoDist[999]->Fill(mcmEvent - fRefEv);
 560
 561         Int_t mcm = data->GetMCM(istack, ilink, imcm);
 562         Int_t rob = data->GetROB(istack, ilink, imcm);
 563
 564         // create a structure for an MCM if needed
 565         Int_t mcmIndex = det * (kMCM * kROB) + rob * kMCM + mcm;
 566         if (fCreateFull && !fFullSignal[mcmIndex])
 567           fFullSignal[mcmIndex] =
 568             new TH2S(Form("mcm_%d_%d_%d_%d_%d", sm, istack, ilink/2, rob, mcm),
 569                      Form("mcm-%d-%d-%d-%d-%d;ADC;time bin", sm, istack, ilink/2, rob, mcm),
 570                      21, -0.5, 20.5, 30, -0.5, 29.5);
 571
 572
 573         //Int_t zsADC[21][40];
 574         /*
 575           for(Int_t ina=0; ina<21; ina++)
 576           for(Int_t th=0; th<40; th++)
 577           zsADC[ina][th] = 0;
 578         */
 579
 580         // first loop over ADC chanels
 581
 582         for (Int_t iadc=0; iadc < data->GetADCcount(istack, ilink, imcm); iadc++) {
 583
 584           //printf("ADC = %d\n", iadc);
 585
 586
 587           // fill ADC error bits
 588           Int_t nErrAdc = FillBits(fErrorADC, data->GetADCErrorCode(), 0);
 589           if (!nErrAdc) fErrorADC->Fill(-1);
 590
 591           fnErrorADC[0]++;
 592           if (nErrAdc > 0) {
 593             fnErrorADC[1]++;
 594             fErrorSMADC->Fill(sm);
 595           }
 596
 597           // ADC protection
 598           if ((data->GetADCErrorCode(istack,ilink,imcm,iadc))) continue;
 599
 600           Int_t minV = 1024;
 601           Int_t maxV = 0;
 602
 603           Int_t *sig = data->GetSignalDirect(istack, ilink, imcm, iadc);
 604
 605           //Int_t adc = data->GetADCnumber(istack, ilink, imcm, iadc);
 606           Int_t row = data->GetRow(istack, ilink, imcm);
 607           Int_t col = data->GetCol(istack, ilink, imcm, iadc);
 608
 609           // loop over Time Bins and fill histograms
 610           for(Int_t k=0; k < data->GetNumberOfTimeBins(istack, ilink); k++) {
 611
 612             //fSignal[det]->Fill(sig[k]);
 613             //fData[det]->Fill(row, col, sig[k]); // slow
 614
 615             if ((sig[k] >=0) && (sig[k] < kCH)) {
 616               fSignalDirect[det][sig[k]]++;
 617               fDataDirect[det][row][col][sig[k]]++; // direct data
 618             }
 619
 620             // peak-peak
 621             minV = (minV < sig[k]) ? minV : sig[k];
 622             maxV = (maxV > sig[k]) ? maxV : sig[k];
 623
 624             // check for active MCMs
 625             if (fCreateFull && fFullSignal[mcmIndex]) {
 626               if (sig[k] > fThresh || sig[k] < 0) fFullCounter[mcmIndex]++;
 627               //if (sm == 0 && istack == 0 && ilink/2 == 1 && rob == 1 && mcm == 15) fFullCounter[mcmIndex]++; // special
 628               //fFullSignal[mcmIndex]->Fill(adc, k, sig[k]); // slow
 629             }
 630
 631             // zero suppresion tests
 632             /*
 633               for(Int_t th=0; th<40; th++) {
 634               if (sig[k] > th) {
 635               zsADC[iadc][th] = 1;
 636               if (iadc > 0) zsADC[iadc-1][th] = 1;
 637               if (iadc < 10) zsADC[iadc+1][th] = 1;
 638               }
 639               }
 640             */
 641
 642           } // tb
 643
 644           if (maxV > 0) {
 645             fnADCinSM[sm]++;
 646             fnADCinSM[kSM]++;
 647           }
 648
 649           Int_t adcPP = maxV - minV;
 650           //if (adcPP == 100) fFullCounter[mcmIndex] += 10;
 651
 652           fPP->Fill(adcPP);
 653           fChPP[det]->Fill(adcPP);
 654           fSmPP[sm]->Fill(adcPP);
 655
 656           for(Int_t i=0; i<3; i++) {
 657             if ((adcPP) > fppThresh[i]) fnPP[i]++;
 658           }
 659
 660
 661         } // adc
 662
 663           // fill ZS histos
 664           /*
 665             for(Int_t th=0; th<40; th++) {
 666             Int_t nnADC = 0;
 667             for(Int_t ins=0; ins<21; ins++)
 668             nnADC += zsADC[ins][th];
 669             fZSsize->Fill(th, nnADC);
 670             }
 671           */
 672
 673         // fill active MCMs
 674
 675         if (fCreateFull && fFullSignal[mcmIndex] && (fFullCounter[mcmIndex] > fCount)) {
 676
 677           for (Int_t iadc=0; iadc < data->GetADCcount(istack, ilink, imcm); iadc++) {
 678
 679             // ADC protection
 680             if ((data->GetADCErrorCode(istack,ilink,imcm,iadc))) continue;
 681
 682             //Int_t row = data->GetRow(istack, ilink, imcm);
 683             //Int_t col = data->GetCol(istack, ilink, imcm, iadc);
 684             Int_t adc = data->GetADCnumber(istack, ilink, imcm, iadc);
 685
 686             Int_t *sig = data->GetSignalDirect(istack, ilink, imcm, iadc);
 687
 688             // loop over Time Bins and fill histograms
 689             for(Int_t k=0; k < data->GetNumberOfTimeBins(istack, ilink); k++) {
 690               fFullSignal[mcmIndex]->Fill(adc, k, sig[k]); // slow
 691             }
 692           } // tb
 693         }
 694
 695       } // mcm
 696     } // link
 697   } // stack
 698
 699   // printf("end of loops\n");
 700 }
 701
 702
 703 ///////////////////////////////////////////////////////////////////////////////////////////////////
 704
 705 void AliTRDqaBlackEvents::FinishEvent()
 706 {
 707   //
 708   // Processing at the end of the current event
 709   //
 710
 711   for(Int_t i=0; i<3; i++) {
 712     fGraphPP[i]->SetPoint(fnEvents, fnEvents, fnPP[i]);
 713   }
 714
 715   // trend of the number of links
 716   for(Int_t i=0; i<3; i++) {
 717     fGrLink[i]->SetPoint(fnEvents, fnEvents, fnLink[i]);
 718   }
 719
 720   // save interesting histos
 721   Int_t mcmTrackCandidate = 0;
 722   for(Int_t i = 0; i < kDET * kROB * kMCM; i++) {
 723     if ((fFullCounter[i] > fCount) && fFullSignal[i] && CheckMCM(i) )  {
 724
 725       fMcmTracks->AddLast(fFullSignal[i]->Clone(Form("event_%d_%s", fnEvents, fFullSignal[i]->GetName())));
 726       mcmTrackCandidate++;
 727
 728       Int_t mcmTrackletDet = i/(kROB * kMCM);
 729       Int_t mcmTrackletMcm = i%(kROB * kMCM);
 730       fMapMCM->Fill(mcmTrackletDet, mcmTrackletMcm);
 731     }
 732   }
 733
 734   fGraphMCM->SetPoint(fnEvents, fnEvents, mcmTrackCandidate);
 735   AliInfo(Form("Number of MCM track candidates = %d\n", mcmTrackCandidate));
 736
 737
 738   // update fraction of error graphs
 739   Double_t err;
 740
 741   err = (fnErrorHC[0] > 0)? 100.*fnErrorHC[1]/fnErrorHC[0] : -1;
 742   fErrorGraphHC->SetPoint(fnEvents, fnEvents, err);
 743
 744   err = (fnErrorMCM[0] > 0)? 100.*fnErrorMCM[1]/fnErrorMCM[0] : -1;
 745   fErrorGraphMCM->SetPoint(fnEvents, fnEvents, err);
 746
 747   err = (fnErrorADC[0] > 0)? 100.*fnErrorADC[1]/fnErrorADC[0] : -1;
 748   fErrorGraphADC->SetPoint(fnEvents, fnEvents, err);
 749
 750   // number of fired ADC per SM
 751   for(Int_t sm=0; sm<kSM+1; sm++)
 752     fNumberADC[sm]->SetPoint(fnEvents, fnEvents, fnADCinSM[sm]);
 753
 754   fnEvents++;
 755 }
 756
 757 ///////////////////////////////////////////////////////////////////////////////////////////////////
 758
 759 void AliTRDqaBlackEvents::Process(const char *filename)
 760 {
 761   //
 762   // Process something
 763   //
 764
 765   //char fn[256];
 766   //strncpy(fn,filename,256);
 767
 768   //AliInfo(Form("FILENAME = %s (%s)\n", filename, fn));
 769
 770   Int_t map[kDET];
 771
 772   TH1D *hist = new TH1D("fitSignal", "", 50, -0.5, 49.5);
 773   TF1 *fit = new TF1("fit", "gaus(0)", 0, 20);
 774   fit->SetParameters(1e3, 10, 1);
 775
 776   for(Int_t det=0; det<kDET; det++) {
 777
 778     //AliInfo(Form("processing chamber %d\n", det));
 779
 780     map[det] = 0;
 781     //if (fData[det]->GetSum() < 10) continue;
 782     //if (fDataDirect[det][10][10][10] < 20) continue;
 783     //map[det] = 1;
 784
 785
 786     // rewrite signal-direct
 787     for(Int_t ch=0; ch<kCH; ch++) {
 788       fSignal[det]->Fill(ch, fSignalDirect[det][ch]);
 789     }
 790
 791     // read reference distributions
 792     ReadRefHists(det);
 793
 794     //for(Int_t row=0; row<fData[det]->GetXaxis()->GetNbins(); row++) {
 795     //for(Int_t pad=0; pad<fData[det]->GetYaxis()->GetNbins(); pad++) {
 796
 797     for(Int_t row=0; row<kROW; row++) {
 798       for(Int_t pad=0; pad<kPAD; pad++) {
 799
 800         // project the histogramm
 801         hist->Reset();
 802         //for(Int_t bb=0; bb<50; bb++) {
 803         for(Int_t bb=0; bb<kCH; bb++) {
 804           //Int_t dataBin = fData[det]->FindBin(row, pad, bb);
 805           //Double_t v = fData[det]->GetBinContent(dataBin);
 806           hist->SetBinContent(bb+1, fDataDirect[det][row][pad][bb]);
 807         }
 808
 809         Int_t bin = fChPed[det]->FindBin(row, pad);
 810
 811         if (hist->GetSum() > 1) {
 812
 813           map[det] = 1;
 814           Double_t ped = 0, noise = 0;
 815
 816           if (fFitType == 0) {
 817             fit->SetParameters(1e3, 10, 1);
 818             hist->Fit(fit, "q0", "goff", 0, 20);
 819             TF1 *f = hist->GetFunction("fit");
 820             ped = TMath::Abs(f->GetParameter(1));
 821             noise = TMath::Abs(f->GetParameter(2));
 822             fSmNoiseFit[det/30]->Fill(noise);
 823           } else {
 824             ped = hist->GetMean();
 825             noise = hist->GetRMS();
 826             fSmNoiseRms[det/30]->Fill(noise);
 827             //if (pad == 0)
 828             //  AliInfo(Form("data %f %f %f\n", hist->GetSum(), ped, noise));
 829           }
 830
 831           fChPed[det]->SetBinContent(bin, ped);
 832           fChNoise[det]->SetBinContent(bin, noise);
 833           fNoiseTotal->Fill(noise);
 834
 835           // subtract reference values
 836           Double_t refped = 0;
 837           Double_t refnoise = 0;
 838
 839           if (fRefHistPed)   refped   = fRefHistPed->GetBinContent(bin);
 840           if (fRefHistPed)   refnoise = fRefHistPed->GetBinContent(bin);
 841           // Original code, should it not be fRefHistNoise->GetBinContent(bin)
 842           // instead of fRefHistPed->GetBinContent(bin)  (CBL) ???
 843           //if (fRefHistNoise) refnoise = fRefHistPed->GetBinContent(bin);
 844
 845           fChPedRes[det]->SetBinContent(bin, ped-refped);
 846           fChNoiseRes[det]->SetBinContent(bin, noise-refnoise);
 847
 848           fPed[det]->Fill(ped);
 849           fNoise[det]->Fill(noise);
 850
 851           // fill SM-HC plot
 852           Int_t sm = det / 30;
 853           Int_t hc = (pad < kPAD/2) ? 2* (det % 30) : 2* (det % 30) + 1;
 854           if (ped > 9. && ped < 11) fSMHCped->Fill(sm, hc, 1./1152.); // number of pads in HC
 855
 856         } else {
 857
 858           // not enought data found
 859           fChPed[det]->SetBinContent(bin, 0);
 860           fChNoise[det]->SetBinContent(bin, 0);
 861           fChPedRes[det]->SetBinContent(bin, 0);
 862           fChNoiseRes[det]->SetBinContent(bin, 0);
 863         }
 864
 865         //delete hist;
 866       }
 867     }
 868   }
 869
 870
 871   //AliInfo(Form("Number of events = %d\n", fnEvents));
 872
 873   // normalize number of entries histos
 874   Int_t max = 0;
 875   for(Int_t i=0; i<kDET; i++) {
 876     if (!map[i]) continue;
 877     for(Int_t j=0; j<fNPoint[i]->GetXaxis()->GetNbins(); j++) {
 878       for(Int_t k=0; k<fNPoint[i]->GetYaxis()->GetNbins(); k++) {
 879         Int_t dataBin = fNPoint[i]->FindBin(j, k);
 880         Double_t v = fNPoint[i]->GetBinContent(dataBin);
 881         if (v > max) max = (Int_t)v;
 882       }
 883     }
 884   }
 885
 886   char entriesDistName[100];
 887
 888   for(Int_t i=0; i<kDET; i++) {
 889
 890     if (!map[i]) continue;
 891
 892     snprintf(entriesDistName,100,"entriesDist_%d",i);
 893     fNPointDist[i] = new TH1D(entriesDistName, ";number of events", max+2, -0.5, max+1.5);
 894
 895     for(Int_t j=0; j<fNPoint[i]->GetXaxis()->GetNbins(); j++) {
 896       for(Int_t k=0; k<fNPoint[i]->GetYaxis()->GetNbins(); k++) {
 897         Int_t dataBin = fNPoint[i]->FindBin(j, k);
 898         Double_t v = fNPoint[i]->GetBinContent(dataBin);
 899         //if (v > fnEvents) AliInfo(Form("N = %d V = %lf\n", fnEvents, v));
 900         fNPointDist[i]->Fill(v);
 901       }
 902     }
 903
 904     fNPoint[i]->Scale(1./fnEvents);
 905   }
 906
 907
 908   for(Int_t i=0; i<kDET; i++) {
 909     fnEntriesRM[i]->SetMaximum(fnEvents * 1.5);
 910   }
 911
 912   // save histograms
 913
 914   //AliInfo(Form("FILENAME 2 = %s (%d)\n", fn, fn));
 915   TFile *file = new TFile(filename, "recreate");
 916   for(Int_t det = 0; det < kDET; det++) {
 917     if (!map[det]) continue;
 918     fChPed[det]->Write();
 919     fChNoise[det]->Write();
 920     fNPoint[det]->Write();
 921     fNPointDist[det]->Write();
 922     fPed[det]->Write();
 923     fNoise[det]->Write();
 924     fSignal[det]->Write();
 925     fnEntriesRM[det]->Write();
 926     fChPP[det]->Write();
 927
 928     fChPedRes[det]->Write();
 929     fChNoiseRes[det]->Write();
 930
 931     // save error hists
 932     fErrorLocMCM[det]->SetMinimum(0);
 933     fErrorLocMCM[det]->SetMaximum(fnEvents);
 934     fErrorLocMCM[det]->Write();
 935
 936     fErrorLocADC[det]->SetMinimum(0);
 937     fErrorLocADC[det]->SetMaximum(fnEvents);
 938     fErrorLocADC[det]->Write();
 939   }
 940
 941   for(Int_t sm=0; sm<kSM; sm++) {
 942     fSmNoiseRms[sm]->Write();
 943     fSmNoiseFit[sm]->Write();
 944     fSmPP[sm]->Write();
 945   }
 946
 947
 948
 949   Int_t nMcm = 0;
 950   for(Int_t i=0; i < kDET * kROB * kMCM; i++) {
 951     if (fFullSignal[i] && fFullCounter[i] > fCount) {
 952       fFullSignal[i]->Write();
 953       nMcm++;
 954     }
 955   }
 956
 957   AliInfo(Form("Number of saved MCMs = %d\n", nMcm));
 958
 959   fMcmTracks->Write();
 960   AliInfo(Form("Number of tracks = %d\n", fMcmTracks->GetEntries()));
 961
 962   // permanently problematic MCMs
 963   for(Int_t det=0; det<kDET; det++) {
 964     for(Int_t mcm=0; mcm<kROB*kMCM; mcm++) {
 965
 966       Int_t mRob = mcm / kMCM;
 967       Int_t mMcm = mcm % kMCM;
 968       Int_t bin = fMapMCM->FindBin(det, mcm);
 969       Double_t frac = 1. * fMapMCM->GetBinContent(bin) / fnEvents;
 970       fFracMCM->Fill(frac);
 971
 972       if (frac > 0.7) {
 973         AliInfo(Form("{%d, %d, %d, %f}, \n", det, mRob, mMcm, frac));
 974       }
 975     }
 976   }
 977
 978
 979   fOccupancy->Write();
 980   fDetRob->Write();
 981   fTBEvent->Write();
 982
 983   // error hists
 984   fErrorHC->Write();
 985   fErrorMCM->Write();
 986   fErrorADC->Write();
 987
 988   fErrorSMHC->Write();
 989   fErrorSMMCM->Write();
 990   fErrorSMADC->Write();
 991
 992   // write graphs
 993   fErrorGraphHC->Write("trendErrorHC");
 994   fErrorGraphMCM->Write("trendErrorMCM");
 995   fErrorGraphADC->Write("trendErrorADC");
 996
 997   fGraphMCM->Write("trendMCM");
 998
 999   for(Int_t i=0; i<3; i++) {
1000     fGraphPP[i]->Write(Form("fracPP_%d", i));
1001   }
1002
1003   //fZSsize->Scale(1./fnEvents);
1004   //fZSsize->Write();
1005
1006   fMapMCM->SetMaximum(fnEvents);
1007   fMapMCM->Write();
1008   fFracMCM->Write();
1009
1010   fSMHCped->Write();
1011
1012   for(Int_t i=0; i<3; i++ ) {
1013     fSMLink[i]->Write();
1014     fGrLink[i]->Write();
1015   }
1016
1017   for(Int_t sm=0; sm<kSM; sm++)
1018     fNumberADC[sm]->Write(Form("nADCinSM%d",sm));
1019
1020   fNumberADC[kSM]->Write("nADCinEvent");
1021
1022   fNoiseTotal->Write();
1023   fPP->Write();
1024
1025   for(Int_t i=0; i<1000; i++) {
1026     if (fEvNoDist[i]->GetSum() > 0) fEvNoDist[i]->Write();
1027   }
1028
1029
1030   file->Close();
1031   delete file;
1032 }
1033
1034 ///////////////////////////////////////////////////////////////////////////////////////////////////
1035
1036 Int_t AliTRDqaBlackEvents::CheckMCM(Int_t /*index*/) const {
1037   //
1038   // Checks a single MCM
1039   //
1040
1041   return 1;
1042
1043   // static Int_t data[21][3] = {
1044   //   {1, 0, 1},
1045   //   {242, 0, 0},
1046   //   {242, 0, 1},
1047   //   {242, 0, 2},
1048   //   {242, 0, 4},
1049   //   {242, 0, 5},
1050   //   {242, 0, 6},
1051   //   {242, 0, 8},
1052   //   {242, 0, 12},
1053   //   {251, 7, 7},
1054   //   {254, 3, 11},
1055   //   {259, 3, 14},
1056   //   {260, 1, 9},
1057   //   {260, 3, 15},
1058   //   {273, 1, 7},
1059   //   {273, 1, 15},
1060   //   {276, 5, 11},
1061   //   {280, 6, 2},
1062   //   {299, 6, 4},
1063   //   {511, 2, 9},
1064   //   {517, 7, 15}
1065   // };
1066
1067   // for(Int_t i=0; i<21; i++) {
1068   //   Int_t wIndex = data[i][0] * kROB*kMCM + data[i][1] * kMCM + data[i][2];
1069   //   if (index == wIndex) return 0;
1070   // }
1071
1072   return 1;
1073 }
1074
1075 ///////////////////////////////////////////////////////////////////////////////////////////////////
1076
1077
1078
1079
1080
1081
1082
1083 void AliTRDqaBlackEvents::DrawChamber(const char *filename, Int_t det, Int_t w, Int_t h)
1084 {
1085   //
1086   // Draw raport for one chamber:
1087   // pedestal map, noise map, distribution of pedestal and noise
1088   //
1089   // input:
1090   // name of the file with histograms (created with Process())
1091   // detector Id (0 - 539)
1092   //
1093
1094   // setup global style
1095   gStyle->SetPalette(1);
1096   gStyle->SetOptStat(0);
1097   gStyle->SetPadTopMargin(0.02);
1098   gStyle->SetPadBottomMargin(0.05);
1099
1100   TFile *file = new TFile(filename, "READ");
1101
1102   TCanvas *c = new TCanvas("blackEvents",Form("blackEvents %d",det), w, h);
1103   c->SetVertical(kFALSE);
1104   c->Divide(3,1, 0.01, 0.01);
1105   c->cd(3);
1106
1107   TPad *mPad = (TPad*) gPad;
1108   mPad->Divide(1,2,0.01,0.01);
1109
1110   c->cd(1);
1111   TH2D *h2 = (TH2D*)file->Get(Form("ped_%d",det));
1112   h2->SetMinimum(5);
1113   h2->SetMaximum(15);
1114   h2->SetTitle(";Z direction;#phi direction");
1115   h2->Draw("colz");
1116
1117   c->cd(2);
1118   h2 = (TH2D*)file->Get(Form("noise_%d",det));
1119   h2->SetMinimum(fMinNoise);
1120   h2->SetMaximum(fMaxNoise);
1121   h2->SetTitle(";Z direction;#phi direction");
1122   h2->Draw("colz");
1123
1124   mPad->cd(1);
1125   //gPad->SetLogy();
1126   TH1D *h1 = (TH1D*)file->Get(Form("pedDist_%d", det));
1127   h1->Draw();
1128
1129   mPad->cd(2);
1130   gPad->SetLogy();
1131   h1 = (TH1D*)file->Get(Form("noiseDist_%d", det));
1132   h1->Draw();
1133
1134   h1->Fit("gaus");
1135   TF1 *f = h1->GetFunction("gaus");
1136   const char *tt = Form("#mu = %.2f #sigma = %0.2f ", f->GetParameter(1),f->GetParameter(2));
1137   TLatex *ll = new TLatex(2, 100, tt);
1138   ll->SetTextSize(0.06);
1139   ll->Draw();
1140 }
1141
1142 ///////////////////////////////////////////////////////////////////////////////////////////////////
1143
1144 void AliTRDqaBlackEvents::DrawSm(const char *filename, Int_t sm, Int_t w, Int_t h)
1145 {
1146   //
1147   // ????????????
1148   //
1149
1150   gStyle->SetPalette(1);
1151   gStyle->SetOptStat(0);
1152
1153   gStyle->SetPadTopMargin(0.02);
1154   //gStyle->SetPadBottomMargin(0.05);
1155   //gStyle->SetPadLeftMargin(0.02);
1156   //gStyle->SetPadRightMargin(0.02);
1157
1158   TFile *file = new TFile(filename, "READ");
1159
1160   TCanvas *c = new TCanvas("blackEventsSM",Form("blackEvents SM %d",sm), w, h);
1161   c->SetVertical(kFALSE);
1162   c->Divide(5, 6, 0.001, 0.01);
1163
1164   for(Int_t i=0; i<30; i++) {
1165
1166     TH2D *h2 = (TH2D*)file->Get(Form("noise_%d",i+30*sm));
1167     if (!h2) continue;
1168     h2->SetMinimum(fMinNoise);
1169     h2->SetMaximum(fMaxNoise);
1170
1171     // to be replaced by the official calculation
1172     Int_t stack = i/6;
1173     Int_t layer = i%6;
1174     Int_t index = (5-layer)*5 + stack + 1;
1175     //AliInfo(Form("%d %d %d %d\n", i, stack, layer, index));
1176     c->cd(index);
1177     gPad->SetBottomMargin(0.02);
1178     gPad->SetTopMargin(0.02);
1179
1180     h2->Draw("col");
1181   }
1182 }
1183
1184 ///////////////////////////////////////////////////////////////////////////////////////////////////
1185
1186 Int_t AliTRDqaBlackEvents::FillBits(TH1D *hist, Int_t code, Int_t offset) {
1187   //
1188   // Fill bits
1189   //
1190
1191   Int_t nb = 0;
1192   UInt_t test = 1;
1193   for(Int_t i=0; i<8; i++) {
1194     if (code & test) {
1195       hist->Fill(i+offset);
1196       nb++;
1197     }
1198     test *= 2;
1199   }
1200
1201   return nb;
1202 }
1203
1204 ///////////////////////////////////////////////////////////////////////////////////////////////////