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   if(fHdetIso) delete fHdetIso ;
 127   if(fHsmIso)  delete fHsmIso ;
 128   if(fHadcIso) delete fHadcIso ;
 129   delete fCalibGain;
 130   delete fCalibPed;
 131 }
 132 // ------------------------------------------------------------------------ //
 133
 134 void AliPMDCalibrator::Exec()
 135 {
 136   // reads parameters and does the calibration
 137   CalculateIsoCell() ;
 138 }
 139 // ------------------------------------------------------------------------ //
 140
 141 void AliPMDCalibrator::Init()
 142 {
 143   // intializes everything
 144   char hname2[kDet];
 145   char htitle2[2];
 146   char hname[kMaxSMN];
 147   char hname24[kMaxSMN];
 148   char hnameiso[120];
 149   char htitle1[120];
 150
 151   for(Int_t d = 0; d < kDet; d++) {
 152     sprintf(hname2,"Isolated cell adc for Det Plane %d",d);
 153     fHdetIso[d]= new TH1F(hname2,htitle2,100,0,1000);
 154     for(Int_t i1 = 0; i1 < kMaxSMN; i1++) {
 155       sprintf(hname,"det_%d_iso_sm_%2d",d,i1);
 156       sprintf(hname24,"det_%d_iso_sm_%2d",d,i1);
 157       fHsmIso[d][i1]= new TH1F(hname,hname24,100,0,1000);
 158       for(Int_t j1 = 0; j1 < kMaxRow; j1++) {
 159         for(Int_t k1 = 0; k1 < kMaxCol; k1++) {
 160           sprintf(hnameiso,"Isolated Cell ADC for det_%d_cell_sm%d_row%d_col%d"
 161                   ,d,i1,j1,k1);
 162           sprintf(htitle1,"Isolated Cell ADC for det_%d_cell_sm%d_row%d_col%d"
 163                   ,d,i1,j1,k1);
 164
 165           TObject *old=gDirectory->GetList()->FindObject(hnameiso);
 166           if (old) gDirectory->GetList()->Remove(old);
 167           fHadcIso[d][i1][j1][k1] = new TH1F(hnameiso,htitle1,100,0.,1000.);
 168         }
 169       }
 170     }
 171   }
 172
 173 }
 174
 175 // ------------------------------------------------------------------------ //
 176
 177 void AliPMDCalibrator::CalculateIsoCell()
 178 {
 179   // Calculates the ADC of isolated cell
 180
 181   TObjArray pmdddlcont;
 182
 183   const Int_t kCellNeighbour = 6;
 184
 185   Int_t neibx[6] = {1,0,-1,-1,0,1};
 186   Int_t neiby[6] = {0,1,1,0,-1,-1};
 187
 188   Int_t id1,jd1;            //neighbour row/col
 189   Int_t countisocell = 0 ;  //number of isilated cell
 190   Int_t isocount;           //number of neighbours with 0 signal
 191   Int_t d1[kDet][kMaxSMN][kMaxRow][kMaxCol];
 192   Int_t maxhit;
 193   Int_t nhit[kDet][kMaxSMN];
 194   Int_t nhitcell[kDet][kMaxSMN][kMaxRow][kMaxCol];
 195
 196   for(Int_t idet = 0; idet < kDet; idet++)
 197     {
 198       for(Int_t ismn = 0; ismn < kMaxSMN; ismn++)
 199         {
 200           nhit[idet][ismn] = 0;
 201           for(Int_t irow = 0; irow < kMaxRow; irow++)
 202             {
 203               for(Int_t icol = 0; icol < kMaxCol; icol++)
 204                 {
 205                   d1[idet][ismn][irow][icol] = 0;
 206                   nhitcell[idet][ismn][irow][icol] = 0;
 207                 }
 208             }
 209         }
 210     }
 211
 212   Float_t tempnhit1[kDet][kMaxSMN];
 213   Float_t tempnhit2[kDet][kMaxSMN];
 214   Float_t meannhit[kDet][kMaxSMN];
 215   Float_t meanSqnhit[kDet][kMaxSMN];
 216   Float_t sigmanhit[kDet][kMaxSMN];
 217   Float_t count[kDet][kMaxSMN];
 218   Float_t nhitcut[kDet][kMaxSMN];
 219
 220   for (Int_t idet = 0; idet < kDet; idet++)
 221   {
 222       for (Int_t ismn = 0; ismn < kMaxSMN; ismn++)
 223       {
 224           tempnhit1[idet][ismn]  = 0.;
 225           tempnhit2[idet][ismn]  = 0.;
 226           meannhit[idet][ismn]   = 0.;
 227           meanSqnhit[idet][ismn] = 0.;
 228           sigmanhit[idet][ismn]  = 0.;
 229           count[idet][ismn]      = 0.;
 230           nhitcut[idet][ismn]    = 0.;
 231       }
 232   }
 233
 234
 235   //accessing raw data
 236   AliRawReaderFile reader(".");
 237   AliPMDRawStream stream(&reader);
 238   while(reader.NextEvent())
 239     {
 240       // New PMD Reader is plugged in
 241         Int_t iddl = -1;
 242         while ((iddl = stream.DdlData(&pmdddlcont)) >=0) {
 243
 244             //reader.Select("PMD", iddl, iddl);
 245             //stream.DdlData(iddl,&pmdddlcont);
 246           Int_t ientries = pmdddlcont.GetEntries();
 247           for (Int_t ient = 0; ient < ientries; ient++)
 248             {
 249               AliPMDddldata *pmdddl =
 250                 (AliPMDddldata*)pmdddlcont.UncheckedAt(ient);
 251               Int_t idet = pmdddl->GetDetector();
 252               Int_t ismn = pmdddl->GetSMN();
 253               Int_t irow = pmdddl->GetRow();
 254               Int_t icol = pmdddl->GetColumn();
 255               Float_t isig1 = pmdddl->GetSignal();
 256               // Pedestal Subtraction
 257               Int_t   pedmeanrms =
 258                 fCalibPed->GetPedMeanRms(idet,ismn,irow,icol);
 259               Int_t   pedrms1    = (Int_t) pedmeanrms%100;
 260               Float_t pedrms     = (Float_t)pedrms1/10.;
 261               Float_t pedmean    = (Float_t) (pedmeanrms - pedrms1)/1000.0;
 262               Float_t isig = isig1 - (pedmean + 3.0*pedrms);
 263               if (isig>0)
 264                 {
 265                   d1[idet][ismn][irow][icol] = (Int_t)isig;
 266                   nhitcell[idet][ismn][irow][icol] += 1;
 267                 }
 268             }//ient loop
 269           pmdddlcont.Clear();
 270         }//iddl loop
 271       maxhit = 0;
 272       for(Int_t idet=0; idet < kDet; idet++)
 273         {
 274           for(Int_t ismn = 0; ismn < kMaxSMN; ismn++)
 275             {
 276               for(Int_t irow = 0; irow < kMaxRow; irow++)
 277                 {
 278                   for(Int_t icol = 0; icol < kMaxCol; icol++)
 279                     {
 280                       if(d1[idet][ismn][irow][icol] > 0)
 281                         {
 282                           isocount = 0;
 283                           for(Int_t ii = 0; ii < kCellNeighbour; ii++)
 284                             {
 285                               id1 = irow + neibx[ii];
 286                               jd1 = icol + neiby[ii];
 287                               if(d1[idet][ismn][id1][jd1] == 0)
 288                                 {
 289                                   isocount++;
 290                                   if(isocount == kCellNeighbour)
 291                                     {
 292                                       countisocell++;
 293                                       maxhit++;
 294                                       fHdetIso[idet]->
 295                                         Fill(d1[idet][ismn][irow][icol]);
 296                                       fHsmIso[idet][ismn]->
 297                                         Fill(d1[idet][ismn][irow][icol]);
 298                                       fHadcIso[idet][ismn][irow][icol]->
 299                                         Fill(d1[idet][ismn][irow][icol]);
 300                                     }
 301                                 }
 302                             }  // neigh cell cond.
 303                         }     // d>0 cond.
 304                     }
 305                 }
 306             }
 307         } //det
 308     }//event loop
 309
 310   // Mean and Sigma Calculations
 311   for(Int_t idet=0; idet < kDet; idet++){
 312       for(Int_t ismn = 0; ismn < kMaxSMN; ismn++){
 313           for(Int_t irow = 0; irow < kMaxRow; irow++){
 314               for(Int_t icol = 0; icol < kMaxCol; icol++){
 315                   if(nhitcell[idet][ismn][irow][icol]>0){
 316                       count[idet][ismn] += 1;
 317                       tempnhit1[idet][ismn] += nhitcell[idet][ismn][irow][icol];
 318                       tempnhit2[idet][ismn] += nhitcell[idet][ismn][irow][icol]
 319                           *nhitcell[idet][ismn][irow][icol];
 320                   }
 321               }
 322           }
 323       }
 324   }//det loop
 325
 326   //cout<<"nhit cell = "<<idet<<"  "<<ismn<<"  "
 327   //  <<irow<<"  "<<icol<<"  "<<nhitcell[idet][ismn][irow][icol]<<endl;
 328   //count[idet][ismn] += 1;
 329
 330   for(Int_t i=0; i < kDet; i++)
 331     {
 332       for(Int_t j=0; j < kMaxSMN; j++)
 333         {
 334           if(count[i][j] > 0.)
 335             {
 336                 meannhit[i][j]   = tempnhit1[i][j]/count[i][j];
 337                 meanSqnhit[i][j] = tempnhit2[i][j]/count[i][j];
 338                 sigmanhit[i][j]  = sqrt(meanSqnhit[i][j]-
 339                                       (meannhit[i][j]*meannhit[i][j]));
 340                 nhitcut[i][j]    = 3*sigmanhit[i][j] + meannhit[i][j];
 341             }
 342         }
 343     }
 344
 345
 346   Double_t histdetMean[kDet];
 347   Double_t histMean[kDet][kMaxSMN];
 348   Double_t isoMean[kDet][kMaxSMN][kMaxRow][kMaxCol];
 349   Double_t smNormFactor[kDet][kMaxSMN];
 350
 351   for(Int_t det1 = 0; det1 < kDet; det1++)
 352     {
 353       histdetMean[det1]= fHdetIso[det1]->GetMean();
 354       for(Int_t i1 = 0; i1 < kMaxSMN; i1++)
 355         {
 356           histMean[det1][i1]= fHsmIso[det1][i1]->GetMean();
 357           if(histMean[det1][i1]>0.0 && histdetMean[det1]>0.0)
 358             {
 359               smNormFactor[det1][i1]= histdetMean[det1]/histMean[det1][i1];
 360             }
 361           for(Int_t j1 = 0; j1 < kMaxRow; j1++)
 362             {
 363               for(Int_t k1 = 0; k1 < kMaxCol; k1++)
 364                 {
 365                  if(nhitcell[det1][i1][j1][k1]> nhitcut[det1][i1]) fGainFact[det1][i1][j1][k1]=-99.0;
 366                    if(nhitcell[det1][i1][j1][k1]< nhitcut[det1][i1])
 367                      {
 368                       isoMean[det1][i1][j1][k1]=fHadcIso[det1][i1][j1][k1]->
 369                         GetMean();
 370                       if(isoMean[det1][i1][j1][k1]>0.0 && histMean[det1][i1]>0.0)
 371                         {
 372                           fGainFact[det1][i1][j1][k1]=
 373                             isoMean[det1][i1][j1][k1]/(histMean[det1][i1]*
 374                                                      smNormFactor[det1][i1]);
 375
 376                         }
 377                     }
 378                    Float_t gain = fGainFact[det1][i1][j1][k1];
 379                    fCalibGain->SetGainFact(det1,i1,j1,k1,gain);
 380                 }
 381
 382
 383             }
 384         }
 385     }
 386 }//CalculateIsoCell()
 387
 388 // ------------------------------------------------------------------------ //
 389
 390 Bool_t AliPMDCalibrator::Store()
 391 {
 392   AliCDBManager *man = AliCDBManager::Instance();
 393   //man->SetDefaultStorage("local://$ALICE_ROOT/OCDB");
 394   if(!man->IsDefaultStorageSet()) return kFALSE;
 395   AliCDBId id("PMD/Calib/Gain",0,999999999);
 396   AliCDBMetaData md;
 397   md.SetBeamPeriod(0);
 398   md.SetComment("Test");
 399
 400   printf("\n\n\n fCalibData\n");
 401   //fCalibData->Print(0);
 402   //printf("\n\n\n fCalibData\n");
 403
 404   Bool_t result = man->Put(fCalibGain,id,&md);
 405
 406   return result;
 407 }
 408