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   Int_t   nPed = 0;
 149   Double_t baseLine = 1.0;
 150   const Int_t nPreSamples = 10;
 151
 152   while ( in->Next() ) {
 153
 154     // Fit the full sample
 155     if(in->IsNewHWAddress() && iBin!=fSamples->GetSize()) {
 156
 157       Double_t pedestal =0. ;
 158       if(fPedSubtract){
 159         if (nPed > 0)
 160           pedestal = (Double_t)(pedMean/nPed);
 161         else
 162           return kFALSE;
 163       }
 164 //      for(Int_t i=0; i<fSamples->GetSize(); i++){
 165 //        h->SetBinContent(i+1,fSamples->At(i)) ;
 166 //      }
 167
 168       //calculate time and energy
 169       Int_t maxBin=0 ;
 170       Int_t maxAmp=0 ;
 171       for(Int_t i=iBin; i<fSamples->GetSize(); i++){
 172         Double_t de=fSamples->At(i)-pedestal ;
 173         if(fEnergy<de){
 174           fEnergy=de ;
 175           maxBin=i ;
 176           maxAmp=fSamples->At(i) ;
 177         }
 178       }
 179       fOverflow =0 ;  //look for plato on the top of sample
 180       if(fEnergy>700 &&  //this is not fluctuation of soft sample
 181          maxBin<fSamples->GetSize()-1 && fSamples->At(maxBin+1)==maxAmp){ //and there is a plato
 182          fOverflow = kTRUE ;
 183       }
 184
 185 //    if(fEnergy>1020.){
 186 //    printf("fE=%f \n",fEnergy) ;
 187 //    cs->cd() ;
 188 //    h->Draw() ;
 189 //    cs->Update() ;
 190 //    getchar() ;
 191 //    }
 192
 193       if(fOverflow)
 194         return kTRUE ; //do not calculate energy and time for overflowed channels
 195
 196       if(fEnergy<baseLine){ //do not evaluate time, drop this sample
 197         fEnergy=0. ;
 198         fTime=-999.;
 199         return kTRUE ;
 200       }
 201
 202       //else calculate time
 203       fTime=0. ;
 204       Double_t tRMS = 0. ;
 205       Double_t tW = 0. ;
 206       Int_t cnts=0 ;
 207       Double_t a,b,c ;
 208       if(fLowGainFlag){
 209         a=fLGpar[0] ;
 210         b=fLGpar[1] ;
 211         c=fLGpar[2] ;
 212       }
 213       else{
 214         a=fHGpar[0] ;
 215         b=fHGpar[1] ;
 216         c=fHGpar[2] ;
 217       }
 218
 219       for(Int_t i=iBin+1; i<fSamples->GetSize()&& cnts<fNtimeSamples; i++){
 220         if(fSamples->At(i)<pedestal)
 221           continue ;
 222 //Presently we do not use channels with overflow
 223 //        if(fOverflow && (fSamples->At(i)==maxAmp ||fSamples->At(i-1)==maxAmp)) //can not calculate time in overflow bin
 224 //          continue ;
 225         if(fTimes->At(i)-fTimes->At(i-1)!=1) //do not use samples with non-consequtive points
 226           continue ;
 227         Double_t de=fSamples->At(i)-pedestal ;
 228         Double_t av = de+fSamples->At(i-1)-pedestal ;
 229         if(av<=0.) //this is fluctuation around pedestal, scip
 230           continue ;
 231         Double_t ds = fSamples->At(i)-fSamples->At(i-1) ;
 232         Double_t ti = ds/av ;     // calculate log. derivative
 233         ti=a/(ti+b)-c*ti ;        // and compare with parameterization
 234         ti=Double_t(fTimes->At(i))-ti ;
 235         Double_t wi = TMath::Abs(ds) ;
 236         fTime+=ti*wi ;
 237         tW+=wi;
 238         tRMS+=ti*ti*wi ;
 239         cnts++ ;
 240       }
 241
 242       if(tW>0.){
 243         fTime/=tW ;
 244         tRMS/=tW ;
 245         tRMS-=fTime*fTime ;
 246 //printf("t0=%f, RMS=%f, cut=%f \n",fTime,tRMS,fRMScut) ;
 247         if(tRMS>=fRMScut){ //bad sample
 248           fTime=-999. ;
 249           fEnergy=0. ;
 250         }
 251       }
 252       else{
 253         fTime=-999. ;
 254       }
 255
 256       fTime*=fPulseGenerator->GetRawFormatTimeTrigger() ;
 257
 258       return kTRUE ; //will scan further
 259     }
 260
 261     fLowGainFlag = in->IsLowGain();
 262     fModule = in->GetModule()+1;
 263     fRow    = in->GetRow()   +1;
 264     fColumn = in->GetColumn()+1;
 265
 266
 267     // Fill array with samples
 268     iBin--;
 269     if(fPedSubtract && (in->GetTime() < nPreSamples)) {
 270       pedMean += in->GetSignal();
 271       nPed++;
 272     }
 273     fSamples->AddAt(in->GetSignal(),iBin);
 274     fTimes->AddAt(in->GetTime(),iBin);
 275   } // in.Next()
 276
 277   return kFALSE;
 278 }