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