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);//sid
 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   const Int_t kDDL           = AliDAQ::NumberOfDdls("PMD");
 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       for (Int_t iddl = 0; iddl < kDDL; iddl++)
 242         {
 243           reader.Select("PMD", iddl, iddl);
 244           stream.DdlData(iddl,&pmdddlcont);
 245           Int_t ientries = pmdddlcont.GetEntries();
 246           for (Int_t ient = 0; ient < ientries; ient++)
 247             {
 248               AliPMDddldata *pmdddl =
 249                 (AliPMDddldata*)pmdddlcont.UncheckedAt(ient);
 250               Int_t idet = pmdddl->GetDetector();
 251               Int_t ismn = pmdddl->GetSMN();
 252               Int_t irow = pmdddl->GetRow();
 253               Int_t icol = pmdddl->GetColumn();
 254               Float_t isig1 = pmdddl->GetSignal();
 255               // Pedestal Subtraction
 256               Int_t   pedmeanrms =
 257                 fCalibPed->GetPedMeanRms(idet,ismn,irow,icol);
 258               Int_t   pedrms1    = (Int_t) pedmeanrms%1000;
 259               Float_t pedrms     = (Float_t)pedrms1/10.;
 260               Float_t pedmean    = (Float_t) (pedmeanrms - pedrms1)/1000.0;
 261               Float_t isig = isig1 - (pedmean + 3.0*pedrms);
 262               if (isig>0)
 263                 {
 264                   d1[idet][ismn][irow][icol] = (Int_t)isig;
 265                   nhitcell[idet][ismn][irow][icol] += 1; //sid
 266                 }
 267             }//ient loop
 268           pmdddlcont.Clear();
 269         }//iddl loop
 270       maxhit = 0;
 271       for(Int_t idet=0; idet < kDet; idet++)
 272         {
 273           for(Int_t ismn = 0; ismn < kMaxSMN; ismn++)
 274             {
 275               for(Int_t irow = 0; irow < kMaxRow; irow++)
 276                 {
 277                   for(Int_t icol = 0; icol < kMaxCol; icol++)
 278                     {
 279                       if(d1[idet][ismn][irow][icol] > 0)
 280                         {
 281                           isocount = 0;
 282                           for(Int_t ii = 0; ii < kCellNeighbour; ii++)
 283                             {
 284                               id1 = irow + neibx[ii];
 285                               jd1 = icol + neiby[ii];
 286                               if(d1[idet][ismn][id1][jd1] == 0)
 287                                 {
 288                                   isocount++;
 289                                   if(isocount == kCellNeighbour)
 290                                     {
 291                                       countisocell++;
 292                                       maxhit++;
 293                                       fHdetIso[idet]->
 294                                         Fill(d1[idet][ismn][irow][icol]);
 295                                       fHsmIso[idet][ismn]->
 296                                         Fill(d1[idet][ismn][irow][icol]);
 297                                       fHadcIso[idet][ismn][irow][icol]->
 298                                         Fill(d1[idet][ismn][irow][icol]);
 299                                     }
 300                                 }
 301                             }  // neigh cell cond.
 302                         }     // d>0 cond.
 303                     }
 304                 }
 305             }
 306         } //det
 307     }//event loop
 308
 309   // Mean and Sigma Calculations
 310   for(Int_t idet=0; idet < kDet; idet++){
 311       for(Int_t ismn = 0; ismn < kMaxSMN; ismn++){
 312           for(Int_t irow = 0; irow < kMaxRow; irow++){
 313               for(Int_t icol = 0; icol < kMaxCol; icol++){
 314                   if(nhitcell[idet][ismn][irow][icol]>0){
 315                       count[idet][ismn] += 1;
 316                       tempnhit1[idet][ismn] += nhitcell[idet][ismn][irow][icol];
 317                       tempnhit2[idet][ismn] += nhitcell[idet][ismn][irow][icol]
 318                           *nhitcell[idet][ismn][irow][icol];
 319                   }
 320               }
 321           }
 322       }
 323   }//det loop
 324
 325   //cout<<"nhit cell = "<<idet<<"  "<<ismn<<"  "
 326   //  <<irow<<"  "<<icol<<"  "<<nhitcell[idet][ismn][irow][icol]<<endl;
 327   //count[idet][ismn] += 1;
 328
 329   for(Int_t i=0; i < kDet; i++)
 330     {
 331       for(Int_t j=0; j < kMaxSMN; j++)
 332         {
 333           if(count[i][j] > 0.)
 334             {
 335                 meannhit[i][j]   = tempnhit1[i][j]/count[i][j];
 336                 meanSqnhit[i][j] = tempnhit2[i][j]/count[i][j];
 337                 sigmanhit[i][j]  = sqrt(meanSqnhit[i][j]-
 338                                       (meannhit[i][j]*meannhit[i][j]));
 339                 nhitcut[i][j]    = 3*sigmanhit[i][j] + meannhit[i][j];
 340             }
 341         }
 342     }
 343
 344
 345   Double_t histdetMean[kDet];
 346   Double_t histMean[kDet][kMaxSMN];
 347   Double_t isoMean[kDet][kMaxSMN][kMaxRow][kMaxCol];
 348   Double_t smNormFactor[kDet][kMaxSMN];
 349
 350   for(Int_t det1 = 0; det1 < kDet; det1++)
 351     {
 352       histdetMean[det1]= fHdetIso[det1]->GetMean();
 353       for(Int_t i1 = 0; i1 < kMaxSMN; i1++)
 354         {
 355           histMean[det1][i1]= fHsmIso[det1][i1]->GetMean();
 356           if(histMean[det1][i1]>0.0 && histdetMean[det1]>0.0)
 357             {
 358               smNormFactor[det1][i1]= histdetMean[det1]/histMean[det1][i1];
 359             }
 360           for(Int_t j1 = 0; j1 < kMaxRow; j1++)
 361             {
 362               for(Int_t k1 = 0; k1 < kMaxCol; k1++)
 363                 {
 364                   if(nhitcell[det1][i1][j1][k1]< nhitcut[det1][i1])//sid
 365                     {
 366                       isoMean[det1][i1][j1][k1]=fHadcIso[det1][i1][j1][k1]->
 367                         GetMean();
 368                       if(isoMean[det1][i1][j1][k1]>0.0 && histMean[det1][i1]>0.0)
 369                         {
 370                           fGainFact[det1][i1][j1][k1]=
 371                             isoMean[det1][i1][j1][k1]/(histMean[det1][i1]*
 372                                                      smNormFactor[det1][i1]);
 373                         }
 374                     }
 375                 }
 376             }
 377         }
 378     }
 379 }//CalculateIsoCell()
 380
 381 // ------------------------------------------------------------------------ //
 382
 383 Bool_t AliPMDCalibrator::Store()
 384 {
 385   AliCDBManager *man = AliCDBManager::Instance();
 386   //man->SetDefaultStorage("local://$ALICE_ROOT");
 387   if(!man->IsDefaultStorageSet()) return kFALSE;
 388   AliCDBId id("PMD/Calib/Gain",0,0);
 389   AliCDBMetaData md;
 390   md.SetBeamPeriod(0);
 391   md.SetComment("Test");
 392
 393   printf("\n\n\n fCalibData\n");
 394   //fCalibData->Print(0);
 395   //printf("\n\n\n fCalibData\n");
 396
 397   Bool_t result = man->Put(fCalibGain,id,&md);
 398
 399   return result;
 400 }
 401