TRD/AliTRDCalibPadStatus.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  * without fee, provided that the above 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$ */
  17
  18 ////////////////////////////////////////////////////////////////////////////
  19 //                                                                        //
  20 // Example: fill pedestal with Gaussian noise                             //
  21 //                                                                        //
  22 //  AliTRDCalibPadStatus ped;                                             //
  23 //  ped.TestEvent(numberofevent);                                         //
  24 //                                                                        //
  25 //  // Method without histo                                               //
  26 //  ped.Analyse();                                                        //
  27 //                                                                        //
  28 //  // Create the histo of the AliTRDCalROC                               //
  29 //  TH2F * histo2dm = ped.GetCalRocMean(0,kFALSE)->MakeHisto2D();         //
  30 //  histo2dm->Scale(10.0);                                                //
  31 //  TH1F * histo1dm = ped.GetCalRocMean(0,kFALSE)->MakeHisto1D();         //
  32 //  histo1dm->Scale(10.0);                                                //
  33 //  TH2F * histo2ds = ped.GetCalRocSquares(0,kFALSE)->MakeHisto2D();      //
  34 //  histo2ds->Scale(10.0);                                                //
  35 //  TH1F * histo1ds = ped.GetCalRocSquares(0,kFALSE)->MakeHisto1D();      //
  36 //  histo1ds->Scale(10.0)                                                 //
  37 //                                                                        //
  38 //  // Draw output                                                        //
  39 //  TCanvas* c1 = new TCanvas;                                            //
  40 //  c1->Divide(2,2);                                                      //
  41 //  c1->cd(1);                                                            //
  42 //  histo2dm->Draw("colz");                                               //
  43 //  c1->cd(2);                                                            //
  44 //  histo1dm->Draw();                                                     //
  45 //  c1->cd(3);                                                            //
  46 //  histo2ds->Draw("colz");                                               //
  47 //  c1->cd(4);                                                            //
  48 //  histo1ds->Draw();                                                     //
  49 //                                                                        //
  50 //  // Method with histo                                                  //
  51 //  ped.AnalyseHisto();                                                   //
  52 //                                                                        //
  53 //  // Take the histo                                                     //
  54 //  TH1F *histo = ped.GetHisto(31);                                       //
  55 //  histo->SetEntries(1);                                                 //
  56 //  histo->Draw();                                                        //
  57 //
  58 // Authors:
  59 //   R. Bailhache (R.Bailhache@gsi.de, rbailhache@ikf.uni-frankfurt.de)
  60 //   J. Book (jbook@ikf.uni-frankfurt.de)
  61 //                                                                                                    //
  62 ////////////////////////////////////////////////////////////////////////////
  63
  64
  65 //Root includes
  66 #include <TObjArray.h>
  67 #include <TH2F.h>
  68 #include <TString.h>
  69 #include <TMath.h>
  70 #include <TRandom.h>
  71
  72 //#include <TRandom.h>
  73 #include <TDirectory.h>
  74 #include <TFile.h>
  75
  76 //AliRoot includes
  77 #include <AliMathBase.h>
  78 #include "AliRawReader.h"
  79 #include "AliRawReaderRoot.h"
  80 #include "AliRawReaderDate.h"
  81
  82 //header file
  83 #include "AliLog.h"
  84 #include "AliTRDCalibPadStatus.h"
  85 #include "AliTRDrawStreamBase.h"
  86 #include "AliTRDgeometry.h"
  87 #include "AliTRDCommonParam.h"
  88 #include "./Cal/AliTRDCalROC.h"
  89 #include "./Cal/AliTRDCalPadStatus.h"
  90 #include "./Cal/AliTRDCalDet.h"
  91 #include "./Cal/AliTRDCalPad.h"
  92 #include "./Cal/AliTRDCalSingleChamberStatus.h"
  93
  94 #include "AliTRDdigitsManager.h"
  95 #include "AliTRDdigitsParam.h"
  96 #include "AliTRDSignalIndex.h"
  97 #include "AliTRDarraySignal.h"
  98 #include "AliTRDarrayADC.h"
  99 #include "AliTRDfeeParam.h"
 100
 101 #include "AliTRDrawStream.h"
 102
 103 #ifdef ALI_DATE
 104 #include "event.h"
 105 #endif
 106
 107 ClassImp(AliTRDCalibPadStatus) /*FOLD00*/
 108
 109 //_____________________________________________________________________
 110 AliTRDCalibPadStatus::AliTRDCalibPadStatus() : /*FOLD00*/
 111   TObject(),
 112   fGeo(0),
 113   fAdcMin(0),
 114   fAdcMax(21),
 115   fDetector(-1),
 116   fNumberOfTimeBins(0),
 117   fCalRocArrayMean(540),
 118   fCalRocArrayRMS(540),
 119   fCalRocArrayMeand(540),
 120   fCalRocArrayRMSd(540),
 121   fHistoArray(540)
 122 {
 123     //
 124     // default constructor
 125     //
 126
 127   fGeo = new AliTRDgeometry();
 128
 129 }
 130
 131 //_____________________________________________________________________
 132 AliTRDCalibPadStatus::AliTRDCalibPadStatus(const AliTRDCalibPadStatus &ped) : /*FOLD00*/
 133   TObject(ped),
 134   fGeo(0),
 135   fAdcMin(ped.GetAdcMin()),
 136   fAdcMax(ped.GetAdcMax()),
 137   fDetector(ped.fDetector),
 138   fNumberOfTimeBins(ped.fNumberOfTimeBins),
 139   fCalRocArrayMean(540),
 140   fCalRocArrayRMS(540),
 141   fCalRocArrayMeand(540),
 142   fCalRocArrayRMSd(540),
 143   fHistoArray(540)
 144 {
 145     //
 146     // copy constructor
 147     //
 148     for (Int_t idet = 0; idet < 540; idet++){
 149         const AliTRDCalROC *calRocMean  = (AliTRDCalROC*)ped.fCalRocArrayMean.UncheckedAt(idet);
 150         const AliTRDCalROC *calRocRMS   = (AliTRDCalROC*)ped.fCalRocArrayRMS.UncheckedAt(idet);
 151         const AliTRDCalROC *calRocMeand = (AliTRDCalROC*)ped.fCalRocArrayMeand.UncheckedAt(idet);
 152         const AliTRDCalROC *calRocRMSd  = (AliTRDCalROC*)ped.fCalRocArrayRMSd.UncheckedAt(idet);
 153         const TH2F         *hped        = (TH2F*)ped.fHistoArray.UncheckedAt(idet);
 154
 155         if ( calRocMean != 0x0 ) fCalRocArrayMean.AddAt(new AliTRDCalROC(*calRocMean), idet);
 156         if ( calRocRMS != 0x0 )  fCalRocArrayRMS.AddAt(new AliTRDCalROC(*calRocRMS), idet);
 157
 158         if ( calRocMeand != 0x0 ) fCalRocArrayMeand.AddAt(new AliTRDCalROC(*calRocMeand), idet);
 159         if ( calRocRMSd != 0x0 )  fCalRocArrayRMSd.AddAt(new AliTRDCalROC(*calRocRMSd), idet);
 160
 161         if ( hped != 0x0 ){
 162           TH2F *hNew = new TH2F(*hped);
 163           hNew->SetDirectory(0);
 164           fHistoArray.AddAt(hNew,idet);
 165         }
 166
 167     }
 168     if (fGeo) {
 169       delete fGeo;
 170     }
 171     fGeo = new AliTRDgeometry();
 172 }
 173
 174 //_____________________________________________________________________
 175 AliTRDCalibPadStatus& AliTRDCalibPadStatus::operator = (const  AliTRDCalibPadStatus &source)
 176 {
 177   //
 178   // assignment operator
 179   //
 180   if (&source == this) return *this;
 181   new (this) AliTRDCalibPadStatus(source);
 182
 183   return *this;
 184 }
 185 //_____________________________________________________________________
 186 AliTRDCalibPadStatus::~AliTRDCalibPadStatus() /*FOLD00*/
 187 {
 188   //
 189   // destructor
 190   //
 191   fCalRocArrayMean.Delete();
 192   fCalRocArrayRMS.Delete();
 193   fCalRocArrayMeand.Delete();
 194   fCalRocArrayRMSd.Delete();
 195   fHistoArray.Delete();
 196   if (fGeo) {
 197     delete fGeo;
 198   }
 199 }
 200 //_____________________________________________________________________
 201 void AliTRDCalibPadStatus::Destroy()
 202 {
 203   //
 204   // Destroy
 205   //
 206   fCalRocArrayMean.Delete();
 207   fCalRocArrayRMS.Delete();
 208   fCalRocArrayMeand.Delete();
 209   fCalRocArrayRMSd.Delete();
 210   fHistoArray.Delete();
 211 }
 212 //_____________________________________________________________________
 213 Int_t AliTRDCalibPadStatus::UpdateHisto(const Int_t icdet, /*FOLD00*/
 214                                         const Int_t icRow,
 215                                         const Int_t icCol,
 216                                         const Int_t csignal,
 217                                         const Int_t crowMax,
 218                                         const Int_t ccold,
 219                                         const Int_t icMcm)
 220 {
 221   //
 222   // Signal filling methode
 223   //
 224   Int_t nbchannel = icRow+icCol*crowMax;
 225
 226   // now the case of double read channel
 227   if(ccold > 0){
 228     nbchannel = (((ccold-1)*8+ icMcm)*crowMax+icRow)+144*crowMax;
 229     //printf("nbchannel %d, ccold %d, icMcm %d, crowMax %d, icRow %d\n",nbchannel,ccold,icMcm,crowMax,icRow);
 230   }
 231
 232   // fast filling methode.
 233   // Attention: the entry counter of the histogram is not increased
 234   //            this means that e.g. the colz draw option gives an empty plot
 235
 236   Int_t bin = 0;
 237
 238   if ( !(((Int_t)csignal>=fAdcMax ) || ((Int_t)csignal<fAdcMin)) )
 239     bin = (nbchannel+1)*(fAdcMax-fAdcMin+2)+((Int_t)csignal-fAdcMin+1);
 240
 241   //GetHisto(icdet,kTRUE)->Fill(csignal,nbchannel);
 242
 243   GetHisto(icdet,kTRUE)->GetArray()[bin]++;
 244
 245   return 0;
 246 }
 247 //_____________________________________________________________________
 248 Int_t AliTRDCalibPadStatus::UpdateHisto2(const Int_t icdet, /*FOLD00*/
 249                                          const Int_t icRow,
 250                                          const Int_t icCol,
 251                                          const Int_t csignal,
 252                                          const Int_t crowMax,
 253                                          const Int_t ccold,
 254                                          const Int_t icMcm,
 255                                          const Int_t icRob
 256                                          )
 257 {
 258   //
 259   // Signal filling methode
 260   //
 261   Int_t nbchannel = icRow+icCol*crowMax;
 262   Int_t mCm = icMcm%4;
 263   Int_t rOb = icRob%2;
 264
 265   // now the case of double read channel
 266   if(ccold > 0){
 267     nbchannel = (((ccold-1)*8+ (mCm+rOb*4))*crowMax+icRow)+144*crowMax;
 268     //printf("nbchannel %d, ccold %d, icMcm %d, crowMax %d, icRow %d\n",nbchannel,ccold,icMcm,crowMax,icRow);
 269   }
 270
 271   // fast filling methode.
 272   // Attention: the entry counter of the histogram is not increased
 273   //            this means that e.g. the colz draw option gives an empty plot
 274
 275   Int_t bin = 0;
 276
 277   if ( !(((Int_t)csignal>=fAdcMax ) || ((Int_t)csignal<fAdcMin)) )
 278     bin = (nbchannel+1)*(fAdcMax-fAdcMin+2)+((Int_t)csignal-fAdcMin+1);
 279
 280   //GetHisto(icdet,kTRUE)->Fill(csignal,nbchannel);
 281
 282   GetHisto(icdet,kTRUE)->GetArray()[bin]++;
 283
 284   return 0;
 285 }
 286 //_____________________________________________________________________
 287 Int_t AliTRDCalibPadStatus::ProcessEvent(AliTRDrawStreamBase *rawStream, Bool_t nocheck)
 288 {
 289   //
 290   // Event Processing loop - AliTRDRawStreamCosmic
 291   // 0 time bin problem or zero suppression
 292   // 1 no input
 293   // 2 input
 294   //
 295
 296   //
 297   // Raw version number:
 298   // [3,31] non zero suppressed
 299   // 2,4 and [32,63] zero suppressed
 300   //
 301
 302   Int_t withInput = 1;
 303
 304   rawStream->SetSharedPadReadout(kTRUE);
 305
 306   if(!nocheck) {
 307
 308     // Check the raw version and if all have the same number of timebins.
 309
 310     while (rawStream->Next()) {
 311
 312       Int_t rawversion = rawStream->GetRawVersion();                     //  current raw version
 313       //printf("Raw version is %d\n",rawversion);
 314
 315       // Could eventually change, have to check with time
 316       if((rawversion < 3) || (rawversion > 31)) {
 317         AliInfo(Form("this is not no-zero-suppressed data, the version is %d",rawversion));
 318         return 0;
 319       }
 320       Int_t idetector  = rawStream->GetDet();                            //  current detector
 321       Int_t iRow       = rawStream->GetRow();                            //  current row
 322       Int_t iRowMax    = rawStream->GetMaxRow();                         //  current rowmax
 323       Int_t iCol       = rawStream->GetCol();                            //  current col
 324       Int_t iADC       = 21-rawStream->GetADC();                         //  current ADC
 325
 326       // It goes in the opposite direction
 327       Int_t col        = 0;
 328       if(iADC == 1) col = 1;
 329       else {
 330         col = TMath::Max(0,(Int_t)(iADC-19));
 331         if(col > 0) col++;
 332       }
 333       Int_t mcm        = (Int_t)(iCol/18);                               //  current group of 18 col pads
 334       if(col > 1) mcm -= 1;
 335       if(col ==1) mcm += 1;
 336
 337       // printf to check
 338       //Bool_t shared = rawStream->IsCurrentPadShared();
 339       //printf("ADC %d, iCol %d, col %d, mcm %d, shared %d\n",iADC,iCol,col,mcm,(Int_t)shared);
 340
 341       // Take the signal
 342       Int_t *signal    = rawStream->GetSignals();                        //  current ADC signal
 343       Int_t nbtimebin  = rawStream->GetNumberOfTimeBins();               //  number of time bins read from data
 344
 345       if((fDetector != -1) && (nbtimebin != fNumberOfTimeBins)) {
 346         AliInfo(Form("the number of time bins is %d, is different from the previous one %d",nbtimebin,fNumberOfTimeBins));
 347         return 0;
 348       }
 349       fNumberOfTimeBins = nbtimebin;
 350       fDetector         = idetector;
 351
 352       for(Int_t k = 0; k < fNumberOfTimeBins; k++){
 353         if(signal[k]>0 && iCol != -1) UpdateHisto(idetector,iRow,iCol,signal[k],iRowMax,col,mcm);
 354       }
 355
 356       withInput = 2;
 357     }
 358   }
 359   else {
 360
 361     while (rawStream->Next()) {
 362
 363       Int_t idetector  = rawStream->GetDet();                            //  current detector
 364       Int_t iRow       = rawStream->GetRow();                            //  current row
 365       Int_t iRowMax    = rawStream->GetMaxRow();                         //  current rowmax
 366       Int_t iCol       = rawStream->GetCol();                            //  current col
 367       Int_t iADC       = 21-rawStream->GetADC();                            //  current ADC
 368
 369       // It goes in the opposite direction
 370       Int_t col        = 0;
 371       if(iADC == 1) col = 1;
 372       else {
 373         col = TMath::Max(0,(Int_t)(iADC-19));
 374         if(col > 0) col++;
 375       }
 376       Int_t mcm        = (Int_t)(iCol/18);                               //  current group of 18 col pads
 377       if(col > 1) mcm -= 1;
 378       if(col ==1) mcm += 1;
 379
 380       // Take the signal
 381       Int_t *signal    = rawStream->GetSignals();                        //  current ADC signal
 382       Int_t nbtimebin = rawStream->GetNumberOfTimeBins();               //  number of time bins read from data
 383
 384
 385       //printf("det %d, row %d, signal[0] %d, signal[1] %d, signal [2] %d\n", idetector, iRow, signal[0], signal[1], signal[2]);
 386
 387       for(Int_t k = 0; k < nbtimebin; k++){
 388         if(signal[k]>0 && iCol != -1) {
 389           UpdateHisto(idetector,iRow,iCol,signal[k],iRowMax,col,mcm);
 390           //printf("Update with det %d, row %d, col %d, signal %d, rowmax %d, col %d, mcm %d\n",idetector,iRow,iCol,signal[n],iRowMax,col,mcm);
 391         }
 392       }
 393
 394       withInput = 2;
 395     }
 396   }
 397
 398   return withInput;
 399 }
 400 //_____________________________________________________________________
 401 Int_t AliTRDCalibPadStatus::ProcessEvent(AliRawReader *rawReader, Bool_t nocheck)
 402 {
 403   //
 404   //  Event processing loop - AliRawReader
 405   //
 406
 407   Int_t result;
 408
 409   rawReader->Select("TRD");
 410
 411   AliTRDrawStreamBase *pstream = AliTRDrawStreamBase::GetRawStream(rawReader);
 412
 413   result = ProcessEvent(pstream, nocheck);
 414
 415   delete pstream;
 416
 417   return result;
 418 }
 419
 420 //_________________________________________________________________________
 421 Int_t AliTRDCalibPadStatus::ProcessEvent(
 422 #ifdef ALI_DATE
 423                                           const eventHeaderStruct *event,
 424                                           Bool_t nocheck
 425 #else
 426                                           const eventHeaderStruct* /*event*/,
 427                                           Bool_t /*nocheck*/
 428
 429 #endif
 430                                           )
 431 {
 432   //
 433   //  process date event
 434   //
 435 #ifdef ALI_DATE
 436     AliRawReader *rawReader = new AliRawReaderDate((void*)event);
 437     Bool_t result=ProcessEvent(rawReader, nocheck);
 438     delete rawReader;
 439     return result;
 440 #else
 441     Fatal("AliTRDCalibPadStatus", "this class was compiled without DATE");
 442     return 0;
 443 #endif
 444
 445 }
 446
 447 //_____________________________________________________________________
 448
 449 Int_t AliTRDCalibPadStatus::ProcessEvent2(AliRawReader *rawReader)
 450 {
 451   //
 452   // RawReader = AliTRDrawStream (Jochen Klein)
 453   //
 454   // Event Processing loop - AliTRDRawStreamCosmic
 455   // 0 time bin problem or zero suppression
 456   // 1 no input
 457   // 2 input
 458   // Raw version number:
 459   // [3,31] non zero suppressed
 460   // 2,4 and [32,63] zero suppressed
 461   //
 462
 463
 464
 465   Int_t withInput = 1;
 466
 467   AliTRDdigitsManager *digitsManager = new AliTRDdigitsManager(kTRUE);
 468   digitsManager->CreateArrays();
 469
 470   AliTRDrawStream *rawStream = new AliTRDrawStream(rawReader);
 471   rawStream->SetDigitsManager(digitsManager);
 472   //rawStream->SetNoErrorWarning();
 473   //rawStream->SetSharedPadReadout(kTRUE);
 474
 475   AliTRDfeeParam *feeParam = AliTRDfeeParam::Instance();
 476
 477   Int_t det    = 0;
 478   while ((det = rawStream->NextChamber(digitsManager, NULL, NULL)) >= 0) { //idetector
 479     if (digitsManager->GetIndexes(det)->HasEntry()) {//QA
 480       //      printf("there is ADC data on this chamber!\n");
 481
 482        AliTRDarrayADC *digits = (AliTRDarrayADC *) digitsManager->GetDigits(det); //mod
 483       if (digits->HasData()) { //array
 484
 485         AliTRDSignalIndex   *indexes = digitsManager->GetIndexes(det);
 486         if (indexes->IsAllocated() == kFALSE) {
 487           AliError("Indexes do not exist!");
 488           break;
 489         }
 490         Int_t iRow  = 0;
 491         Int_t iCol  = 0;
 492         indexes->ResetCounters();
 493
 494         while (indexes->NextRCIndex(iRow, iCol)) { //column,row
 495
 496           AliTRDdigitsParam *digitParam = (AliTRDdigitsParam *)digitsManager->GetDigitsParam();
 497
 498           Int_t mcm          = 0;     // MCM from AliTRDfeeParam
 499           Int_t rob          = 0;     // ROB from AliTRDfeeParam
 500           Int_t extCol       = 0;     // extended column from AliTRDfeeParam
 501           mcm = feeParam->GetMCMfromPad(iRow,iCol);
 502           rob = feeParam->GetROBfromPad(iRow,iCol);
 503
 504           Int_t idetector  = det;                            //  current detector
 505           Int_t iRowMax    = 16;                              //  current rowmax
 506           if(GetStack(det) == 2) iRowMax = 12;
 507
 508           Int_t adc        = 20 - (iCol%18) -1;                 //  current adc
 509           Int_t col        = 0;                              //  col!=0 ->Shared Pad
 510           extCol = feeParam->GetExtendedPadColFromADC(rob,mcm,adc);
 511           //printf("  iCol %d  iRow %d  iRowMax %d  rob %d  mcm %d  adc %d  extCol %d\n",iCol,iRow,iRowMax,rob,mcm,adc,extCol);
 512
 513           // Signal for regular pads
 514           Int_t nbtimebin  = digitParam->GetNTimeBins(idetector);  //  number of time bins read from data
 515           for(Int_t k = 0; k < nbtimebin; k++){
 516             Short_t signal = 0;
 517             signal = digits->GetData(iRow,iCol,k);
 518
 519             if(signal>0) {
 520               UpdateHisto2(idetector,iRow,iCol,signal,iRowMax,col,mcm,rob);
 521             }
 522           }
 523
 524
 525
 526           if((adc==3-1 || adc==20-1 || adc==19-1) && (iCol > 1 && iCol <142)  ) { //SHARED PADS
 527
 528             switch(adc) {
 529             case 2:
 530               adc = 20;                                       //shared Pad adc
 531               mcm = feeParam->GetMCMfromSharedPad(iRow,iCol); //shared Pad mcm
 532               col =  1;
 533               break;
 534             case 19:
 535               adc = 1;                                        //shared Pad adc
 536               mcm = feeParam->GetMCMfromSharedPad(iRow,iCol); //shared Pad mcm
 537               col =  2;
 538               break;
 539             case 18:
 540               adc =  0;                                       //shared Pad adc
 541               mcm = feeParam->GetMCMfromSharedPad(iRow,iCol); //shared Pad mcm
 542               col =  3;
 543               break;
 544             }
 545             rob = feeParam->GetROBfromSharedPad(iRow,iCol);     //shared Pad rob
 546
 547
 548             extCol = feeParam->GetExtendedPadColFromADC(rob,mcm,adc);     //extended pad col via the shared pad rob,mcm and adc
 549
 550             //printf("SHARED PAD ---  iCol %d  iRow %d  rob %d  mcm %d  adc %d  extCol %d  col %d\n",iCol,iRow,rob,mcm,adc,extCol,col);
 551             for(Int_t k = 0; k < nbtimebin; k++){
 552               Short_t signal = 0;
 553               signal = digits->GetDataByAdcCol(iRow,extCol,k);
 554
 555               if(signal>0) {
 556                 UpdateHisto2(idetector,iRow,iCol,signal,iRowMax,col,mcm,rob);
 557               }
 558             }
 559           } //shared pads end
 560
 561
 562           withInput = 2;
 563         }//column,row
 564
 565       }//array
 566     }//QA
 567     digitsManager->ClearArrays(det);
 568   }//idetector
 569   delete digitsManager;
 570   delete rawStream;
 571   return withInput;
 572
 573 }
 574
 575
 576 //_____________________________________________________________________
 577 Bool_t AliTRDCalibPadStatus::TestEventHisto(Int_t nevent, Int_t sm, Int_t ch) /*FOLD00*/
 578 {
 579   //
 580   //  Test event loop
 581   // fill one oroc and one iroc with random gaus
 582   //
 583
 584   gRandom->SetSeed(0);
 585
 586     for (Int_t ism=sm; ism<sm+1; ism++){
 587         for (Int_t ich=ch; ich < ch+1; ich++){
 588             for (Int_t ipl=0; ipl < 6; ipl++){
 589               for(Int_t irow = 0; irow < fGeo->GetRowMax(ipl,ich,ism); irow++){
 590                 for(Int_t icol = 0; icol < fGeo->GetColMax(ipl); icol++){
 591                   for (Int_t iTimeBin=0; iTimeBin<(30*nevent); iTimeBin++){
 592                     Int_t signal=TMath::Nint(gRandom->Gaus(10,1.5));
 593                     if ( signal>0 )UpdateHisto((ipl+ich*6+ism*6*5),irow,icol,signal,fGeo->GetRowMax(ipl,ich,ism),0,0);
 594                   }
 595                 }
 596               }
 597             }
 598         }
 599     }
 600     return kTRUE;
 601 }
 602
 603 //_____________________________________________________________________
 604 TH2F* AliTRDCalibPadStatus::GetHisto(Int_t det, TObjArray *arr, /*FOLD00*/
 605                                   Int_t nbinsY, Float_t ymin, Float_t ymax,
 606                                   const Char_t *type, Bool_t force)
 607 {
 608     //
 609     // return pointer to histogram
 610     // if force is true create a new histogram if it doesn't exist allready
 611     //
 612     if ( !force || arr->UncheckedAt(det) )
 613         return (TH2F*)arr->UncheckedAt(det);
 614
 615     // if we are forced and histogram doesn't yes exist create it
 616     Char_t name[255], title[255];
 617
 618     snprintf(name,255,"hCalib%s%.3d",type,det);
 619     snprintf(title,255,"%s calibration histogram detector %.2d;ADC channel;Channel (pad)",type,det);
 620
 621
 622     Int_t nbchannels = fGeo->GetRowMax(GetLayer(det),GetStack(det),GetSector(det))*fGeo->GetColMax(GetLayer(det));
 623
 624     // we will add 3*8*rowMax channels at the end for the double counted
 625     nbchannels += 3*8*(fGeo->GetRowMax(GetLayer(det),GetStack(det),GetSector(det)));
 626
 627
 628     // new histogram with calib information. One value for each pad!
 629     TH2F* hist = new TH2F(name,title,
 630                           nbinsY, ymin, ymax,
 631                           nbchannels,0,nbchannels
 632                           );
 633     hist->SetDirectory(0);
 634     arr->AddAt(hist,det);
 635     return hist;
 636 }
 637
 638 //_____________________________________________________________________
 639 TH2F* AliTRDCalibPadStatus::GetHisto(Int_t det, Bool_t force) /*FOLD00*/
 640 {
 641     //
 642     // return pointer to histogram
 643     // if force is true create a new histogram if it doesn't exist allready
 644     //
 645     TObjArray *arr = &fHistoArray;
 646     return GetHisto(det, arr, fAdcMax-fAdcMin, fAdcMin-0.5, fAdcMax-0.5, "Pedestal", force);
 647 }
 648
 649 //_____________________________________________________________________
 650 AliTRDCalROC* AliTRDCalibPadStatus::GetCalRoc(Int_t det, TObjArray* arr, Bool_t force) /*FOLD00*/
 651 {
 652     //
 653     // return pointer to ROC Calibration
 654     // if force is true create a new AliTRDCalROC if it doesn't exist allready
 655     //
 656     if ( !force || arr->UncheckedAt(det) )
 657         return (AliTRDCalROC*)arr->UncheckedAt(det);
 658
 659     // if we are forced and histogram doesn't yes exist create it
 660
 661     // new AliTRDCalROC. One value for each pad!
 662     AliTRDCalROC *croc = new AliTRDCalROC(GetLayer(det),GetStack(det));
 663     arr->AddAt(croc,det);
 664     return croc;
 665 }
 666 //_____________________________________________________________________
 667 AliTRDCalROC* AliTRDCalibPadStatus::GetCalRocMean(Int_t det, Bool_t force) /*FOLD00*/
 668 {
 669     //
 670     // return pointer to Carge ROC Calibration
 671     // if force is true create a new histogram if it doesn't exist allready
 672     //
 673     TObjArray *arr = &fCalRocArrayMean;
 674     return GetCalRoc(det, arr, force);
 675 }
 676
 677 //_____________________________________________________________________
 678 AliTRDCalROC* AliTRDCalibPadStatus::GetCalRocRMS(Int_t det, Bool_t force) /*FOLD00*/
 679 {
 680     //
 681     // return pointer to Carge ROC Calibration
 682     // if force is true create a new histogram if it doesn't exist allready
 683     //
 684     TObjArray *arr = &fCalRocArrayRMS;
 685     return GetCalRoc(det, arr, force);
 686 }
 687 //_____________________________________________________________________
 688 AliTRDCalROC* AliTRDCalibPadStatus::GetCalRocMeand(Int_t det, Bool_t force) /*FOLD00*/
 689 {
 690     //
 691     // return pointer to Carge ROC Calibration
 692     // if force is true create a new histogram if it doesn't exist allready
 693     //
 694     TObjArray *arr = &fCalRocArrayMeand;
 695     return GetCalRoc(det, arr, force);
 696 }
 697
 698 //_____________________________________________________________________
 699 AliTRDCalROC* AliTRDCalibPadStatus::GetCalRocRMSd(Int_t det, Bool_t force) /*FOLD00*/
 700 {
 701     //
 702     // return pointer to Carge ROC Calibration
 703     // if force is true create a new histogram if it doesn't exist allready
 704     //
 705     TObjArray *arr = &fCalRocArrayRMSd;
 706     return GetCalRoc(det, arr, force);
 707 }
 708
 709 //_____________________________________________________________________
 710 void AliTRDCalibPadStatus::AnalyseHisto() /*FOLD00*/
 711 {
 712     //
 713     //  Calculate calibration constants
 714     //
 715
 716     Int_t nbinsAdc = fAdcMax-fAdcMin;
 717
 718     TVectorD param(4);
 719     TMatrixD dummy(3,3);
 720
 721     Float_t *arrayHP=0;
 722
 723
 724     for (Int_t idet=0; idet<540; idet++){
 725         TH2F *hP = GetHisto(idet);
 726         if ( !hP ) {
 727           continue;
 728         }
 729
 730         //printf("Entries for %d\n",idet);
 731
 732         AliTRDCalROC *rocMean     = GetCalRocMean(idet,kTRUE);
 733         AliTRDCalROC *rocRMS      = GetCalRocRMS(idet,kTRUE);
 734
 735         arrayHP = hP->GetArray();
 736         Int_t nChannels = rocMean->GetNchannels();
 737
 738         for (Int_t iChannel=0; iChannel<nChannels; iChannel++){
 739             Int_t offset = (nbinsAdc+2)*(iChannel+1)+1;
 740             Double_t ret = AliMathBase::FitGaus(arrayHP+offset,nbinsAdc,fAdcMin-0.5,fAdcMax-0.5,&param,&dummy);
 741             // if the fitting failed set noise and pedestal to 0
 742             if ((ret==-4) || (ret==-1) || (ret==-2)) {
 743                 param[1]=0.0;
 744                 param[2]=0.0;
 745             }
 746             if((param[1]/10.0) > 65534.0) param[1] = 0.0;
 747             if((param[2]/10.0) > 65534.0) param[2] = 0.0;
 748             rocMean->SetValue(iChannel,param[1]/10.0);
 749             rocRMS->SetValue(iChannel,param[2]/10.0);
 750         }
 751
 752         // here we analyse doubled read channels
 753
 754         AliTRDCalROC *rocMeand     = GetCalRocMeand(idet,kTRUE);
 755         AliTRDCalROC *rocRMSd      = GetCalRocRMSd(idet,kTRUE);
 756
 757         Int_t nrows = rocMeand->GetNrows();
 758         Int_t shift = 144*nrows;
 759         Int_t total = shift+3*8*nrows;
 760
 761         for (Int_t iChannel=shift; iChannel<total; iChannel++){
 762             Int_t offset = (nbinsAdc+2)*(iChannel+1)+1;
 763             Double_t ret = AliMathBase::FitGaus(arrayHP+offset,nbinsAdc,fAdcMin-0.5,fAdcMax-0.5,&param,&dummy);
 764             // if the fitting failed set noise and pedestal to 0
 765             if ((ret==-4) || (ret==-1) || (ret==-2)) {
 766                 param[1]=0.0;
 767                 param[2]=0.0;
 768             }
 769             if((param[1]/10.0) > 65534.0) param[1] = 0.0;
 770             if((param[2]/10.0) > 65534.0) param[2] = 0.0;
 771
 772             // here we have to recalculate backward
 773             Int_t nb   = iChannel-shift;
 774             Int_t row  = nb%nrows;
 775             Int_t j    = (Int_t)(nb/nrows);
 776             Int_t imcm = j%8;
 777             Int_t icol = (Int_t)(j/8);
 778
 779             Int_t finalcol = 18*imcm;
 780             if(icol > 0) icol += 17;
 781             else icol = -1;
 782             finalcol += icol;
 783
 784             Int_t channel = row+finalcol*nrows;
 785
 786             //printf("iChannel %d, nrows %d, finalcol %d, row %d, channel %d\n",iChannel,nrows,finalcol,row,channel);
 787             if((finalcol < 0) || (finalcol >= 144)) continue;
 788
 789             rocMeand->SetValue(channel,param[1]/10.0);
 790             rocRMSd->SetValue(channel,param[2]/10.0);
 791         }
 792
 793     }
 794
 795 }
 796
 797 //_______________________________________________________________________________________
 798 AliTRDCalPadStatus* AliTRDCalibPadStatus::CreateCalPadStatus()
 799 {
 800   //
 801   // Create Pad Status out of Mean and RMS values
 802   // The chamber without data are masked, this is the corrected in the preprocessor
 803   //
 804
 805   AliTRDCalPadStatus* obj = new AliTRDCalPadStatus("padstatus", "padstatus");
 806
 807   for (Int_t idet=0; idet<540; ++idet)
 808     {
 809       AliTRDCalSingleChamberStatus *calROC = obj->GetCalROC(idet);
 810
 811
 812       if ( !GetCalRocMean(idet)) {
 813         for(Int_t k = 0; k < calROC->GetNchannels(); k++){
 814           calROC->SetStatus(k,AliTRDCalPadStatus::kNotConnected);
 815         }
 816         continue;
 817       }
 818
 819
 820       //Take the stuff
 821       AliTRDCalROC *calRocMean    = new AliTRDCalROC(*( (AliTRDCalROC *) GetCalRocMean(idet)));
 822       AliTRDCalROC *calRocRMS     = new AliTRDCalROC(*( (AliTRDCalROC *) GetCalRocRMS(idet)));
 823
 824       //Take the stuff second chance
 825       AliTRDCalROC *calRocMeand    = new AliTRDCalROC(*( (AliTRDCalROC *) GetCalRocMeand(idet)));
 826       AliTRDCalROC *calRocRMSd     = new AliTRDCalROC(*( (AliTRDCalROC *) GetCalRocRMSd(idet)));
 827
 828       calRocRMS->Unfold();
 829       calRocRMSd->Unfold();
 830
 831
 832       //Range
 833       Int_t row      = calROC->GetNrows();
 834       Int_t col      = calROC->GetNcols();
 835
 836       Double_t rmsmean       = calRocMean->GetRMS()*10.0;
 837       Double_t meanmean      = calRocMean->GetMean()*10.0;
 838       Double_t meansquares   = calRocRMS->GetMean();
 839
 840
 841       for(Int_t irow = 0; irow < row; irow++){
 842
 843         // for bridged pads
 844         Float_t meanprevious = 0.0;
 845         Float_t rmsprevious  = 0.0;
 846         Float_t mean         = 0.0;
 847         Float_t rms          = 0.0;
 848
 849         for(Int_t icol = 0; icol < col; icol++){
 850
 851           mean     = calRocMean->GetValue(icol,irow)*10.0;
 852           rms      = calRocRMS->GetValue(icol,irow);
 853
 854           if(icol > 0) {
 855             meanprevious     = calRocMean->GetValue((icol -1),irow)*10.0;
 856             rmsprevious      = calRocRMS->GetValue((icol - 1),irow);
 857           }
 858
 859           Bool_t pb = kFALSE;
 860           // masked if two noisy
 861           if((rms <= 0.0001) || (TMath::Abs(mean-meanmean)>(5*rmsmean)) || (TMath::Abs(rms)>(5.0*TMath::Abs(meansquares)))) {
 862
 863             pb = kTRUE;
 864             // look at second chance
 865             Float_t meand     = calRocMeand->GetValue(icol,irow)*10.0;
 866             Float_t rmsd      = calRocRMSd->GetValue(icol,irow);
 867
 868             if((rmsd <= 0.0001) || (TMath::Abs(meand-meanmean)>(5*rmsmean)) || (TMath::Abs(rmsd)>(5.0*TMath::Abs(meansquares)))) {
 869               if((rmsd <= 0.0001) && (rms <= 0.0001)) {
 870                 calROC->SetStatus(icol,irow,AliTRDCalPadStatus::kNotConnected);
 871               }
 872               else {
 873                 calROC->SetStatus(icol, irow, AliTRDCalPadStatus::kMasked);
 874               }
 875             }
 876             else {
 877               calROC->SetStatus(icol, irow, AliTRDCalPadStatus::kReadSecond);
 878             }
 879           }
 880
 881
 882           // bridge if previous pad found something
 883           if(!pb) {
 884             if((meanprevious == mean) && (rmsprevious == rms) && (mean > 0.0001)) {
 885               //printf("mean previous %f, mean %f, rms %f, rmsprevious %f, col %d\n",meanprevious,mean,rms,rmsprevious,icol);
 886               calROC->SetStatus(icol -1 ,irow, AliTRDCalPadStatus::kPadBridgedRight);
 887               calROC->SetStatus(icol ,irow, AliTRDCalPadStatus::kPadBridgedLeft);
 888             }
 889           }
 890
 891         }
 892       }
 893
 894       delete calRocMean;
 895       delete calRocRMS;
 896       delete calRocMeand;
 897       delete calRocRMSd;
 898
 899
 900     }
 901
 902   return obj;
 903
 904 }
 905 //_______________________________________________________________________________________
 906 AliTRDCalPad* AliTRDCalibPadStatus::CreateCalPad()
 907 {
 908   //
 909   // Create Pad Noise out of RMS values
 910   //
 911
 912   AliTRDCalPad* obj = new AliTRDCalPad("PadNoise", "PadNoise");
 913
 914
 915   for (Int_t det=0; det<AliTRDgeometry::kNdet; ++det)  {
 916
 917     AliTRDCalROC *calROC22 = obj->GetCalROC(det);
 918
 919     AliTRDCalROC *calRocRMS     = ((AliTRDCalROC *)GetCalRocRMS(det,kTRUE));
 920
 921     for(Int_t k = 0; k < calROC22->GetNchannels(); k++){
 922       calROC22->SetValue(k,calRocRMS->GetValue(k));
 923     }
 924
 925   }
 926
 927   return obj;
 928
 929 }
 930
 931 //_______________________________________________________________________________________
 932 AliTRDCalDet* AliTRDCalibPadStatus::CreateCalDet() const
 933 {
 934   //
 935   // Create Det Noise correction factor
 936   //
 937
 938   AliTRDCalDet* obj = new AliTRDCalDet("DetNoise", "DetNoise (correction factor)");
 939
 940   for(Int_t l = 0; l < 540; l++){
 941     obj->SetValue(l,10.0);
 942   }
 943
 944   return obj;
 945
 946 }
 947
 948 //_____________________________________________________________________
 949 void AliTRDCalibPadStatus::DumpToFile(const Char_t *filename, const Char_t *dir, Bool_t append) /*FOLD00*/
 950 {
 951     //
 952     //  Write class to file
 953     //
 954
 955     TString sDir(dir);
 956     TString option;
 957
 958     if ( append )
 959         option = "update";
 960     else
 961         option = "recreate";
 962
 963     TDirectory *backup = gDirectory;
 964     TFile f(filename,option.Data());
 965     f.cd();
 966     if ( !sDir.IsNull() ){
 967         f.mkdir(sDir.Data());
 968         f.cd(sDir);
 969     }
 970     this->Write();
 971     f.Close();
 972
 973     if ( backup ) backup->cd();
 974 }
 975
 976 //_____________________________________________________________________
 977 void AliTRDCalibPadStatus::SetCalRocMean(AliTRDCalROC *mean, Int_t det) /*FOLD00*/
 978 {
 979     //
 980     //  Put the AliTRDCalROC in the array fCalRocArrayMean
 981     //
 982
 983
 984   AliTRDCalROC *rocMean = GetCalRocMean(det,kTRUE);
 985
 986   Int_t nChannels = rocMean->GetNchannels();
 987
 988   for (Int_t iChannel=0; iChannel<nChannels; iChannel++){
 989
 990     rocMean->SetValue(iChannel,mean->GetValue(iChannel));
 991
 992   }
 993
 994 }
 995
 996 //_____________________________________________________________________
 997 void AliTRDCalibPadStatus::SetCalRocRMS(AliTRDCalROC *rms, Int_t det) /*FOLD00*/
 998 {
 999     //
1000     //  Put the AliTRDCalROC in the array fCalRocArrayRMS
1001     //
1002
1003
1004   AliTRDCalROC *rocRms = GetCalRocRMS(det,kTRUE);
1005
1006   Int_t nChannels = rocRms->GetNchannels();
1007
1008   for (Int_t iChannel=0; iChannel<nChannels; iChannel++){
1009
1010     rocRms->SetValue(iChannel,rms->GetValue(iChannel));
1011
1012   }
1013
1014 }
1015 //_____________________________________________________________________
1016 void AliTRDCalibPadStatus::SetCalRocMeand(AliTRDCalROC *mean, Int_t det) /*FOLD00*/
1017 {
1018     //
1019     //  Put the AliTRDCalROC in the array fCalRocArrayMean
1020     //
1021
1022
1023   AliTRDCalROC *rocMean = GetCalRocMeand(det,kTRUE);
1024
1025   Int_t nChannels = rocMean->GetNchannels();
1026
1027   for (Int_t iChannel=0; iChannel<nChannels; iChannel++){
1028
1029     rocMean->SetValue(iChannel,mean->GetValue(iChannel));
1030
1031   }
1032
1033 }
1034
1035 //_____________________________________________________________________
1036 void AliTRDCalibPadStatus::SetCalRocRMSd(AliTRDCalROC *rms, Int_t det) /*FOLD00*/
1037 {
1038     //
1039     //  Put the AliTRDCalROC in the array fCalRocArrayRMS
1040     //
1041
1042
1043   AliTRDCalROC *rocRms = GetCalRocRMSd(det,kTRUE);
1044
1045   Int_t nChannels = rocRms->GetNchannels();
1046
1047   for (Int_t iChannel=0; iChannel<nChannels; iChannel++){
1048
1049     rocRms->SetValue(iChannel,rms->GetValue(iChannel));
1050
1051   }
1052
1053 }
1054 //_____________________________________________________________________________
1055 Int_t AliTRDCalibPadStatus::GetLayer(Int_t d) const
1056 {
1057   //
1058   // Reconstruct the layer number from the detector number
1059   //
1060
1061   return ((Int_t) (d % 6));
1062
1063 }
1064
1065 //_____________________________________________________________________________
1066 Int_t AliTRDCalibPadStatus::GetStack(Int_t d) const
1067 {
1068   //
1069   // Reconstruct the chamber number from the detector number
1070   //
1071
1072   return ((Int_t) (d % 30) / 6);
1073
1074 }
1075
1076 //_____________________________________________________________________________
1077 Int_t AliTRDCalibPadStatus::GetSector(Int_t d) const
1078 {
1079   //
1080   // Reconstruct the sector number from the detector number
1081   //
1082
1083   return ((Int_t) (d / 30));
1084
1085 }
1086
1087