PHOS/AliPHOSRawDecoderv2.cxx

   1 /**************************************************************************
   2  * Copyright(c) 2007, ALICE Experiment at CERN, All rights reserved.      *
   3  *                                                                        *
   4  * Author: The ALICE Off-line Project.                                    *
   5  * Contributors are mentioned in the code where appropriate.              *
   6  *                                                                        *
   7  * Permission to use, copy, modify and distribute this software and its   *
   8  * documentation strictly for non-commercial purposes is hereby granted   *
   9  * without fee, provided that the above copyright notice appears in all   *
  10  * copies and that both the copyright notice and this permission notice   *
  11  * appear in the supporting documentation. The authors make no claims     *
  12  * about the suitability of this software for any purpose. It is          *
  13  * provided "as is" without express or implied warranty.                  *
  14  **************************************************************************/
  15
  16 /* $Id$ */
  17
  18 // This class decodes the stream of ALTRO samples to extract
  19 // the PHOS "digits" of current event. Uses fitting procedure
  20 // to separate reasonable samples
  21 //
  22 // Typical use case:
  23 //     AliRawReader* rf = new AliRawReaderDate("2006run2211.raw");
  24 //     AliPHOSRawDecoder dc(rf);
  25 //     while (rf->NextEvent()) {
  26 //       dc.SetOldRCUFormat(kTRUE);
  27 //       dc.SubtractPedestals(kTRUE);
  28 //       while ( dc.NextDigit() ) {
  29 //         Int_t module = dc.GetModule();
  30 //         Int_t column = dc.GetColumn();
  31 //         Int_t row = dc.GetRow();
  32 //         Double_t amplitude = dc.GetEnergy();
  33 //         Double_t time = dc.GetTime();
  34 //         Bool_t IsLowGain = dc.IsLowGain();
  35 //            ..........
  36 //       }
  37 //     }
  38
  39 // Author: Dmitri Peressounko using idea of Y.Kucheryaev
  40
  41 // --- ROOT system ---
  42 #include "TList.h"
  43 #include "TMath.h"
  44 #include "TMinuit.h"
  45
  46 #include "TCanvas.h"
  47 #include "TH1.h"
  48 #include "TH2.h"
  49 #include "TF1.h"
  50 #include "TROOT.h"
  51
  52 // --- AliRoot header files ---
  53 //#include "AliLog.h"
  54 #include "AliPHOSRawDecoderv2.h"
  55 #include "AliPHOSPulseGenerator.h"
  56
  57
  58 ClassImp(AliPHOSRawDecoderv2)
  59
  60 //-----------------------------------------------------------------------------
  61   AliPHOSRawDecoderv2::AliPHOSRawDecoderv2():AliPHOSRawDecoder(),
  62 fNtimeSamples(25),fRMScut(11.)
  63 {
  64   //Default constructor.
  65   fLGpar[0]=0.971 ;
  66   fLGpar[1]=0.0465;
  67   fLGpar[2]=1.56  ;
  68   fHGpar[0]=0.941 ;
  69   fHGpar[1]=0.0436;
  70   fHGpar[2]=1.65  ;
  71 }
  72
  73 //-----------------------------------------------------------------------------
  74 AliPHOSRawDecoderv2::AliPHOSRawDecoderv2(AliRawReader* rawReader,  AliAltroMapping **mapping):
  75   AliPHOSRawDecoder(rawReader,mapping),
  76 fNtimeSamples(25),fRMScut(11.)
  77 {
  78   //Construct a decoder object.
  79   //Is is user responsibility to provide next raw event
  80   //using AliRawReader::NextEvent().
  81   fLGpar[0]=0.971 ;
  82   fLGpar[1]=0.0465;
  83   fLGpar[2]=1.56  ;
  84   fHGpar[0]=0.941 ;
  85   fHGpar[1]=0.0436;
  86   fHGpar[2]=1.65  ;
  87 }
  88
  89 //-----------------------------------------------------------------------------
  90 AliPHOSRawDecoderv2::~AliPHOSRawDecoderv2()
  91 {
  92   //Destructor.
  93   //Nothing to delete
  94 }
  95
  96 //-----------------------------------------------------------------------------
  97 AliPHOSRawDecoderv2::AliPHOSRawDecoderv2(const AliPHOSRawDecoderv2 &phosDecoder ):
  98   AliPHOSRawDecoder(phosDecoder),
  99 fNtimeSamples(25),fRMScut(11.)
 100 {
 101   //Copy constructor.
 102   fNtimeSamples=phosDecoder.fNtimeSamples ;
 103   for(Int_t i=0; i<3;i++){
 104     fLGpar[i]=phosDecoder.fLGpar[i] ;
 105     fHGpar[i]=phosDecoder.fHGpar[i] ;
 106   }
 107   fRMScut=phosDecoder.fRMScut ;
 108 }
 109
 110 //-----------------------------------------------------------------------------
 111 AliPHOSRawDecoderv2& AliPHOSRawDecoderv2::operator = (const AliPHOSRawDecoderv2 &phosDecoder)
 112 {
 113   //Assignment operator.
 114
 115   fNtimeSamples=phosDecoder.fNtimeSamples ;
 116   for(Int_t i=0; i<3;i++){
 117     fLGpar[i]=phosDecoder.fLGpar[i] ;
 118     fHGpar[i]=phosDecoder.fHGpar[i] ;
 119   }
 120   fRMScut=phosDecoder.fRMScut ;
 121   return *this;
 122 }
 123
 124 //-----------------------------------------------------------------------------
 125 Bool_t AliPHOSRawDecoderv2::NextDigit()
 126 {
 127   //Extract an energy deposited in the crystal,
 128   //crystal' position (module,column,row),
 129   //time and gain (high or low).
 130   //First collects sample, then evaluates it and if it has
 131   //reasonable shape, fits it with Gamma2 function and extracts
 132   //energy and time.
 133
 134   TCanvas * cs = (TCanvas*)gROOT->FindObjectAny("CSample") ;
 135   if(!cs)
 136     cs = new TCanvas("CSample","CSample") ;
 137
 138   TH1D * h = (TH1D*)gROOT->FindObjectAny("hSample") ;
 139   if(!h)
 140     h=new TH1D("hSample","",200,0.,200.) ;
 141
 142
 143   AliCaloRawStream* in = fCaloStream;
 144
 145   Int_t    iBin     = fSamples->GetSize();
 146   fEnergy = 0;
 147   Double_t pedMean = 0;
 148   Double_t pedRMS = 0;
 149   Int_t   nPed = 0;
 150   Double_t baseLine = 1.0;
 151   const Int_t nPreSamples = 10;
 152   fQuality = 0. ;
 153
 154   while ( in->Next() ) {
 155
 156     // Fit the full sample
 157     if(in->IsNewHWAddress() && iBin!=fSamples->GetSize()) {
 158
 159       Double_t pedestal =0. ;
 160       if(fPedSubtract){
 161         if (nPed > 0)
 162           pedestal = (Double_t)(pedMean/nPed);
 163         else
 164           return kFALSE;
 165       }
 166       for(Int_t i=0; i<fSamples->GetSize(); i++){
 167         h->SetBinContent(i+1,fSamples->At(i)) ;
 168       }
 169
 170       //calculate time and energy
 171       Int_t maxBin=0 ;
 172       Int_t maxAmp=0 ;
 173       for(Int_t i=iBin; i<fSamples->GetSize(); i++){
 174         if(maxAmp<fSamples->At(i)){
 175           maxBin=i ;
 176           maxAmp=fSamples->At(i) ;
 177         }
 178       }
 179       if(maxBin==fSamples->GetSize()-1){//bad sample
 180         fEnergy=0. ;
 181         fTime=-999.;
 182         return kTRUE ;
 183       }
 184       fEnergy=Double_t(maxAmp)-pedestal ;
 185       fOverflow =0 ;  //look for plato on the top of sample
 186       if(fEnergy>500 &&  //this is not fluctuation of soft sample
 187          maxBin<fSamples->GetSize()-1 && fSamples->At(maxBin+1)==maxAmp){ //and there is a plato
 188          fOverflow = kTRUE ;
 189       }
 190
 191 //    if(fEnergy>500.){
 192 // if(fRow==54 && fColumn==24){
 193 //    printf("fE=%f, ped=%f, row=%d, col=%d \n",fEnergy,pedestal,fRow,fColumn) ;
 194 //    if(fOverflow)printf(" Overflow \n") ;
 195 //    else printf("iBin=%d, maxBin=%d, maxAmp=%d,Amp(+1)=%d,Amp(-1)=%d  \n",iBin,maxBin,maxAmp,fSamples->At(maxBin+1),fSamples->At(maxBin-1)) ;
 196 //    cs->cd() ;
 197 //    h->Draw() ;
 198 //    cs->Update() ;
 199 //    getchar() ;
 200 // }
 201
 202       if(fOverflow)
 203         return kTRUE ; //do not calculate energy and time for overflowed channels
 204
 205       if(fEnergy<baseLine){ //do not evaluate time, drop this sample
 206         fEnergy=0. ;
 207         fTime=-999.;
 208         return kTRUE ;
 209       }
 210
 211       //else calculate time
 212       fTime=0. ;
 213       Double_t tRMS = 0. ;
 214       Double_t tW = 0. ;
 215       Int_t cnts=0 ;
 216       Double_t a,b,c ;
 217       if(fLowGainFlag){
 218         a=fLGpar[0] ;
 219         b=fLGpar[1] ;
 220         c=fLGpar[2] ;
 221       }
 222       else{
 223         a=fHGpar[0] ;
 224         b=fHGpar[1] ;
 225         c=fHGpar[2] ;
 226       }
 227
 228       for(Int_t i=iBin+1; i<fSamples->GetSize()&& cnts<fNtimeSamples; i++){
 229         if(fSamples->At(i)<pedestal)
 230           continue ;
 231 //Presently we do not use channels with overflow
 232 //        if(fOverflow && (fSamples->At(i)==maxAmp ||fSamples->At(i-1)==maxAmp)) //can not calculate time in overflow bin
 233 //          continue ;
 234         if(fTimes->At(i)-fTimes->At(i-1)!=1) //do not use samples with non-consequtive points
 235           continue ;
 236         Double_t de=fSamples->At(i)-pedestal ;
 237         Double_t av = de+fSamples->At(i-1)-pedestal ;
 238         if(av<=0.) //this is fluctuation around pedestal, scip
 239           continue ;
 240         Double_t ds = fSamples->At(i)-fSamples->At(i-1) ;
 241         Double_t ti = ds/av ;     // calculate log. derivative
 242         ti=a/(ti+b)-c*ti ;        // and compare with parameterization
 243         ti=Double_t(fTimes->At(i))-ti ;
 244         Double_t wi = TMath::Abs(ds) ;
 245         fTime+=ti*wi ;
 246         tW+=wi;
 247         tRMS+=ti*ti*wi ;
 248         cnts++ ;
 249       }
 250
 251       if(tW>0.){
 252         fTime/=tW ;
 253         fQuality = tRMS/tW-fTime*fTime ;
 254         //Normalize quality
 255 //printf("t0=%f, RMS=%f, cut=%f \n",fTime,tRMS,fRMScut) ;
 256 //        if(tRMS>=fRMScut){ //bad sample
 257 //          fTime=-999. ;
 258 //          fEnergy=0. ;
 259 //        }
 260       }
 261       else{
 262         fTime=-999. ;
 263         fQuality=999. ;
 264       }
 265
 266       Bool_t isBad = 0 ;
 267       for(Int_t i=iBin+1; i<fSamples->GetSize()-1&&!isBad; i++){
 268         if(fSamples->At(i)>fSamples->At(i-1)+5 && fSamples->At(i)>fSamples->At(i+1)+5) { //single jump
 269           isBad=1 ;
 270         }
 271       }
 272       if(pedestal<10.)
 273         isBad=1 ;
 274
 275       pedRMS=pedRMS/nPed-pedestal*pedestal ;
 276       if(pedRMS>0.1)
 277         isBad=1 ;
 278
 279       for(Int_t i=iBin+1; i<fSamples->GetSize()-1&&!isBad; i++){
 280          if(fSamples->At(i)<pedestal-1)
 281            isBad=1 ;
 282       }
 283
 284       //two maxima
 285
 286
 287     if(fEnergy>10. && !isBad ){
 288     printf("fE=%f, ped=%f, fQuality=%f, pedRMS=%f \n",fEnergy,pedestal,fQuality,pedRMS) ;
 289     if(fOverflow)printf(" Overflow \n") ;
 290     if(isBad)printf("bad") ;
 291 //    else printf("iBin=%d, maxBin=%d, maxAmp=%d,Amp(+1)=%d,Amp(-1)=%d  \n",iBin,maxBin,maxAmp,fSamples->At(maxBin+1),fSamples->At(maxBin-1)) ;
 292     cs->cd() ;
 293     h->Draw() ;
 294     cs->Update() ;
 295     getchar() ;
 296  }
 297
 298
 299       return kTRUE ; //will scan further
 300     }
 301
 302     fLowGainFlag = in->IsLowGain();
 303     fModule = in->GetModule()+1;
 304     fRow    = in->GetRow()   +1;
 305     fColumn = in->GetColumn()+1;
 306
 307
 308     // Fill array with samples
 309     iBin--;
 310     if(fPedSubtract && (in->GetTime() < nPreSamples)) {
 311       pedMean += in->GetSignal();
 312       pedRMS += in->GetSignal()*in->GetSignal();
 313       nPed++;
 314     }
 315     fSamples->AddAt(in->GetSignal(),iBin);
 316     fTimes->AddAt(in->GetTime(),iBin);
 317   } // in.Next()
 318
 319   return kFALSE;
 320 }