PMD/AliPMDCalibGain.cxx

   1 //////////////////////////////////////////////////////////////////////////////
   2 //
   3 // * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
   4 // *                                                                        *
   5 // * Author: The ALICE Off-line Project.                                    *
   6 // * Contributors are mentioned in the code where appropriate.              *
   7 // *                                                                        *
   8 // * Permission to use, copy, modify and distribute this software and its   *
   9 // * documentation strictly for non-commercial purposes is hereby granted   *
  10 // * without fee, provided that the above copyright notice appears in all   *
  11 // * copies and that both the copyright notice and this permission notice   *
  12 // * appear in the supporting documentation. The authors make no claims     *
  13 // * about the suitability of this software for any purpose. It is          *
  14 // * provided "as is" without express or implied warranty.                  *
  15 // **************************************************************************/
  16 //
  17 //////////////////////////////////////////////////////////////////////////////
  18
  19 #include "TF1.h"
  20 #include "TFile.h"
  21 #include "TObjString.h"
  22 #include "TROOT.h"
  23 #include "TClonesArray.h"
  24 #include "TH1F.h"
  25 #include "TObjArray.h"
  26 #include "TTree.h"
  27 #include "TMath.h"
  28
  29 // --- Standard library ---
  30
  31 // --- AliRoot header files ---
  32 #include "AliDAQ.h"
  33 #include "AliLog.h"
  34 #include "AliRawReader.h"
  35 #include "AliPMDRawStream.h"
  36 #include "AliPMDddldata.h"
  37 #include "AliPMDCalibGain.h"
  38
  39 ClassImp(AliPMDCalibGain)
  40
  41 AliPMDCalibGain::AliPMDCalibGain():
  42   TObject(),
  43   fpw(NULL)
  44 {
  45   // Standard Constructor
  46   for(Int_t idet = 0; idet < kDet; idet++)
  47     {
  48       for(Int_t ismn = 0; ismn < kMaxSMN; ismn++)
  49         {
  50           fSMIso[idet][ismn]   = 0.;
  51           fSMCount[idet][ismn] = 0.;
  52           fCountSm[idet][ismn]=0.;
  53           fTempnhit[idet][ismn]=0.;
  54           fTempnhitSq[idet][ismn]=0.;
  55           for(Int_t jrow = 0; jrow < kMaxRow; jrow++)
  56             {
  57               for(Int_t kcol = 0; kcol < kMaxCol; kcol++)
  58                 {
  59                   fCellIso[idet][ismn][jrow][kcol]    = 0.;
  60                   fCellCount[idet][ismn][jrow][kcol]  = 0.;
  61                   fNhitCell[idet][ismn][jrow][kcol]   = 0.;
  62                   fPedMeanRMS[idet][ismn][jrow][kcol] = 0.;
  63                   fHotFlag[idet][ismn][jrow][kcol]    = 0.;
  64
  65                 }
  66             }
  67         }
  68     }
  69
  70 }
  71 // ------------------------------------------------------------------------ //
  72 AliPMDCalibGain::AliPMDCalibGain(const AliPMDCalibGain &pmdcalibgain):
  73   TObject(pmdcalibgain),
  74   fpw(NULL)
  75 {
  76   for(Int_t idet = 0; idet < kDet; idet++)
  77     {
  78       for(Int_t ismn = 0; ismn < kMaxSMN; ismn++)
  79         {
  80           fSMIso[idet][ismn] = pmdcalibgain.fSMIso[idet][ismn] ;
  81           fSMCount[idet][ismn] = pmdcalibgain.fSMCount[idet][ismn] ;
  82           fCountSm[idet][ismn] = pmdcalibgain.fCountSm[idet][ismn];
  83           fTempnhit[idet][ismn] = pmdcalibgain.fTempnhit[idet][ismn];
  84           fTempnhitSq[idet][ismn] = pmdcalibgain.fTempnhitSq[idet][ismn];
  85           for(Int_t jrow = 0; jrow < kMaxRow; jrow++)
  86             {
  87               for(Int_t kcol = 0; kcol < kMaxCol; kcol++)
  88                 {
  89                   fCellIso[idet][ismn][jrow][kcol] = pmdcalibgain.fCellIso[idet][ismn][jrow][kcol];
  90                   fCellCount[idet][ismn][jrow][kcol] = pmdcalibgain.fCellCount[idet][ismn][jrow][kcol];
  91                   fNhitCell[idet][ismn][jrow][kcol] = pmdcalibgain.fNhitCell[idet][ismn][jrow][kcol];
  92                   fPedMeanRMS[idet][ismn][jrow][kcol] = pmdcalibgain.fPedMeanRMS[idet][ismn][jrow][kcol];
  93                   fHotFlag[idet][ismn][jrow][kcol] = pmdcalibgain.fHotFlag[idet][ismn][jrow][kcol];
  94
  95                 }
  96             }
  97         }
  98     }
  99
 100 }
 101 // ------------------------------------------------------------------------ //
 102 AliPMDCalibGain &AliPMDCalibGain::operator=(const AliPMDCalibGain &pmdcalibgain)
 103 {
 104   if(this != &pmdcalibgain)
 105     {
 106       this->fpw = pmdcalibgain.fpw;
 107       for(Int_t idet = 0; idet < kDet; idet++)
 108         {
 109           for(Int_t ismn = 0; ismn < kMaxSMN; ismn++)
 110             {
 111               fSMIso[idet][ismn]      = pmdcalibgain.fSMIso[idet][ismn];
 112               fSMCount[idet][ismn]    = pmdcalibgain.fSMCount[idet][ismn];
 113               fCountSm[idet][ismn]    = pmdcalibgain.fCountSm[idet][ismn];
 114               fTempnhit[idet][ismn]   = pmdcalibgain.fTempnhit[idet][ismn];
 115               fTempnhitSq[idet][ismn] = pmdcalibgain.fTempnhitSq[idet][ismn];
 116               for(Int_t jrow = 0; jrow < kMaxRow;jrow++)
 117                 {
 118                   for(Int_t kcol = 0; kcol < kMaxCol; kcol++)
 119                     {
 120                       fCellIso[idet][ismn][jrow][kcol] = pmdcalibgain.fCellIso[idet][ismn][jrow][kcol];
 121                       fCellCount[idet][ismn][jrow][kcol] = pmdcalibgain.fCellCount[idet][ismn][jrow][kcol];
 122                       fNhitCell[idet][ismn][jrow][kcol] = pmdcalibgain.fNhitCell[idet][ismn][jrow][kcol]; //za
 123                       fPedMeanRMS[idet][ismn][jrow][kcol] = pmdcalibgain.fPedMeanRMS[idet][ismn][jrow][kcol];
 124                       fHotFlag[idet][ismn][jrow][kcol] = pmdcalibgain.fHotFlag[idet][ismn][jrow][kcol];
 125
 126                     }
 127                 }
 128             }
 129         }
 130     }
 131   return *this;
 132 }
 133 // ------------------------------------------------------------------------ //
 134 AliPMDCalibGain::~AliPMDCalibGain()
 135 {
 136   // dtor
 137 }
 138
 139 // ------------------------------------------------------------------------ //
 140
 141 Int_t AliPMDCalibGain::ExtractPedestal(const Char_t *rootFile)
 142 {
 143   // Pedestal extraction from the PMD_PED.root file
 144   // To be called once at the beginning
 145
 146   Int_t   det, sm, row, col;
 147   Float_t mean, rms;
 148
 149   TFile *pedfile = new TFile(rootFile);
 150
 151   if(!pedfile)
 152     {
 153       printf("ERROR --- NO PEDESTAL (PMD_PED1.root) FILE IS FOUND --- STOP GAIN DA\n");
 154       return -3;
 155     }
 156
 157
 158   TTree *ped =(TTree*)pedfile->Get("ped");
 159
 160   ped->SetBranchAddress("det",&det);
 161   ped->SetBranchAddress("sm",&sm);
 162   ped->SetBranchAddress("row",&row);
 163   ped->SetBranchAddress("col",&col);
 164   ped->SetBranchAddress("mean",&mean);
 165   ped->SetBranchAddress("rms",&rms);
 166
 167   Int_t nentries = (Int_t)ped->GetEntries();
 168
 169   for (Int_t ient = 0; ient < nentries; ient++)
 170     {
 171       ped->GetEntry(ient);
 172       fPedMeanRMS[det][sm][row][col] = mean + 3.*rms;
 173       //printf("Mean= %f, RMS= %f, PedMeanRMS=%f\n",mean,rms,fPedMeanRMS[det][sm][row][col]);
 174
 175     }
 176
 177   pedfile->Close();
 178   delete pedfile;
 179   pedfile = 0x0;
 180
 181   return 1;
 182 }
 183 //------------------------------------------------------------------------------------------------
 184
 185 Int_t AliPMDCalibGain::ExtractHotChannel(const Char_t *rootFile)
 186 {
 187   // HotChannel extraction from the PMD_HOT.root file
 188   // To be called once at the beginning
 189
 190   Int_t   det, sm, row, col;
 191   Float_t flag;
 192
 193   TFile *hotmapfile = new TFile(rootFile);
 194
 195   if(!hotmapfile)
 196     {
 197       printf(" NO HOTCHANNEL MAP FOUND (PMD_HOT.root) FILE IS FOUND \n");
 198       fHotFlag[kDet][kMaxSMN][kMaxRow][kMaxCol] = 0.;
 199     }
 200
 201
 202   TTree *hot =(TTree*)hotmapfile->Get("hot");
 203
 204   hot->SetBranchAddress("det",&det);
 205   hot->SetBranchAddress("sm",&sm);
 206   hot->SetBranchAddress("row",&row);
 207   hot->SetBranchAddress("col",&col);
 208   hot->SetBranchAddress("flag",&flag);
 209
 210   Int_t nentries = (Int_t)hot->GetEntries();
 211
 212   for (Int_t ient = 0; ient < nentries; ient++)
 213     {
 214       hot->GetEntry(ient);
 215       fHotFlag[det][sm][row][col] = flag;
 216
 217       //printf(" HotFlag=%f\n",fHotFlag[det][sm][row][col]);
 218     }
 219
 220   hotmapfile->Close();
 221   delete hotmapfile;
 222   hotmapfile = 0x0;
 223
 224   return 1;
 225 }
 226
 227
 228 // ------------------------------------------------------------------------ //
 229
 230 void AliPMDCalibGain::ReadTempFile(const Char_t *tempFile)
 231 {
 232   // Read the variables from the file
 233
 234   fpw = fopen(tempFile,"r");
 235
 236   Float_t smcount, smiso;
 237   Float_t cellcount, celliso;
 238
 239
 240   for (Int_t idet = 0; idet < kDet; idet++)
 241     {
 242       for (Int_t ism = 0; ism < kMaxSMN; ism++)
 243         {
 244           fscanf(fpw,"%d %d %f %f",&idet,&ism,&smcount,&smiso);
 245
 246           fSMCount[idet][ism] = smcount;
 247           fSMIso[idet][ism]   = smiso;
 248         }
 249     }
 250
 251   for (Int_t idet = 0; idet < kDet; idet++)
 252     {
 253       for (Int_t ism = 0; ism < kMaxSMN; ism++)
 254         {
 255           for (Int_t irow = 0; irow < kMaxRow; irow++)
 256             {
 257               for (Int_t icol = 0; icol < kMaxCol; icol++)
 258                 {
 259                   fscanf(fpw,"%d %d %d %d %f %f",&idet,&ism,&irow,&icol,
 260                          &cellcount,&celliso);
 261
 262                   fCellCount[idet][ism][irow][icol] = cellcount;
 263                   fCellIso[idet][ism][irow][icol]   = celliso;
 264                 }
 265             }
 266         }
 267     }
 268
 269   fclose(fpw);
 270
 271 }
 272 // ------------------------------------------------------------------------ //
 273 void AliPMDCalibGain::WriteTempFile(const Char_t *tempFile)
 274 {
 275   // Write the Temporary file if the required statics is not achieved
 276
 277
 278   /*
 279     Following variables to be written to a file
 280     fDetIso[idet] ;
 281     fSMIso[idet][ismn];
 282     fCellIso[idet][ismn][irow][icol];
 283
 284     fDetCount[idet];
 285     fSMCount[idet][ismn];
 286     fCellCount[idet][ismn][irow][icol];
 287   */
 288
 289
 290   fpw = fopen(tempFile,"w+");
 291
 292   for (Int_t idet = 0; idet < kDet; idet++)
 293     {
 294      // fprintf(fpw,"%d %f %f\n",idet,fDetCount[idet],fDetIso[idet]);
 295     }
 296
 297   for (Int_t idet = 0; idet < kDet; idet++)
 298     {
 299       for (Int_t ism = 0; ism < kMaxSMN; ism++)
 300         {
 301           fprintf(fpw,"%d %d %f %f\n",idet,ism, fSMCount[idet][ism],fSMIso[idet][ism]);
 302         }
 303     }
 304
 305   for (Int_t idet = 0; idet < kDet; idet++)
 306     {
 307       for (Int_t ism = 0; ism < kMaxSMN; ism++)
 308         {
 309           for (Int_t irow = 0; irow < kMaxRow; irow++)
 310             {
 311               for (Int_t icol = 0; icol < kMaxCol; icol++)
 312                 {
 313                   fprintf(fpw,"%d %d %d %d %f %f\n",idet,ism,irow,icol,
 314                           fCellCount[idet][ism][irow][icol],
 315                           fCellIso[idet][ism][irow][icol]);
 316                 }
 317             }
 318         }
 319     }
 320
 321   fclose(fpw);
 322
 323 }
 324
 325 // ------------------------------------------------------------------------ //
 326
 327 Bool_t AliPMDCalibGain::ProcessEvent(AliRawReader *rawReader, TObjArray *pmdddlcont)
 328 {
 329   // Calculates the ADC of isolated cell
 330
 331   const Int_t kDDL           = AliDAQ::NumberOfDdls("PMD");
 332   const Int_t kCellNeighbour = 6;
 333   Int_t neibx[6] = {1,0,-1,-1,0,1};
 334   Int_t neiby[6] = {0,1,1,0,-1,-1};
 335
 336   Int_t id1,jd1;  //neighbour row/col
 337   Int_t isocount; //number of neighbours with 0 signal
 338
 339   Float_t d1[kDet][kMaxSMN][kMaxRow][kMaxCol];
 340
 341   for(Int_t idet = 0; idet < kDet; idet++)
 342     {
 343       for(Int_t ismn = 0; ismn < kMaxSMN; ismn++)
 344         {
 345           for(Int_t irow = 0; irow < kMaxRow; irow++)
 346             {
 347               for(Int_t icol = 0; icol < kMaxCol; icol++)
 348                 {
 349                   d1[idet][ismn][irow][icol] = 0.;
 350                 }
 351             }
 352         }
 353     }
 354
 355   AliPMDRawStream rawStream(rawReader);
 356
 357   Int_t iddl = -1;
 358
 359   Int_t numberofDDLs = 0;
 360
 361   while ((iddl = rawStream.DdlData(pmdddlcont)) >=0) {
 362     numberofDDLs++;
 363
 364     Int_t ientries = pmdddlcont->GetEntries();
 365
 366     for (Int_t ient = 0; ient < ientries; ient++)
 367       {
 368         AliPMDddldata *pmdddl = (AliPMDddldata*)pmdddlcont->UncheckedAt(ient);
 369
 370         Int_t idet = pmdddl->GetDetector();
 371         Int_t ismn = pmdddl->GetSMN();
 372         Int_t mcm = pmdddl->GetMCM();
 373         //Int_t ichno = pmdddl->GetChannel();
 374         Int_t irow = pmdddl->GetRow();
 375         Int_t icol = pmdddl->GetColumn();
 376         Int_t isig = pmdddl->GetSignal();
 377
 378         // This is the protection not to crash the code
 379
 380         if(mcm == 0) continue;
 381         if (irow < 0 || icol < 0 || irow > 47 || icol > 95) continue;
 382
 383         // Pedestal subtraction
 384         if(fHotFlag[idet][ismn][irow][icol] == 1.0) isig = 0;
 385
 386         if (isig>0)
 387           {
 388             d1[idet][ismn][irow][icol] =
 389               (Float_t) isig - fPedMeanRMS[idet][ismn][irow][icol];
 390
 391             //printf("Signal_ped_subtracted=%f, pedestal=%f\n",d1[idet][ismn][irow][icol]),fPedMeanRMS[idet][ismn][irow][icol];
 392
 393             fNhitCell[idet][ismn][irow][icol]++;     //cell hit frequency
 394           }
 395
 396
 397       }
 398     pmdddlcont->Delete();
 399     }
 400
 401   for(Int_t idet=0; idet < kDet; idet++)
 402     {
 403       for(Int_t ismn = 0; ismn < kMaxSMN; ismn++)
 404         {
 405           for(Int_t irow = 0; irow < kMaxRow; irow++)
 406             {
 407               for(Int_t icol = 0; icol < kMaxCol; icol++)
 408                 {
 409                   if(d1[idet][ismn][irow][icol] > 0)
 410                     {
 411                       isocount = 0;
 412                       for(Int_t ii = 0; ii < kCellNeighbour; ii++)
 413                         {
 414                           id1 = irow + neibx[ii];
 415                           jd1 = icol + neiby[ii];
 416                           if (id1 < 0) id1 = 0;
 417                           if (id1 > kMaxRow-1) id1 = kMaxRow - 1;
 418                           if (jd1 < 0) jd1 = 0;
 419                           if (jd1 > kMaxCol-1) jd1 = kMaxCol - 1;
 420                           if(d1[idet][ismn][id1][jd1] == 0)
 421                             {
 422                               isocount++;
 423                               if(isocount == kCellNeighbour)
 424                                 {
 425                                   //fDetIso[idet] += d1[idet][ismn][irow][icol];
 426                                   fSMIso[idet][ismn] += d1[idet][ismn][irow][icol];
 427                                   fCellIso[idet][ismn][irow][icol] += d1[idet][ismn][irow][icol];
 428                                   //fDetCount[idet]++;
 429                                   fSMCount[idet][ismn]++;
 430                                   fCellCount[idet][ismn][irow][icol]++;
 431
 432                                 }
 433                             }
 434                         }  // neigh cell cond.
 435                     }     // d>0 cond.
 436                 }
 437             }
 438         }
 439     }
 440
 441   for(Int_t idet=0; idet < kDet; idet++)
 442     {
 443       for(Int_t ismn = 0; ismn < kMaxSMN; ismn++)
 444         {
 445           for(Int_t irow = 0; irow < kMaxRow; irow++)
 446             {
 447               for(Int_t icol = 0; icol < kMaxCol; icol++)
 448                 {
 449                   if(fNhitCell[idet][ismn][irow][icol]>0)
 450                     {
 451                       fCountSm[idet][ismn]  += 1;
 452                       fTempnhit[idet][ismn] += fNhitCell[idet][ismn][irow][icol];
 453                       fTempnhitSq[idet][ismn] += fNhitCell[idet][ismn][irow][icol]
 454                         *fNhitCell[idet][ismn][irow][icol];
 455                     }
 456                 }
 457             }
 458         }
 459     }
 460
 461   if (numberofDDLs < kDDL) return kFALSE;
 462   return kTRUE;
 463
 464 }
 465 // ------------------------------------------------------------------------ //
 466 void AliPMDCalibGain::Analyse(TTree *gaintree, TTree *meantree)
 467 {
 468   // Calculates the mean
 469   Int_t   det, sm, row, col;
 470   Float_t gain;
 471   Float_t cellmean = 0.;
 472
 473   Float_t modmean[2][24];
 474
 475   for (Int_t idet=0; idet < 2; idet++)
 476     {
 477       for (Int_t ism = 0; ism < 24; ism++)
 478         {
 479           modmean[idet][ism] = 0.;
 480         }
 481     }
 482
 483   gaintree->Branch("det",&det,"det/I");
 484   gaintree->Branch("sm",&sm,"sm/I");
 485   gaintree->Branch("row",&row,"row/I");
 486   gaintree->Branch("col",&col,"col/I");
 487   gaintree->Branch("gain",&gain,"gain/F");
 488
 489   for(Int_t idet = 0; idet < kDet; idet++)
 490     {
 491       for(Int_t ism = 0; ism < kMaxSMN; ism++)
 492         {
 493           if (fSMCount[idet][ism] > 0)
 494             modmean[idet][ism] = fSMIso[idet][ism]/fSMCount[idet][ism];
 495           for(Int_t irow = 0; irow < kMaxRow; irow++)
 496             {
 497               for(Int_t icol = 0; icol < kMaxCol; icol++)
 498                 {
 499                   if (fCellCount[idet][ism][irow][icol] > 0.)
 500                     {
 501                       cellmean = fCellIso[idet][ism][irow][icol]/fCellCount[idet][ism][irow][icol];
 502                     }
 503                   det      = idet;
 504                   sm       = ism;
 505                   row      = irow;
 506                   col      = icol;
 507                   if (cellmean > 0.0 && fCellCount[idet][ism][irow][icol]>0.)
 508                     {
 509                       gain = cellmean/modmean[idet][ism];
 510                     }
 511                   else
 512                     {
 513                       gain = 0.;
 514                     }
 515                   //if(fCellCount[idet][ism][irow][icol]>0.) printf("CellCount =%f, gain= %f\n",fCellCount[idet][ism][irow][icol],gain);
 516                   gaintree->Fill();
 517                 }
 518             }
 519         }
 520     }
 521
 522   Float_t smmean;
 523
 524   // Writing each module mean value
 525   meantree->Branch("det",&det,"det/I");
 526   meantree->Branch("sm",&sm,"sm/I");
 527   meantree->Branch("smmean",&smmean,"row/F");
 528
 529   for(Int_t idet = 0; idet < kDet; idet++)
 530     {
 531       for (Int_t ism = 0; ism < kMaxSMN; ism++)
 532         {
 533           det    = idet;
 534           sm     = ism;
 535           smmean = modmean[idet][ism];
 536           meantree->Fill();
 537         }
 538     }
 539
 540 }
 541 // ------------------------------------------------------------------------ //
 542 void AliPMDCalibGain::FindHotCell(TTree *hottree, Float_t xvar)
 543 {
 544   // Calculates the mean
 545   Int_t   det, sm, row, col;
 546   Float_t flag;
 547   Float_t meannhit;
 548   Float_t meanSqnhit;
 549   Float_t sigmanhit,nhitcut;
 550
 551   //Float_t xvar = 5.;
 552
 553   hottree->Branch("det",&det,"det/I");
 554   hottree->Branch("sm",&sm,"sm/I");
 555   hottree->Branch("row",&row,"row/I");
 556   hottree->Branch("col",&col,"col/I");
 557   hottree->Branch("flag",&flag,"flag/F");
 558
 559   for(Int_t idet = 0; idet < kDet; idet++)
 560     {
 561       for(Int_t ism = 0; ism < kMaxSMN; ism++)
 562         {
 563           if (fCountSm[idet][ism]> 0)
 564             {
 565               meannhit   = fTempnhit[idet][ism]/fCountSm[idet][ism];
 566               meanSqnhit = fTempnhitSq[idet][ism]/fCountSm[idet][ism];
 567               sigmanhit  = sqrt(meanSqnhit-(meannhit*meannhit));
 568               nhitcut    = meannhit + xvar*sigmanhit;
 569
 570               for(Int_t irow = 0; irow < kMaxRow; irow++)
 571                 {
 572                   for(Int_t icol = 0; icol < kMaxCol; icol++)
 573                     {
 574                       det      = idet;
 575                       sm       = ism;
 576                       row      = irow;
 577                       col      = icol;
 578
 579                       if(fNhitCell[idet][ism][irow][icol] > nhitcut)
 580                         {
 581                           flag  = 1.0;
 582                         }
 583                       else
 584                         {
 585                           flag = 0.;
 586                         }
 587                       hottree->Fill();
 588                     }
 589
 590                 }
 591             }
 592         }
 593     }
 594 }
 595