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][kMaxSMN][kMaxRow][kMaxCol] = 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   fpw = fopen(tempFile,"r");
 249
 250   Int_t iddet = 0, issm = 0, irrow = 0, iccol = 0;
 251   Float_t smcount = 0., smiso = 0.;
 252   Float_t cellcount = 0., celliso = 0.;
 253
 254
 255   for (Int_t idet = 0; idet < kDet; idet++)
 256     {
 257       for (Int_t ism = 0; ism < kMaxSMN; ism++)
 258         {
 259           fscanf(fpw,"%d %d %f %f",&iddet,&issm,&smcount,&smiso);
 260           fSMCount[idet][ism] = smcount;
 261           fSMIso[idet][ism]   = smiso;
 262         }
 263     }
 264
 265   for (Int_t idet = 0; idet < kDet; idet++)
 266     {
 267       for (Int_t ism = 0; ism < kMaxSMN; ism++)
 268         {
 269           for (Int_t irow = 0; irow < kMaxRow; irow++)
 270             {
 271               for (Int_t icol = 0; icol < kMaxCol; icol++)
 272                 {
 273                   fscanf(fpw,"%d %d %d %d %f %f",&iddet,&issm,&irrow,&iccol,
 274                          &cellcount,&celliso);
 275                   fCellCount[idet][ism][irow][icol] = cellcount;
 276                   fCellIso[idet][ism][irow][icol]   = celliso;
 277                 }
 278             }
 279         }
 280     }
 281
 282   fclose(fpw);
 283
 284 }
 285 // ------------------------------------------------------------------------ //
 286 void AliPMDCalibGain::WriteTempFile(const Char_t *tempFile)
 287 {
 288   // Write the Temporary file if the required statics is not achieved
 289
 290
 291   /*
 292     Following variables to be written to a file
 293     fDetIso[idet] ;
 294     fSMIso[idet][ismn];
 295     fCellIso[idet][ismn][irow][icol];
 296
 297     fDetCount[idet];
 298     fSMCount[idet][ismn];
 299     fCellCount[idet][ismn][irow][icol];
 300   */
 301
 302
 303   fpw = fopen(tempFile,"w+");
 304
 305   for (Int_t idet = 0; idet < kDet; idet++)
 306     {
 307      // fprintf(fpw,"%d %f %f\n",idet,fDetCount[idet],fDetIso[idet]);
 308     }
 309
 310   for (Int_t idet = 0; idet < kDet; idet++)
 311     {
 312       for (Int_t ism = 0; ism < kMaxSMN; ism++)
 313         {
 314           fprintf(fpw,"%d %d %f %f\n",idet,ism, fSMCount[idet][ism],fSMIso[idet][ism]);
 315         }
 316     }
 317
 318   for (Int_t idet = 0; idet < kDet; idet++)
 319     {
 320       for (Int_t ism = 0; ism < kMaxSMN; ism++)
 321         {
 322           for (Int_t irow = 0; irow < kMaxRow; irow++)
 323             {
 324               for (Int_t icol = 0; icol < kMaxCol; icol++)
 325                 {
 326                   fprintf(fpw,"%d %d %d %d %f %f\n",idet,ism,irow,icol,
 327                           fCellCount[idet][ism][irow][icol],
 328                           fCellIso[idet][ism][irow][icol]);
 329                 }
 330             }
 331         }
 332     }
 333
 334   fclose(fpw);
 335
 336 }
 337
 338 // ------------------------------------------------------------------------ //
 339
 340 Bool_t AliPMDCalibGain::ProcessEvent(AliRawReader *rawReader, TObjArray *pmdddlcont)
 341 {
 342   // Calculates the ADC of isolated cell
 343
 344   const Int_t kDDL           = AliDAQ::NumberOfDdls("PMD");
 345   const Int_t kCellNeighbour = 6;
 346   Int_t neibx[6] = {1,0,-1,-1,0,1};
 347   Int_t neiby[6] = {0,1,1,0,-1,-1};
 348
 349   Int_t id1 = 0,jd1 = 0;  //neighbour row/col
 350   Int_t isocount = 0;     //number of neighbours with 0 signal
 351
 352   Float_t d1[kDet][kMaxSMN][kMaxRow][kMaxCol];
 353
 354   for(Int_t idet = 0; idet < kDet; idet++)
 355     {
 356       for(Int_t ismn = 0; ismn < kMaxSMN; ismn++)
 357         {
 358           for(Int_t irow = 0; irow < kMaxRow; irow++)
 359             {
 360               for(Int_t icol = 0; icol < kMaxCol; icol++)
 361                 {
 362                   d1[idet][ismn][irow][icol] = 0.;
 363                 }
 364             }
 365         }
 366     }
 367
 368   AliPMDRawStream rawStream(rawReader);
 369
 370   Int_t iddl = -1;
 371
 372   Int_t numberofDDLs = 0;
 373
 374   while ((iddl = rawStream.DdlData(pmdddlcont)) >=0) {
 375     numberofDDLs++;
 376
 377     Int_t ientries = pmdddlcont->GetEntries();
 378
 379     for (Int_t ient = 0; ient < ientries; ient++)
 380       {
 381         AliPMDddldata *pmdddl = (AliPMDddldata*)pmdddlcont->UncheckedAt(ient);
 382
 383         Int_t idet = pmdddl->GetDetector();
 384         Int_t ismn = pmdddl->GetSMN();
 385         Int_t mcm = pmdddl->GetMCM();
 386         //Int_t ichno = pmdddl->GetChannel();
 387         Int_t irow = pmdddl->GetRow();
 388         Int_t icol = pmdddl->GetColumn();
 389         Int_t isig = pmdddl->GetSignal();
 390
 391         // This is the protection not to crash the code
 392
 393         if(mcm == 0) continue;
 394         if (irow < 0 || icol < 0 || irow > 47 || icol > 95) continue;
 395
 396         // Pedestal subtraction
 397         if(fHotFlag[idet][ismn][irow][icol] == 1.0) isig = 0;
 398
 399         if (isig>0)
 400           {
 401             d1[idet][ismn][irow][icol] =
 402               (Float_t) isig - fPedMeanRMS[idet][ismn][irow][icol];
 403
 404             //printf("Signal_ped_subtracted=%f, pedestal=%f\n",d1[idet][ismn][irow][icol]),fPedMeanRMS[idet][ismn][irow][icol];
 405
 406             fNhitCell[idet][ismn][irow][icol]++;     //cell hit frequency
 407           }
 408
 409
 410       }
 411     pmdddlcont->Delete();
 412     }
 413
 414   for(Int_t idet=0; idet < kDet; idet++)
 415     {
 416       for(Int_t ismn = 0; ismn < kMaxSMN; ismn++)
 417         {
 418           for(Int_t irow = 0; irow < kMaxRow; irow++)
 419             {
 420               for(Int_t icol = 0; icol < kMaxCol; icol++)
 421                 {
 422                   if(d1[idet][ismn][irow][icol] > 0)
 423                     {
 424                       isocount = 0;
 425                       for(Int_t ii = 0; ii < kCellNeighbour; ii++)
 426                         {
 427                           id1 = irow + neibx[ii];
 428                           jd1 = icol + neiby[ii];
 429                           if (id1 < 0) id1 = 0;
 430                           if (id1 > kMaxRow-1) id1 = kMaxRow - 1;
 431                           if (jd1 < 0) jd1 = 0;
 432                           if (jd1 > kMaxCol-1) jd1 = kMaxCol - 1;
 433                           if(d1[idet][ismn][id1][jd1] == 0)
 434                             {
 435                               isocount++;
 436                               if(isocount == kCellNeighbour)
 437                                 {
 438                                   //fDetIso[idet] += d1[idet][ismn][irow][icol];
 439                                   fSMIso[idet][ismn] += d1[idet][ismn][irow][icol];
 440                                   fCellIso[idet][ismn][irow][icol] += d1[idet][ismn][irow][icol];
 441                                   //fDetCount[idet]++;
 442                                   fSMCount[idet][ismn]++;
 443                                   fCellCount[idet][ismn][irow][icol]++;
 444
 445                                 }
 446                             }
 447                         }  // neigh cell cond.
 448                     }     // d>0 cond.
 449                 }
 450             }
 451         }
 452     }
 453
 454   for(Int_t idet=0; idet < kDet; idet++)
 455     {
 456       for(Int_t ismn = 0; ismn < kMaxSMN; ismn++)
 457         {
 458           for(Int_t irow = 0; irow < kMaxRow; irow++)
 459             {
 460               for(Int_t icol = 0; icol < kMaxCol; icol++)
 461                 {
 462                   if(fNhitCell[idet][ismn][irow][icol]>0)
 463                     {
 464                       fCountSm[idet][ismn]  += 1;
 465                       fTempnhit[idet][ismn] += fNhitCell[idet][ismn][irow][icol];
 466                       fTempnhitSq[idet][ismn] += fNhitCell[idet][ismn][irow][icol]
 467                         *fNhitCell[idet][ismn][irow][icol];
 468                     }
 469                 }
 470             }
 471         }
 472     }
 473
 474   if (numberofDDLs < kDDL) return kFALSE;
 475   return kTRUE;
 476
 477 }
 478 // ------------------------------------------------------------------------ //
 479 void AliPMDCalibGain::Analyse(TTree *gaintree, TTree *meantree)
 480 {
 481   // Calculates the mean
 482   Int_t   det = 0, sm = 0, row = 0, col = 0;
 483   Float_t gain = 0.;
 484   Float_t cellmean = 0.;
 485
 486   Float_t modmean[2][24];
 487
 488   for (Int_t idet=0; idet < 2; idet++)
 489     {
 490       for (Int_t ism = 0; ism < 24; ism++)
 491         {
 492           modmean[idet][ism] = 0.;
 493         }
 494     }
 495
 496   gaintree->Branch("det",&det,"det/I");
 497   gaintree->Branch("sm",&sm,"sm/I");
 498   gaintree->Branch("row",&row,"row/I");
 499   gaintree->Branch("col",&col,"col/I");
 500   gaintree->Branch("gain",&gain,"gain/F");
 501
 502   for(Int_t idet = 0; idet < kDet; idet++)
 503     {
 504       for(Int_t ism = 0; ism < kMaxSMN; ism++)
 505         {
 506           if (fSMCount[idet][ism] > 0)
 507             modmean[idet][ism] = fSMIso[idet][ism]/fSMCount[idet][ism];
 508           for(Int_t irow = 0; irow < kMaxRow; irow++)
 509             {
 510               for(Int_t icol = 0; icol < kMaxCol; icol++)
 511                 {
 512                   if (fCellCount[idet][ism][irow][icol] > 0.)
 513                     {
 514                       cellmean = fCellIso[idet][ism][irow][icol]/fCellCount[idet][ism][irow][icol];
 515                     }
 516                   det      = idet;
 517                   sm       = ism;
 518                   row      = irow;
 519                   col      = icol;
 520                   if (cellmean > 0.0 && fCellCount[idet][ism][irow][icol]>0.)
 521                     {
 522                       gain = cellmean/modmean[idet][ism];
 523                     }
 524                   else
 525                     {
 526                       gain = 0.;
 527                     }
 528                   //if(fCellCount[idet][ism][irow][icol]>0.) printf("CellCount =%f, gain= %f\n",fCellCount[idet][ism][irow][icol],gain);
 529                   gaintree->Fill();
 530                 }
 531             }
 532         }
 533     }
 534
 535   Float_t smmean;
 536
 537   // Writing each module mean value
 538   meantree->Branch("det",&det,"det/I");
 539   meantree->Branch("sm",&sm,"sm/I");
 540   meantree->Branch("smmean",&smmean,"row/F");
 541
 542   for(Int_t idet = 0; idet < kDet; idet++)
 543     {
 544       for (Int_t ism = 0; ism < kMaxSMN; ism++)
 545         {
 546           det    = idet;
 547           sm     = ism;
 548           smmean = modmean[idet][ism];
 549           meantree->Fill();
 550         }
 551     }
 552
 553 }
 554 // ------------------------------------------------------------------------ //
 555 void AliPMDCalibGain::FindHotCell(TTree *hottree, Float_t xvar)
 556 {
 557   // Calculates the mean
 558   Int_t   det = 0, sm = 0, row = 0, col = 0;
 559   Float_t flag = 0.;
 560   Float_t meannhit = 0.;
 561   Float_t meanSqnhit = 0.;
 562   Float_t sigmanhit = 0.,nhitcut = 0.;
 563
 564   //Float_t xvar = 5.;
 565
 566   hottree->Branch("det",&det,"det/I");
 567   hottree->Branch("sm",&sm,"sm/I");
 568   hottree->Branch("row",&row,"row/I");
 569   hottree->Branch("col",&col,"col/I");
 570   hottree->Branch("flag",&flag,"flag/F");
 571
 572   for(Int_t idet = 0; idet < kDet; idet++)
 573     {
 574       for(Int_t ism = 0; ism < kMaxSMN; ism++)
 575         {
 576           if (fCountSm[idet][ism]> 0)
 577             {
 578               meannhit   = fTempnhit[idet][ism]/fCountSm[idet][ism];
 579               meanSqnhit = fTempnhitSq[idet][ism]/fCountSm[idet][ism];
 580               sigmanhit  = sqrt(meanSqnhit-(meannhit*meannhit));
 581               nhitcut    = meannhit + xvar*sigmanhit;
 582
 583               for(Int_t irow = 0; irow < kMaxRow; irow++)
 584                 {
 585                   for(Int_t icol = 0; icol < kMaxCol; icol++)
 586                     {
 587                       det      = idet;
 588                       sm       = ism;
 589                       row      = irow;
 590                       col      = icol;
 591
 592                       if(fNhitCell[idet][ism][irow][icol] > nhitcut)
 593                         {
 594                           flag  = 1.0;
 595                         }
 596                       else
 597                         {
 598                           flag = 0.;
 599                         }
 600                       hottree->Fill();
 601                     }
 602
 603                 }
 604             }
 605         }
 606     }
 607 }
 608