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