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