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 "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 "TTree.h"
  27
  28 // --- Standard library ---
  29
  30 // --- AliRoot header files ---
  31 #include "AliDAQ.h"
  32 #include "AliLog.h"
  33 #include "AliRawReader.h"
  34 #include "AliPMDRawStream.h"
  35 #include "AliPMDddldata.h"
  36 #include "AliPMDCalibGain.h"
  37
  38 ClassImp(AliPMDCalibGain)
  39
  40 AliPMDCalibGain::AliPMDCalibGain():
  41   TObject(),
  42   fpw(NULL)
  43 {
  44   // Standard Constructor
  45     for(Int_t idet = 0; idet < kDet; idet++)
  46     {
  47      fDetCount[kDet] =0.;
  48         for(Int_t ismn = 0; ismn < kMaxSMN; ismn++)
  49         {
  50             fSMIso[idet][ismn]   = 0.;
  51             fSMCount[idet][ismn] = 0.;
  52             for(Int_t jrow = 0; jrow < kMaxRow; jrow++)
  53             {
  54                 for(Int_t kcol = 0; kcol < kMaxCol; kcol++)
  55                 {
  56                     fCellIso[idet][ismn][jrow][kcol]   = 0.;
  57                     fCellCount[idet][ismn][jrow][kcol] = 0.;
  58                     fPedMeanRMS[idet][ismn][jrow][kcol] = 0.;
  59
  60                 }
  61             }
  62         }
  63     }
  64
  65
  66
  67 }
  68 // ------------------------------------------------------------------------ //
  69 AliPMDCalibGain::AliPMDCalibGain(const AliPMDCalibGain &pmdcalibgain):
  70   TObject(pmdcalibgain),
  71   fpw(NULL)
  72 {
  73     for(Int_t idet = 0; idet < kDet; idet++)
  74     {
  75      fDetCount[idet] = pmdcalibgain.fDetCount[idet];
  76      fDetIso[idet] = pmdcalibgain.fDetIso[idet];
  77         for(Int_t ismn = 0; ismn < kMaxSMN; ismn++)
  78         {
  79             fSMIso[idet][ismn] = pmdcalibgain.fSMIso[idet][ismn] ;
  80             fSMCount[idet][ismn] = pmdcalibgain.fSMCount[idet][ismn] ;
  81             for(Int_t jrow = 0; jrow < kMaxRow; jrow++)
  82             {
  83                 for(Int_t kcol = 0; kcol < kMaxCol; kcol++)
  84                 {
  85                   fCellIso[idet][ismn][jrow][kcol]    =
  86                     pmdcalibgain.fCellIso[idet][ismn][jrow][kcol];
  87                   fCellCount[idet][ismn][jrow][kcol]  =
  88                     pmdcalibgain.fCellCount[idet][ismn][jrow][kcol];
  89                   fPedMeanRMS[idet][ismn][jrow][kcol] =
  90                     pmdcalibgain.fPedMeanRMS[idet][ismn][jrow][kcol];
  91
  92                 }
  93             }
  94         }
  95     }
  96
  97 }
  98 // ------------------------------------------------------------------------ //
  99 AliPMDCalibGain &AliPMDCalibGain::operator=(const AliPMDCalibGain &pmdcalibgain)
 100 {
 101     if(this != &pmdcalibgain)
 102     {
 103       this->fpw = pmdcalibgain.fpw;
 104         for(Int_t idet = 0; idet < kDet; idet++)
 105         {
 106             for(Int_t ismn = 0; ismn < kMaxSMN; ismn++)
 107             {
 108                 fSMIso[idet][ismn] = pmdcalibgain.fSMIso[idet][ismn];
 109                 fSMCount[idet][ismn] = pmdcalibgain.fSMCount[idet][ismn];
 110                 for(Int_t jrow = 0; jrow < kMaxRow;jrow++)
 111                 {
 112                     for(Int_t kcol = 0; kcol < kMaxCol; kcol++)
 113                       {
 114                         fCellIso[idet][ismn][jrow][kcol]  =
 115                           pmdcalibgain.fCellIso[idet][ismn][jrow][kcol];
 116                         fCellCount[idet][ismn][jrow][kcol]  =
 117                           pmdcalibgain.fCellCount[idet][ismn][jrow][kcol];
 118                         fPedMeanRMS[idet][ismn][jrow][kcol] =
 119                           pmdcalibgain.fPedMeanRMS[idet][ismn][jrow][kcol];
 120
 121                     }
 122                 }
 123             }
 124         }
 125     }
 126     return *this;
 127 }
 128 // ------------------------------------------------------------------------ //
 129 AliPMDCalibGain::~AliPMDCalibGain()
 130 {
 131     // dtor
 132
 133 }
 134
 135 // ------------------------------------------------------------------------ //
 136
 137 Int_t AliPMDCalibGain::ExtractPedestal(const Char_t *rootFile)
 138 {
 139   // Pedestal extraction from the PMD_PED.root file
 140   // To be called once at the beginning
 141
 142   Int_t   det, sm, row, col;
 143   Float_t mean, rms;
 144
 145   TFile *pedfile = new TFile(rootFile);
 146
 147   if(!pedfile)
 148     {
 149       printf("ERROR --- NO PEDESTAL (PMD_PED1.root) FILE IS FOUND --- STOP GAIN DA\n");
 150       return -3;
 151     }
 152
 153
 154   TTree *ped =(TTree*)pedfile->Get("ped");
 155
 156   ped->SetBranchAddress("det",&det);
 157   ped->SetBranchAddress("sm",&sm);
 158   ped->SetBranchAddress("row",&row);
 159   ped->SetBranchAddress("col",&col);
 160   ped->SetBranchAddress("mean",&mean);
 161   ped->SetBranchAddress("rms",&rms);
 162
 163   Int_t nentries = (Int_t)ped->GetEntries();
 164
 165   for (Int_t ient = 0; ient < nentries; ient++)
 166     {
 167       ped->GetEntry(ient);
 168       fPedMeanRMS[det][sm][row][col] = mean + 3.*rms;
 169       //printf("Mean= %f, RMS= %f, PedMeanRMS=%f\n",mean,rms,fPedMeanRMS[det][sm][row][col]);
 170
 171     }
 172
 173   pedfile->Close();
 174   delete pedfile;
 175   pedfile = 0x0;
 176
 177   return 1;
 178 }
 179 // ------------------------------------------------------------------------ //
 180
 181 void AliPMDCalibGain::ReadTempFile(const Char_t *tempFile)
 182 {
 183   // Read the variables from the file
 184
 185   fpw = fopen(tempFile,"r");
 186
 187   Float_t detcount, detiso;
 188   Float_t smcount, smiso;
 189   Float_t cellcount, celliso;
 190
 191   for (Int_t idet = 0; idet < kDet; idet++)
 192     {
 193       fscanf(fpw,"%d %f %f",&idet,&detcount,&detiso);
 194       fDetCount[idet] = detcount;
 195       fDetIso[idet]   = detiso;
 196     }
 197
 198   for (Int_t idet = 0; idet < kDet; idet++)
 199     {
 200       for (Int_t ism = 0; ism < kMaxSMN; ism++)
 201         {
 202           fscanf(fpw,"%d %d %f %f",&idet,&ism,&smcount,&smiso);
 203
 204           fSMCount[idet][ism] = smcount;
 205           fSMIso[idet][ism]   = smiso;
 206         }
 207     }
 208
 209   for (Int_t idet = 0; idet < kDet; idet++)
 210     {
 211       for (Int_t ism = 0; ism < kMaxSMN; ism++)
 212         {
 213           for (Int_t irow = 0; irow < kMaxRow; irow++)
 214             {
 215               for (Int_t icol = 0; icol < kMaxCol; icol++)
 216                 {
 217                   fscanf(fpw,"%d %d %d %d %f %f",&idet,&ism,&irow,&icol,
 218                           &cellcount,&celliso);
 219
 220                   fCellCount[idet][ism][irow][icol] = cellcount;
 221                   fCellIso[idet][ism][irow][icol]   = celliso;
 222                 }
 223             }
 224         }
 225     }
 226
 227   fclose(fpw);
 228
 229 }
 230 // ------------------------------------------------------------------------ //
 231 void AliPMDCalibGain::WriteTempFile(const Char_t *tempFile)
 232 {
 233   // Write the Temporary file if the required statics is not achieved
 234
 235
 236   /*
 237     Following variables to be written to a file
 238     fDetIso[idet] ;
 239     fSMIso[idet][ismn];
 240     fCellIso[idet][ismn][irow][icol];
 241
 242     fDetCount[idet];
 243     fSMCount[idet][ismn];
 244     fCellCount[idet][ismn][irow][icol];
 245   */
 246
 247
 248   fpw = fopen(tempFile,"w+");
 249
 250   for (Int_t idet = 0; idet < kDet; idet++)
 251     {
 252       fprintf(fpw,"%d %f %f\n",idet,fDetCount[idet],fDetIso[idet]);
 253     }
 254
 255   for (Int_t idet = 0; idet < kDet; idet++)
 256     {
 257       for (Int_t ism = 0; ism < kMaxSMN; ism++)
 258         {
 259           fprintf(fpw,"%d %d %f %f\n",idet,ism, fSMCount[idet][ism],fSMIso[idet][ism]);
 260         }
 261     }
 262
 263   for (Int_t idet = 0; idet < kDet; idet++)
 264     {
 265       for (Int_t ism = 0; ism < kMaxSMN; ism++)
 266         {
 267           for (Int_t irow = 0; irow < kMaxRow; irow++)
 268             {
 269               for (Int_t icol = 0; icol < kMaxCol; icol++)
 270                 {
 271                   fprintf(fpw,"%d %d %d %d %f %f\n",idet,ism,irow,icol,
 272                           fCellCount[idet][ism][irow][icol],
 273                           fCellIso[idet][ism][irow][icol]);
 274                 }
 275             }
 276         }
 277     }
 278
 279   fclose(fpw);
 280
 281 }
 282
 283 // ------------------------------------------------------------------------ //
 284
 285 Bool_t AliPMDCalibGain::ProcessEvent(AliRawReader *rawReader, TObjArray *pmdddlcont)
 286 {
 287   // Calculates the ADC of isolated cell
 288
 289   const Int_t kDDL           = AliDAQ::NumberOfDdls("PMD");
 290   const Int_t kCellNeighbour = 6;
 291   Int_t neibx[6] = {1,0,-1,-1,0,1};
 292   Int_t neiby[6] = {0,1,1,0,-1,-1};
 293
 294   Int_t id1,jd1;  //neighbour row/col
 295   Int_t isocount; //number of neighbours with 0 signal
 296
 297   Float_t d1[kDet][kMaxSMN][kMaxRow][kMaxCol];
 298
 299   for(Int_t idet = 0; idet < kDet; idet++)
 300     {
 301       for(Int_t ismn = 0; ismn < kMaxSMN; ismn++)
 302         {
 303           for(Int_t irow = 0; irow < kMaxRow; irow++)
 304             {
 305               for(Int_t icol = 0; icol < kMaxCol; icol++)
 306                 {
 307                   d1[idet][ismn][irow][icol] = 0.;
 308                 }
 309             }
 310         }
 311     }
 312
 313   AliPMDRawStream rawStream(rawReader);
 314
 315   Int_t iddl = -1;
 316
 317   Int_t numberofDDLs = 0;
 318
 319     while ((iddl = rawStream.DdlData(pmdddlcont)) >=0) {
 320       numberofDDLs++;
 321
 322       Int_t ientries = pmdddlcont->GetEntries();
 323
 324       for (Int_t ient = 0; ient < ientries; ient++)
 325       {
 326           AliPMDddldata *pmdddl = (AliPMDddldata*)pmdddlcont->UncheckedAt(ient);
 327
 328           Int_t idet = pmdddl->GetDetector();
 329           Int_t ismn = pmdddl->GetSMN();
 330           Int_t mcm = pmdddl->GetMCM();
 331           //Int_t ichno = pmdddl->GetChannel();
 332           Int_t irow = pmdddl->GetRow();
 333           Int_t icol = pmdddl->GetColumn();
 334           Int_t isig = pmdddl->GetSignal();
 335
 336           // This is the protection not to crash the code
 337
 338           if(mcm == 0) continue;
 339           if (irow < 0 || icol < 0 || irow > 47 || icol > 95) continue;
 340
 341           // Pedestal subtraction
 342
 343           if (isig>0)
 344             {
 345               d1[idet][ismn][irow][icol] =
 346                 (Float_t) isig - fPedMeanRMS[idet][ismn][irow][icol];
 347 //printf("Signal_ped_subtracted=%f, pedestal=%f\n",d1[idet][ismn][irow][icol]),fPedMeanRMS[idet][ismn][irow][icol];
 348             }
 349       }
 350       pmdddlcont->Delete();
 351   }
 352
 353   for(Int_t idet=0; idet < kDet; idet++)
 354   {
 355       for(Int_t ismn = 0; ismn < kMaxSMN; ismn++)
 356       {
 357           for(Int_t irow = 0; irow < kMaxRow; irow++)
 358           {
 359               for(Int_t icol = 0; icol < kMaxCol; icol++)
 360               {
 361                   if(d1[idet][ismn][irow][icol] > 0)
 362                   {
 363                       isocount = 0;
 364                       for(Int_t ii = 0; ii < kCellNeighbour; ii++)
 365                       {
 366                           id1 = irow + neibx[ii];
 367                           jd1 = icol + neiby[ii];
 368                           if(d1[idet][ismn][id1][jd1] == 0)
 369                           {
 370                               isocount++;
 371                               if(isocount == kCellNeighbour)
 372                               {
 373                                   fDetIso[idet] += d1[idet][ismn][irow][icol];
 374                                   fSMIso[idet][ismn] += d1[idet][ismn][irow][icol];
 375                                   fCellIso[idet][ismn][irow][icol] += d1[idet][ismn][irow][icol];
 376                                   fDetCount[idet]++;
 377                                   fSMCount[idet][ismn]++;
 378                                   fCellCount[idet][ismn][irow][icol]++;
 379
 380                               }
 381                           }
 382                       }  // neigh cell cond.
 383                   }     // d>0 cond.
 384               }
 385           }
 386       }
 387   }
 388
 389   if (numberofDDLs < kDDL)
 390       return kFALSE;
 391   return kTRUE;
 392
 393 }
 394 // ------------------------------------------------------------------------ //
 395 void AliPMDCalibGain::Analyse(TTree *gaintree)
 396 {
 397     // Calculates the mean
 398     Int_t   det, sm, row, col;
 399     Float_t gain;
 400     Float_t modmean  = 0.;
 401     Float_t cellmean = 0.;
 402     Float_t detmean =0.;
 403
 404     gaintree->Branch("det",&det,"det/I");
 405     gaintree->Branch("sm",&sm,"sm/I");
 406     gaintree->Branch("row",&row,"row/I");
 407     gaintree->Branch("col",&col,"col/I");
 408     gaintree->Branch("gain",&gain,"gain/F");
 409
 410     for(Int_t idet = 0; idet < kDet; idet++)
 411       {
 412         if (fDetCount[idet]>0 )
 413           detmean=fDetIso[idet]/fDetCount[idet];
 414         for(Int_t ism = 0; ism < kMaxSMN; ism++)
 415           {
 416             if (fSMCount[idet][ism] > 0)
 417               modmean = fSMIso[idet][ism]/fSMCount[idet][ism];
 418             for(Int_t irow = 0; irow < kMaxRow; irow++)
 419               {
 420                 for(Int_t icol = 0; icol < kMaxCol; icol++)
 421                   {
 422                     if (fCellCount[idet][ism][irow][icol] > 0.)
 423                       {
 424                         cellmean = fCellIso[idet][ism][irow][icol]/fCellCount[idet][ism][irow][icol];
 425                       }
 426                     det      = idet;
 427                     sm       = ism;
 428                     row      = irow;
 429                     col      = icol;
 430                     if (cellmean > 0.0 && fCellCount[idet][ism][irow][icol]>0.)
 431                       {
 432                         gain = cellmean/detmean;
 433                       }
 434                     else
 435                       {
 436                         gain = -1.;
 437                       }
 438                     //if(fCellCount[idet][ism][irow][icol]>0.) printf("CellCount =%f, gain= %f\n",fCellCount[idet][ism][irow][icol],gain);
 439                     gaintree->Fill();
 440                   }
 441               }
 442           }
 443       }
 444
 445 }
 446 // ------------------------------------------------------------------------ //