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