PHOS/AliPHOSRawDecoderv1.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
  40
  41 // --- ROOT system ---
  42 #include "TMath.h"
  43 #include "TMinuit.h"
  44
  45 // --- AliRoot header files ---
  46 #include "AliPHOSRawDecoderv1.h"
  47 #include "AliPHOSPulseGenerator.h"
  48
  49
  50 ClassImp(AliPHOSRawDecoderv1)
  51
  52 //-----------------------------------------------------------------------------
  53   AliPHOSRawDecoderv1::AliPHOSRawDecoderv1():AliPHOSRawDecoder()
  54 {
  55   //Default constructor.
  56 }
  57
  58 //-----------------------------------------------------------------------------
  59 AliPHOSRawDecoderv1::AliPHOSRawDecoderv1(AliRawReader* rawReader,  AliAltroMapping **mapping):
  60   AliPHOSRawDecoder(rawReader,mapping)
  61 {
  62   //Construct a decoder object.
  63   //Is is user responsibility to provide next raw event
  64   //using AliRawReader::NextEvent().
  65
  66   if(!gMinuit)
  67     gMinuit = new TMinuit(100);
  68
  69 }
  70
  71 //-----------------------------------------------------------------------------
  72 AliPHOSRawDecoderv1::~AliPHOSRawDecoderv1()
  73 {
  74   //Destructor.
  75 }
  76
  77 //-----------------------------------------------------------------------------
  78 AliPHOSRawDecoderv1::AliPHOSRawDecoderv1(const AliPHOSRawDecoderv1 &phosDecoder ):
  79   AliPHOSRawDecoder(phosDecoder)
  80 {
  81   //Copy constructor.
  82 }
  83
  84 //-----------------------------------------------------------------------------
  85 AliPHOSRawDecoderv1& AliPHOSRawDecoderv1::operator = (const AliPHOSRawDecoderv1 &phosDecode)
  86 {
  87   //Assignment operator.
  88
  89   //  if(this != &phosDecode) {
  90   //  }
  91
  92   return *this;
  93 }
  94
  95 //-----------------------------------------------------------------------------
  96 Bool_t AliPHOSRawDecoderv1::NextDigit()
  97 {
  98   //Extract an energy deposited in the crystal,
  99   //crystal' position (module,column,row),
 100   //time and gain (high or low).
 101   //First collects sample, then evaluates it and if it has
 102   //reasonable shape, fits it with Gamma2 function and extracts
 103   //energy and time.
 104
 105   AliCaloRawStream* in = fCaloStream;
 106
 107   Int_t    iBin     = 0;
 108   Int_t    tLength  = 0;
 109   fEnergy = -111;
 110   Float_t pedMean = 0;
 111   Int_t   nPed = 0;
 112   Float_t baseLine = 1.0;
 113   const Float_t nPreSamples = 10;
 114
 115   while ( in->Next() ) {
 116
 117     if(!tLength) {
 118       tLength = in->GetTimeLength();
 119       if(tLength!=fSamples->GetSize()) {
 120         delete fSamples ;
 121         fSamples = new TArrayI(tLength);
 122       }
 123       else{
 124         for(Int_t i=0; i<fSamples->GetSize(); i++){
 125           fSamples->AddAt(0,i) ;
 126         }
 127       }
 128     }
 129
 130     // Fit the full sample
 131     if(in->IsNewHWAddress() && iBin>0) {
 132
 133       Double_t pedestal =0. ;
 134       if(fPedSubtract){
 135         if (nPed > 0)
 136           pedestal = (Double_t)(pedMean/nPed);
 137         else
 138           return kFALSE;
 139       }
 140
 141       //calculate energy
 142       //first estimate if this sample looks like gamma2 function
 143       Double_t aMean=0. ;
 144       Double_t aRMS=0. ;
 145       Int_t tStart = 0 ;
 146       Int_t cnts=0 ;
 147       for(Int_t i=0; i<fSamples->GetSize(); i++){
 148         if(fSamples->At(i)>0){
 149           Double_t de=fSamples->At(i)-pedestal ;
 150           aMean+=de ;
 151           aRMS+=de*de ;
 152           cnts++;
 153           if(de>2 && tStart==0)
 154             tStart=i ;
 155           if(fSamples->At(i)>fEnergy)
 156             fEnergy=fSamples->At(i) ;
 157         }
 158       }
 159       if(cnts>0){
 160         aMean/=cnts;
 161         aRMS=aRMS/cnts-aMean*aMean;
 162       }
 163
 164       if(fPedSubtract){
 165         fEnergy-=pedestal ;
 166       }
 167
 168       //IF sample has reasonable mean and RMS, try to fit it with gamma2
 169       if(fEnergy>2.&& cnts >20 && aMean>0. && aRMS>2.){ //more or less reasonable sample
 170
 171         gMinuit->mncler();                     // Reset Minuit's list of paramters
 172         gMinuit->SetPrintLevel(-1) ;           // No Printout
 173         gMinuit->SetFCN(AliPHOSRawDecoderv1::UnfoldingChiSquare) ;
 174         // To set the address of the minimization function
 175
 176         gMinuit->SetObjectFit((TObject*)fSamples) ;         // To tranfer pointer to UnfoldingChiSquare
 177         Int_t ierflg ;
 178         gMinuit->mnparm(0, "p",  pedestal, 0.1, 0, 0, ierflg) ;
 179         if(ierflg != 0){
 180           //      Warning("NextDigit", "Uunable to set initial value for fit procedure : ped = %f\n", pedestal ) ;
 181           fEnergy=0. ;
 182           fTime=-999. ;
 183           return kTRUE ; //will scan further
 184         }
 185         gMinuit->mnparm(1, "t0",  1.*tStart, 0.1, 0, 0, ierflg) ;
 186         if(ierflg != 0){
 187           //      Warning("NextDigit", "Unable to set initial value for fit procedure : t0\n" ) ;
 188           fEnergy=0. ;
 189           fTime=-999. ;
 190           return kTRUE ; //will scan further
 191         }
 192         gMinuit->mnparm(2, "Energy", fEnergy*0.018 , 0.001*fEnergy, 0, 0, ierflg) ;
 193         if(ierflg != 0){
 194           //      Warning("NextDigit", "Unable to set initial value for fit procedure : E=%f\n", fEnergy*0.018) ;
 195           fEnergy=0. ;
 196           fTime=-999. ;
 197           return kTRUE ; //will scan further
 198         }
 199         gMinuit->mnparm(3, "Slope", 0.09 , 0.001, 0.001, 1., ierflg) ;
 200         if(ierflg != 0){
 201           //      Warning("NextDigit", "Unable to set initial value for fit procedure : Slope\n") ;
 202           fEnergy=0. ;
 203           fTime=-999. ;
 204           return kTRUE ; //will scan further
 205         }
 206
 207         Double_t p0 = 1. ; // "Tolerance" Evaluation stops when EDM = 0.0001*p0 ; The number of function call slightly
 208         //  depends on it.
 209         Double_t p1 = 1.0 ;
 210         Double_t p2 = 0.0 ;
 211         gMinuit->mnexcm("SET STR", &p2, 0, ierflg) ;   // force TMinuit to reduce function calls
 212         gMinuit->mnexcm("SET GRA", &p1, 1, ierflg) ;   // force TMinuit to use my gradient
 213         //      gMinuit->SetMaxIterations(100);
 214         gMinuit->mnexcm("SET NOW", &p2 , 0, ierflg) ;  // No Warnings
 215
 216         gMinuit->mnexcm("MIGRAD", &p0, 0, ierflg) ;    // minimize
 217
 218         Double_t err ;
 219         Double_t p,t0,efit,slope ;
 220         gMinuit->GetParameter(0,p, err) ;
 221         gMinuit->GetParameter(1,t0, err) ;
 222         gMinuit->GetParameter(2,efit, err) ;
 223         gMinuit->GetParameter(3,slope, err) ;
 224
 225         efit=efit*4.*TMath::Exp(-2.)/slope/slope ; //slope can not be zero - par limits
 226
 227         Double_t fmin,fedm,errdef ;
 228         Int_t npari,nparx,istat;
 229         gMinuit->mnstat(fmin,fedm,errdef,npari,nparx,istat) ;
 230
 231         if(fmin < cnts*(efit/10.) ){ //Chi^2 of a good sample
 232           if(efit<1050.&& efit>0.){
 233             fEnergy=efit ;
 234             fTime=t0 ;
 235           }
 236         }
 237         else{
 238           fEnergy=0 ; //bad sample
 239           fTime=-999.;
 240         }
 241         if(fLowGainFlag)
 242           fEnergy *= fPulseGenerator->GetRawFormatHighLowGainFactor(); // *16
 243
 244         if (fEnergy < baseLine) fEnergy = 0;
 245       }
 246       else{ //bad sample
 247         fEnergy=0. ;
 248         fTime=-999. ;
 249       }
 250
 251       return kTRUE;
 252     }
 253
 254     fLowGainFlag = in->IsLowGain();
 255     fTime = fPulseGenerator->GetRawFormatTimeTrigger() * in->GetTime();
 256     fModule = in->GetModule()+1;
 257     fRow    = in->GetRow()   +1;
 258     fColumn = in->GetColumn()+1;
 259
 260
 261     // Fill array with samples
 262     iBin++;
 263     if(tLength-iBin < nPreSamples) {
 264       pedMean += in->GetSignal();
 265       nPed++;
 266     }
 267     fSamples->AddAt(in->GetSignal(),tLength-iBin);
 268
 269   } // in.Next()
 270
 271   return kFALSE;
 272 }
 273 //_____________________________________________________________________________
 274 void AliPHOSRawDecoderv1::UnfoldingChiSquare(Int_t & /*nPar*/, Double_t * Grad, Double_t & fret, Double_t * x, Int_t iflag)
 275 {
 276   // Calculates the Chi square for the samples minimization
 277   // Number of parameters, Gradient, Chi squared, parameters, what to do
 278
 279   TArrayI * samples = (TArrayI*)gMinuit->GetObjectFit() ;
 280
 281   fret = 0. ;
 282   if(iflag == 2)
 283     for(Int_t iparam = 0 ; iparam < 4 ; iparam++)
 284       Grad[iparam] = 0 ; // Will evaluate gradient
 285
 286   Int_t nSamples=samples->GetSize() ;
 287   Double_t p =x[0] ;
 288   Double_t t0=x[1] ;
 289   Double_t en=x[2] ;
 290   Double_t a =x[3] ;
 291
 292   for(Int_t i = 0 ; i < nSamples ; i++) {
 293     if(samples->At(i)==0)
 294       continue ;
 295     Double_t dt=i*1.-t0 ;
 296     Double_t diff=samples->At(i)-Gamma2(dt,p,en,a) ;
 297     Double_t w=1. ; //TMath::Ceil(TMath::Abs(samples->At(i)-ped)/10.) ;
 298     //    if(w==0)w=1. ;
 299     diff/=w ;
 300     if(iflag == 2){  // calculate gradient
 301       Grad[0] += -2.*diff/w ; //derivative over pedestal
 302       if(dt>=0.){
 303         Grad[1] += -2.*en*diff*(a*dt-2.*dt)*dt*TMath::Exp(-a*dt)/w ; //derivative over t0
 304         Grad[2] += -2.*diff*dt*dt*TMath::Exp(-a*dt)/w ;
 305         Grad[3] +=  2.*en*diff*dt*dt*dt*TMath::Exp(-a*dt)/w ;
 306       }
 307     }
 308     fret += diff*diff ;
 309   }
 310   /*
 311   if(nSamples){
 312     fret/=nSamples ;
 313     if(iflag == 2){
 314       for(Int_t iparam = 0 ; iparam < 4 ; iparam++)
 315         Grad[iparam] /= nSamples ;
 316     }
 317   }
 318   */
 319
 320
 321 }
 322 //-----------------------------------------------------------------------------
 323 Double_t AliPHOSRawDecoderv1::Gamma2(Double_t dt,Double_t p,Double_t en,Double_t a){
 324   //Function for fitting samples
 325   //parameters:
 326   //p-pedestal
 327   //en-amplutude
 328   //a-decay time
 329
 330   if(dt<0.)
 331     return p ; //pedestal
 332   else
 333     return p+en*dt*dt*TMath::Exp(-dt*a) ;
 334 }
 335