PMD/AliPMDCalibrator.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 "TMath.h"
  27
  28 // --- Standard library ---
  29
  30 // --- AliRoot header files ---
  31 #include "AliLog.h"
  32 #include "AliRawReaderFile.h"
  33 #include "AliPMDCalibrator.h"
  34 #include "AliRawReaderDate.h"
  35 #include "AliPMDRawStream.h"
  36 #include "AliPMDCalibData.h"
  37 #include "AliPMDddldata.h"
  38 #include "AliCDBManager.h"
  39 #include "AliCDBId.h"
  40 #include "AliCDBMetaData.h"
  41 #include "AliDAQ.h"
  42
  43 //#include "AliPMDCleanNoise.h"
  44 //#include "AliPMDCleaner.h"
  45
  46 #include "AliPMDPedestal.h"
  47 #include "AliCDBManager.h"
  48 #include "AliCDBEntry.h"
  49
  50 ClassImp(AliPMDCalibrator)
  51
  52
  53 AliPMDCalibrator::AliPMDCalibrator():
  54   fCalibGain(new AliPMDCalibData()), fCalibPed(new AliPMDPedestal())
  55 {
  56   // Standard Constructor
  57   for(Int_t idet = 0; idet < kDet; idet++)
  58     {
  59       fHdetIso[idet] = NULL ;
  60       for(Int_t ismn = 0; ismn < kMaxSMN; ismn++)
  61         {
  62           fHsmIso[idet][ismn] = NULL ;
  63           for(Int_t jrow = 0; jrow < kMaxRow; jrow++)
  64             {
  65               for(Int_t kcol = 0; kcol < kMaxCol; kcol++)
  66                 {
  67                   fGainFact[idet][ismn][jrow][kcol] = 0.0;
  68                   fHadcIso[idet][ismn][jrow][kcol]  = NULL;
  69                 }
  70             }
  71         }
  72     }
  73 }
  74 // ------------------------------------------------------------------------ //
  75 AliPMDCalibrator::AliPMDCalibrator(const AliPMDCalibrator &pmdcalibrator):
  76   fCalibGain(new AliPMDCalibData()), fCalibPed(new AliPMDPedestal())
  77 {
  78   for(Int_t idet = 0; idet < 2; idet++)
  79     {
  80       fHdetIso[idet] = pmdcalibrator.fHdetIso[idet] ;
  81       for(Int_t ismn = 0; ismn < kMaxSMN; ismn++)
  82         {
  83           fHsmIso[idet][ismn] = pmdcalibrator.fHsmIso[idet][ismn] ;
  84           for(Int_t jrow = 0; jrow < kMaxRow; jrow++)
  85             {
  86               for(Int_t kcol = 0; kcol < kMaxCol; kcol++)
  87                 {
  88                   fGainFact[idet][ismn][jrow][kcol] = pmdcalibrator.fGainFact[idet][ismn][jrow][kcol];
  89                   fHadcIso[idet][ismn][jrow][kcol]  = pmdcalibrator.fHadcIso[idet][ismn][jrow][kcol];
  90                 }
  91             }
  92         }
  93     }
  94
  95 }
  96 // ------------------------------------------------------------------------ //
  97 AliPMDCalibrator &AliPMDCalibrator::operator=(const AliPMDCalibrator &pmdcalibrator)
  98 {
  99   if(this != &pmdcalibrator)
 100     {
 101       for(Int_t idet = 0; idet < kDet; idet++)
 102         {
 103           fHdetIso[idet] = pmdcalibrator.fHdetIso[idet] ;
 104           for(Int_t ismn = 0; ismn < kMaxSMN; ismn++)
 105             {
 106               fHsmIso[idet][ismn] = pmdcalibrator.fHsmIso[idet][ismn] ;
 107               for(Int_t jrow = 0; jrow < kMaxRow;jrow++)
 108                 {
 109                   for(Int_t kcol = 0; kcol < kMaxCol; kcol++)
 110                     {
 111                       fGainFact[idet][ismn][jrow][kcol] =
 112                         pmdcalibrator.fGainFact[idet][ismn][jrow][kcol];
 113                       fHadcIso[idet][ismn][jrow][kcol]  =
 114                         pmdcalibrator.fHadcIso[idet][ismn][jrow][kcol];
 115                     }
 116                 }
 117             }
 118         }
 119     }
 120   return *this;
 121 }
 122 // ------------------------------------------------------------------------ //
 123 AliPMDCalibrator::~AliPMDCalibrator()
 124 {
 125   // destructor
 126   for (Int_t i=0; i<kDet; i++)
 127     delete fHdetIso[i] ;
 128   for (Int_t i=0; i<kDet; i++)
 129     for (Int_t j=0; j<kMaxSMN; j++)
 130       delete fHsmIso[i][j] ;
 131   for (Int_t i=0; i<kDet; i++)
 132     for (Int_t j=0; j<kMaxSMN; j++)
 133       for (Int_t k=0; k<kMaxRow; k++)
 134         for (Int_t l=0; l<kMaxCol; l++)
 135           delete fHadcIso[i][j][k][l] ;
 136   delete fCalibGain;
 137   delete fCalibPed;
 138 }
 139 // ------------------------------------------------------------------------ //
 140
 141 void AliPMDCalibrator::Exec()
 142 {
 143   // reads parameters and does the calibration
 144   CalculateIsoCell() ;
 145 }
 146 // ------------------------------------------------------------------------ //
 147
 148 void AliPMDCalibrator::Init()
 149 {
 150   // intializes everything
 151   char hname[50];
 152   char hname2[50];
 153   char hname24[50];
 154   char hnameiso[120];
 155   char htitle1[120];
 156   char htitle2[120];
 157
 158   for(Int_t d = 0; d < kDet; d++) {
 159     snprintf(hname2,50,"Isolated cell adc for Det Plane %d",d);
 160     fHdetIso[d]= new TH1F(hname2,htitle2,100,0,1000);
 161     for(Int_t i1 = 0; i1 < kMaxSMN; i1++) {
 162       snprintf(hname,50,"det_%d_iso_sm_%2d",d,i1);
 163       snprintf(hname24,50,"det_%d_iso_sm_%2d",d,i1);
 164       fHsmIso[d][i1]= new TH1F(hname,hname24,100,0,1000);
 165       for(Int_t j1 = 0; j1 < kMaxRow; j1++) {
 166         for(Int_t k1 = 0; k1 < kMaxCol; k1++) {
 167           snprintf(hnameiso,120,"Isolated Cell ADC for det_%d_cell_sm%d_row%d_col%d"
 168                   ,d,i1,j1,k1);
 169           snprintf(htitle1,120,"Isolated Cell ADC for det_%d_cell_sm%d_row%d_col%d"
 170                   ,d,i1,j1,k1);
 171
 172           TObject *old=gDirectory->GetList()->FindObject(hnameiso);
 173           if (old) gDirectory->GetList()->Remove(old);
 174           fHadcIso[d][i1][j1][k1] = new TH1F(hnameiso,htitle1,100,0.,1000.);
 175         }
 176       }
 177     }
 178   }
 179
 180 }
 181
 182 // ------------------------------------------------------------------------ //
 183
 184 void AliPMDCalibrator::CalculateIsoCell()
 185 {
 186   // Calculates the ADC of isolated cell
 187
 188   TObjArray pmdddlcont;
 189
 190   const Int_t kCellNeighbour = 6;
 191
 192   Int_t neibx[6] = {1,0,-1,-1,0,1};
 193   Int_t neiby[6] = {0,1,1,0,-1,-1};
 194
 195   Int_t id1 = 0,jd1 = 0;        //neighbour row/col
 196   Int_t countisocell = 0 ;      //number of isilated cell
 197   Int_t isocount = 0;           //number of neighbours with 0 signal
 198   Int_t maxhit = 0;
 199   Int_t nhit[kDet][kMaxSMN];
 200   Int_t d1[kDet][kMaxSMN][kMaxRow][kMaxCol];
 201   Int_t nhitcell[kDet][kMaxSMN][kMaxRow][kMaxCol];
 202
 203   for(Int_t idet = 0; idet < kDet; idet++)
 204     {
 205       for(Int_t ismn = 0; ismn < kMaxSMN; ismn++)
 206         {
 207           nhit[idet][ismn] = 0;
 208           for(Int_t irow = 0; irow < kMaxRow; irow++)
 209             {
 210               for(Int_t icol = 0; icol < kMaxCol; icol++)
 211                 {
 212                   d1[idet][ismn][irow][icol] = 0;
 213                   nhitcell[idet][ismn][irow][icol] = 0;
 214                 }
 215             }
 216         }
 217     }
 218
 219   Float_t tempnhit1[kDet][kMaxSMN];
 220   Float_t tempnhit2[kDet][kMaxSMN];
 221   Float_t meannhit[kDet][kMaxSMN];
 222   Float_t meanSqnhit[kDet][kMaxSMN];
 223   Float_t sigmanhit[kDet][kMaxSMN];
 224   Float_t count[kDet][kMaxSMN];
 225   Float_t nhitcut[kDet][kMaxSMN];
 226
 227   for (Int_t idet = 0; idet < kDet; idet++)
 228   {
 229       for (Int_t ismn = 0; ismn < kMaxSMN; ismn++)
 230       {
 231           tempnhit1[idet][ismn]  = 0.;
 232           tempnhit2[idet][ismn]  = 0.;
 233           meannhit[idet][ismn]   = 0.;
 234           meanSqnhit[idet][ismn] = 0.;
 235           sigmanhit[idet][ismn]  = 0.;
 236           count[idet][ismn]      = 0.;
 237           nhitcut[idet][ismn]    = 0.;
 238       }
 239   }
 240
 241
 242   //accessing raw data
 243   AliRawReaderFile reader(".");
 244   AliPMDRawStream stream(&reader);
 245   while(reader.NextEvent())
 246     {
 247       // New PMD Reader is plugged in
 248         Int_t iddl = -1;
 249         while ((iddl = stream.DdlData(&pmdddlcont)) >=0) {
 250
 251             //reader.Select("PMD", iddl, iddl);
 252             //stream.DdlData(iddl,&pmdddlcont);
 253           Int_t ientries = pmdddlcont.GetEntries();
 254           for (Int_t ient = 0; ient < ientries; ient++)
 255             {
 256               AliPMDddldata *pmdddl =
 257                 (AliPMDddldata*)pmdddlcont.UncheckedAt(ient);
 258               Int_t idet = pmdddl->GetDetector();
 259               Int_t ismn = pmdddl->GetSMN();
 260               Int_t irow = pmdddl->GetRow();
 261               Int_t icol = pmdddl->GetColumn();
 262               Float_t isig1 = pmdddl->GetSignal();
 263               // Pedestal Subtraction
 264               Int_t   pedmeanrms =
 265                 fCalibPed->GetPedMeanRms(idet,ismn,irow,icol);
 266               Int_t   pedrms1    = (Int_t) pedmeanrms%100;
 267               Float_t pedrms     = (Float_t)pedrms1/10.;
 268               Float_t pedmean    = (Float_t) (pedmeanrms - pedrms1)/1000.0;
 269               Float_t isig = isig1 - (pedmean + 3.0*pedrms);
 270               if (isig>0)
 271                 {
 272                   d1[idet][ismn][irow][icol] = (Int_t)isig;
 273                   nhitcell[idet][ismn][irow][icol] += 1;
 274                 }
 275             }//ient loop
 276           pmdddlcont.Clear();
 277         }//iddl loop
 278       maxhit = 0;
 279       for(Int_t idet=0; idet < kDet; idet++)
 280         {
 281           for(Int_t ismn = 0; ismn < kMaxSMN; ismn++)
 282             {
 283               for(Int_t irow = 0; irow < kMaxRow; irow++)
 284                 {
 285                   for(Int_t icol = 0; icol < kMaxCol; icol++)
 286                     {
 287                       if(d1[idet][ismn][irow][icol] > 0)
 288                         {
 289                           isocount = 0;
 290                           for(Int_t ii = 0; ii < kCellNeighbour; ii++)
 291                             {
 292                               id1 = irow + neibx[ii];
 293                               jd1 = icol + neiby[ii];
 294                               if(d1[idet][ismn][id1][jd1] == 0)
 295                                 {
 296                                   isocount++;
 297                                   if(isocount == kCellNeighbour)
 298                                     {
 299                                       countisocell++;
 300                                       maxhit++;
 301                                       fHdetIso[idet]->
 302                                         Fill(d1[idet][ismn][irow][icol]);
 303                                       fHsmIso[idet][ismn]->
 304                                         Fill(d1[idet][ismn][irow][icol]);
 305                                       fHadcIso[idet][ismn][irow][icol]->
 306                                         Fill(d1[idet][ismn][irow][icol]);
 307                                     }
 308                                 }
 309                             }  // neigh cell cond.
 310                         }     // d>0 cond.
 311                     }
 312                 }
 313             }
 314         } //det
 315     }//event loop
 316
 317   // Mean and Sigma Calculations
 318   for(Int_t idet=0; idet < kDet; idet++){
 319       for(Int_t ismn = 0; ismn < kMaxSMN; ismn++){
 320           for(Int_t irow = 0; irow < kMaxRow; irow++){
 321               for(Int_t icol = 0; icol < kMaxCol; icol++){
 322                   if(nhitcell[idet][ismn][irow][icol]>0){
 323                       count[idet][ismn] += 1;
 324                       tempnhit1[idet][ismn] += nhitcell[idet][ismn][irow][icol];
 325                       tempnhit2[idet][ismn] += nhitcell[idet][ismn][irow][icol]
 326                           *nhitcell[idet][ismn][irow][icol];
 327                   }
 328               }
 329           }
 330       }
 331   }//det loop
 332
 333   //cout<<"nhit cell = "<<idet<<"  "<<ismn<<"  "
 334   //  <<irow<<"  "<<icol<<"  "<<nhitcell[idet][ismn][irow][icol]<<endl;
 335   //count[idet][ismn] += 1;
 336
 337   for(Int_t i=0; i < kDet; i++)
 338     {
 339       for(Int_t j=0; j < kMaxSMN; j++)
 340         {
 341           if(count[i][j] > 0.)
 342             {
 343                 meannhit[i][j]   = tempnhit1[i][j]/count[i][j];
 344                 meanSqnhit[i][j] = tempnhit2[i][j]/count[i][j];
 345                 sigmanhit[i][j]  = sqrt(meanSqnhit[i][j]-
 346                                       (meannhit[i][j]*meannhit[i][j]));
 347                 nhitcut[i][j]    = 3*sigmanhit[i][j] + meannhit[i][j];
 348             }
 349         }
 350     }
 351
 352
 353   Double_t histdetMean[kDet];
 354   Double_t histMean[kDet][kMaxSMN];
 355   Double_t isoMean[kDet][kMaxSMN][kMaxRow][kMaxCol];
 356   Double_t smNormFactor[kDet][kMaxSMN];
 357
 358   for(Int_t det1 = 0; det1 < kDet; det1++)
 359     {
 360       histdetMean[det1]= fHdetIso[det1]->GetMean();
 361       for(Int_t i1 = 0; i1 < kMaxSMN; i1++)
 362         {
 363           histMean[det1][i1]= fHsmIso[det1][i1]->GetMean();
 364           if(histMean[det1][i1]>0.0 && histdetMean[det1]>0.0)
 365             {
 366               smNormFactor[det1][i1]= histdetMean[det1]/histMean[det1][i1];
 367             }
 368           for(Int_t j1 = 0; j1 < kMaxRow; j1++)
 369             {
 370               for(Int_t k1 = 0; k1 < kMaxCol; k1++)
 371                 {
 372                  if(nhitcell[det1][i1][j1][k1]> nhitcut[det1][i1]) fGainFact[det1][i1][j1][k1]=-99.0;
 373                    if(nhitcell[det1][i1][j1][k1]< nhitcut[det1][i1])
 374                      {
 375                       isoMean[det1][i1][j1][k1]=fHadcIso[det1][i1][j1][k1]->
 376                         GetMean();
 377                       if(isoMean[det1][i1][j1][k1]>0.0 && histMean[det1][i1]>0.0)
 378                         {
 379                           fGainFact[det1][i1][j1][k1]=
 380                             isoMean[det1][i1][j1][k1]/(histMean[det1][i1]*
 381                                                      smNormFactor[det1][i1]);
 382
 383                         }
 384                     }
 385                    Float_t gain = fGainFact[det1][i1][j1][k1];
 386                    fCalibGain->SetGainFact(det1,i1,j1,k1,gain);
 387                 }
 388
 389
 390             }
 391         }
 392     }
 393 }//CalculateIsoCell()
 394
 395 // ------------------------------------------------------------------------ //
 396
 397 Bool_t AliPMDCalibrator::Store()
 398 {
 399   AliCDBManager *man = AliCDBManager::Instance();
 400   //man->SetDefaultStorage("local://$ALICE_ROOT/OCDB");
 401   if(!man->IsDefaultStorageSet()) return kFALSE;
 402   AliCDBId id("PMD/Calib/Gain",0,999999999);
 403   AliCDBMetaData md;
 404   md.SetBeamPeriod(0);
 405   md.SetComment("Test");
 406
 407   printf("\n\n\n fCalibData\n");
 408   //fCalibData->Print(0);
 409   //printf("\n\n\n fCalibData\n");
 410
 411   Bool_t result = man->Put(fCalibGain,id,&md);
 412
 413   return result;
 414 }
 415