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 "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 "AliPHOSRawDecoderv1.h"
  55 #include "AliPHOSPulseGenerator.h"
  56
  57
  58 ClassImp(AliPHOSRawDecoderv1)
  59
  60 //-----------------------------------------------------------------------------
  61   AliPHOSRawDecoderv1::AliPHOSRawDecoderv1():AliPHOSRawDecoder(),
  62 fSampleParamsLow(0x0),fSampleParamsHigh(0x0),fToFit(0x0)
  63 {
  64   //Default constructor.
  65 }
  66
  67 //-----------------------------------------------------------------------------
  68 AliPHOSRawDecoderv1::AliPHOSRawDecoderv1(AliRawReader* rawReader,  AliAltroMapping **mapping):
  69   AliPHOSRawDecoder(rawReader,mapping),
  70   fSampleParamsLow(0x0),fSampleParamsHigh(0x0),fToFit(0x0)
  71 {
  72   //Construct a decoder object.
  73   //Is is user responsibility to provide next raw event
  74   //using AliRawReader::NextEvent().
  75
  76   if(!gMinuit)
  77     gMinuit = new TMinuit(100);
  78   fSampleParamsHigh =new TArrayD(5) ;
  79   fSampleParamsHigh->AddAt(4.25,0) ;
  80   fSampleParamsHigh->AddAt(0.094,1) ;
  81   fSampleParamsHigh->AddAt(0.0151,2) ;
  82   fSampleParamsHigh->AddAt(0.0384,3) ;
  83   fSampleParamsLow=new TArrayD(5) ;
  84   fSampleParamsLow->AddAt(5.14,0) ;
  85   fSampleParamsLow->AddAt(0.0970,1) ;
  86   fSampleParamsLow->AddAt(0.0088,2) ;
  87   fSampleParamsLow->AddAt(0.0346,3) ;
  88   fToFit = new TList() ;
  89 }
  90
  91 //-----------------------------------------------------------------------------
  92 AliPHOSRawDecoderv1::~AliPHOSRawDecoderv1()
  93 {
  94   //Destructor.
  95   if(fSampleParamsLow){
  96     delete fSampleParamsLow ;
  97     fSampleParamsLow=0 ;
  98   }
  99   if(fSampleParamsHigh){
 100     delete fSampleParamsHigh ;
 101     fSampleParamsHigh=0;
 102   }
 103   if(fToFit){
 104     delete fToFit ;
 105     fToFit=0 ;
 106   }
 107 }
 108
 109 //-----------------------------------------------------------------------------
 110 AliPHOSRawDecoderv1::AliPHOSRawDecoderv1(const AliPHOSRawDecoderv1 &phosDecoder ):
 111   AliPHOSRawDecoder(phosDecoder),
 112   fSampleParamsLow(0x0),fSampleParamsHigh(0x0),fToFit(0x0)
 113 {
 114   //Copy constructor.
 115   fToFit = new TList() ;
 116   fSampleParamsLow =new TArrayD(*(phosDecoder.fSampleParamsLow)) ;
 117   fSampleParamsHigh=new TArrayD(*(phosDecoder.fSampleParamsHigh)) ;
 118 }
 119
 120 //-----------------------------------------------------------------------------
 121 AliPHOSRawDecoderv1& AliPHOSRawDecoderv1::operator = (const AliPHOSRawDecoderv1 &phosDecoder)
 122 {
 123   //Assignment operator.
 124
 125   //  if(this != &phosDecoder) {
 126   //  }
 127   fToFit = new TList() ;
 128   if(fSampleParamsLow){
 129     fSampleParamsLow = phosDecoder.fSampleParamsLow ;
 130     fSampleParamsHigh= phosDecoder.fSampleParamsHigh ;
 131   }
 132   else{
 133     fSampleParamsLow =new TArrayD(*(phosDecoder.fSampleParamsLow)) ;
 134     fSampleParamsHigh=new TArrayD(*(phosDecoder.fSampleParamsHigh)) ;
 135   }
 136   return *this;
 137 }
 138
 139 //-----------------------------------------------------------------------------
 140 Bool_t AliPHOSRawDecoderv1::NextDigit()
 141 {
 142   //Extract an energy deposited in the crystal,
 143   //crystal' position (module,column,row),
 144   //time and gain (high or low).
 145   //First collects sample, then evaluates it and if it has
 146   //reasonable shape, fits it with Gamma2 function and extracts
 147   //energy and time.
 148
 149 //  TCanvas * c  = (TCanvas *)gROOT->FindObjectAny("canvMy") ;
 150 //  TH1S * h = new TH1S("s","",200,0.5,200.5) ;
 151 //  TF1 * fff = new TF1("fff","[0]+[1]*((x-[2])+[3]*(x-[2])*(x-[2]))*(exp(-(x-[2])*[4])+[5]*exp(-(x-[2])*[6]))",0.,1000.) ;
 152
 153   AliCaloRawStream* in = fCaloStream;
 154
 155   Int_t    iBin     = 0;
 156   Int_t    tLength  = 0;
 157   fEnergy = -111;
 158   Float_t pedMean = 0;
 159   Int_t   nPed = 0;
 160   Float_t baseLine = 1.0;
 161   const Float_t nPreSamples = 10;
 162
 163   while ( in->Next() ) {
 164
 165     if(!tLength) {
 166       tLength = in->GetTimeLength();
 167       if(tLength!=fSamples->GetSize()) {
 168         delete fSamples ;
 169         fSamples = new TArrayI(tLength);
 170       }
 171       else{
 172         for(Int_t i=0; i<fSamples->GetSize(); i++){
 173           fSamples->AddAt(0,i) ;
 174         }
 175       }
 176     }
 177
 178     // Fit the full sample
 179     if(in->IsNewHWAddress() && iBin>0) {
 180
 181       Double_t pedestal =0. ;
 182       if(fPedSubtract){
 183         if (nPed > 0)
 184           pedestal = (Double_t)(pedMean/nPed);
 185         else
 186           return kFALSE;
 187       }
 188
 189       //calculate energy
 190       //first estimate if this sample looks like gamma2 function
 191       Double_t aMean=0. ;
 192       Double_t aRMS=0. ;
 193       Int_t tStart = 0 ;
 194       Int_t cnts=0 ;
 195       for(Int_t i=0; i<fSamples->GetSize(); i++){
 196         if(fSamples->At(i)>0){
 197           Double_t de=fSamples->At(i)-pedestal ;
 198           aMean+=de ;
 199           aRMS+=de*de ;
 200           cnts++;
 201           if(de>2 && tStart==0)
 202             tStart=i ;
 203           if(fSamples->At(i)>fEnergy)
 204             fEnergy=fSamples->At(i) ;
 205         }
 206       }
 207       if(cnts>0){
 208         aMean/=cnts;
 209         aRMS=aRMS/cnts-aMean*aMean;
 210       }
 211
 212       //IF sample has reasonable mean and RMS, try to fit it with gamma2
 213       if(fEnergy>2.&& cnts >20 && aMean>0. && aRMS>2.){ //more or less reasonable sample
 214
 215         gMinuit->mncler();                     // Reset Minuit's list of paramters
 216         gMinuit->SetPrintLevel(-1) ;           // No Printout
 217         gMinuit->SetFCN(AliPHOSRawDecoderv1::UnfoldingChiSquare) ;
 218         // To set the address of the minimization function
 219
 220        fToFit->Clear() ;
 221        if(fLowGainFlag){
 222          fSampleParamsLow->AddAt(pedestal,4) ;
 223          fToFit->AddFirst((TObject*)fSampleParamsLow) ;
 224         }
 225         else{
 226          fSampleParamsHigh->AddAt(pedestal,4) ;
 227          fToFit->AddFirst((TObject*)fSampleParamsHigh) ;
 228         }
 229         fToFit->AddLast((TObject*)fSamples) ;
 230
 231         gMinuit->SetObjectFit((TObject*)fToFit) ;         // To tranfer pointer to UnfoldingChiSquare
 232         Int_t ierflg ;
 233         gMinuit->mnparm(0, "t0",  1.*tStart, 0.1, 0, 0, ierflg) ;
 234         if(ierflg != 0){
 235 //        AliWarning(Form("Unable to set initial value for fit procedure : t0=%e\n",1.*tStart) ) ;
 236           fEnergy=0. ;
 237           fTime=-999. ;
 238           return kTRUE ; //will scan further
 239         }
 240         Double_t amp0=(fEnergy-pedestal)*0.0032;
 241
 242         gMinuit->mnparm(1, "Energy", amp0 , 0.001*amp0, 0, 0, ierflg) ;
 243         if(ierflg != 0){
 244 //        AliWarning(Form("Unable to set initial value for fit procedure : E=%e\n", amp0)) ;
 245           fEnergy=0. ;
 246           fTime=-999. ;
 247           return kTRUE ; //will scan further
 248         }
 249
 250         Double_t p0 = 0.0001 ; // "Tolerance" Evaluation stops when EDM = 0.0001*p0 ; The number of function call slightly
 251         //  depends on it.
 252         Double_t p1 = 1.0 ;
 253         Double_t p2 = 0.0 ;
 254         gMinuit->mnexcm("SET STR", &p2, 0, ierflg) ;   // force TMinuit to reduce function calls
 255         gMinuit->mnexcm("SET GRA", &p1, 1, ierflg) ;   // force TMinuit to use my gradient
 256         //      gMinuit->SetMaxIterations(100);
 257         gMinuit->mnexcm("SET NOW", &p2 , 0, ierflg) ;  // No Warnings
 258
 259         gMinuit->mnexcm("MIGRAD", &p0, 0, ierflg) ;    // minimize
 260
 261         Double_t err,t0err ;
 262         Double_t t0,efit ;
 263         gMinuit->GetParameter(0,t0, t0err) ;
 264         gMinuit->GetParameter(1,efit, err) ;
 265
 266         Double_t a,alpha ;
 267         if(fLowGainFlag){
 268           a=fSampleParamsLow->At(0) ;
 269           alpha=fSampleParamsLow->At(1) ;
 270         }
 271         else{
 272           a=fSampleParamsHigh->At(0) ;
 273           alpha=fSampleParamsHigh->At(1) ;
 274         }
 275
 276 //    c->cd() ;
 277 //    h->Draw() ;
 278 //    if(fLowGainFlag){
 279 //      fff->SetParameters(pedestal,efit,t0,a,alpha,fSampleParamsLow->At(2),fSampleParamsLow->At(3)) ;
 280 //    }
 281 //    else{
 282 //      fff->SetParameters(pedestal,efit,t0,a,alpha,fSampleParamsHigh->At(2),fSampleParamsHigh->At(3)) ;
 283 //    }
 284 //    fff->Draw("same") ;
 285 //    c->Update();
 286
 287         efit*=(2.*a+TMath::Sqrt(4.*a*a+alpha*alpha))/alpha/alpha*TMath::Exp(-1.+(alpha-TMath::Sqrt(4.*a*a+alpha*alpha))/2./a) ;
 288 //printf("efit=%f, t0=%e +- %e, ped=%f \n",efit,t0,t0err,pedestal) ;
 289         Double_t fmin,fedm,errdef ;
 290         Int_t npari,nparx,istat;
 291         gMinuit->mnstat(fmin,fedm,errdef,npari,nparx,istat) ;
 292
 293 //if(fLowGainFlag)
 294 // printf("LowGain \n") ;
 295 //else
 296 // printf("highGain \n") ;
 297
 298 //printf("fmin=%e \n",fmin) ;
 299 //getchar() ;
 300
 301         if(1){ //fmin < 3.+0.3*efit ){ //Chi^2 of a good sample
 302           if(efit>0.){
 303             fEnergy=efit ;
 304             fTime=t0 ;
 305           }
 306         }
 307         else{
 308           fEnergy=0 ; //bad sample
 309           fTime=-999.;
 310         }
 311
 312         fTime*=fPulseGenerator->GetRawFormatTimeTrigger() ;
 313
 314         if (fEnergy < baseLine) fEnergy = 0;
 315       }
 316       else{ //bad sample
 317         fEnergy=0. ;
 318         fTime=-999. ;
 319       }
 320
 321       return kTRUE;
 322     }
 323
 324     fLowGainFlag = in->IsLowGain();
 325     fTime = fPulseGenerator->GetRawFormatTimeTrigger() * in->GetTime();
 326     fModule = in->GetModule()+1;
 327     fRow    = in->GetRow()   +1;
 328     fColumn = in->GetColumn()+1;
 329
 330
 331     // Fill array with samples
 332     iBin++;
 333     if(tLength-iBin < nPreSamples) {
 334       pedMean += in->GetSignal();
 335       nPed++;
 336     }
 337     fSamples->AddAt(in->GetSignal(),tLength-iBin);
 338 //    h->SetBinContent(tLength-iBin+1,in->GetSignal()) ;
 339
 340   } // in.Next()
 341
 342   return kFALSE;
 343 }
 344 //_____________________________________________________________________________
 345 void AliPHOSRawDecoderv1::UnfoldingChiSquare(Int_t & /*nPar*/, Double_t * Grad, Double_t & fret, Double_t * x, Int_t iflag)
 346 {
 347   // Calculates the Chi square for the samples minimization
 348   // Number of parameters, Gradient, Chi squared, parameters, what to do
 349
 350   TList * toFit= (TList*)gMinuit->GetObjectFit() ;
 351   TArrayD * params=(TArrayD*)toFit->At(0) ;
 352   TArrayI * samples = (TArrayI*)toFit->At(1) ;
 353
 354   fret = 0. ;
 355   if(iflag == 2)
 356     for(Int_t iparam = 0 ; iparam < 2 ; iparam++)
 357       Grad[iparam] = 0 ; // Will evaluate gradient
 358
 359   Int_t nSamples=samples->GetSize() ; //Math::Min(70,samples->GetSize()) ;
 360   Double_t t0=x[0] ;
 361   Double_t en=x[1] ;
 362   Double_t a=params->At(0) ;
 363   Double_t alpha=params->At(1) ;
 364   Double_t b=params->At(2) ;
 365   Double_t beta=params->At(3) ;
 366
 367   for(Int_t i = 0 ; i < nSamples ; i++) {
 368     if(samples->At(i)==0 || samples->At(i)==1023) //zero or overflow
 369       continue ;
 370     Double_t dt=i*1.-t0 ;
 371     Double_t diff=float(samples->At(i))-Gamma2(dt,en,params) ;
 372     Double_t w=0.1+0.005*i ; //Mean Pedestal RMS + rising modulation
 373     //    if(w==0)w=1. ;
 374     diff/=w ;
 375     if(iflag == 2){  // calculate gradient
 376       if(dt>=0.){
 377         Grad[0] += -2.*en*diff*((alpha*dt*(1.+a*dt)-1.-2.*a*dt)*TMath::Exp(-alpha*dt)+
 378                               b*(beta*dt*(1.+a*dt)-1.-2.*a*dt)*TMath::Exp(-beta*dt)) /w ; //derivative over t0
 379         Grad[1] += -2.*diff*(dt+a*dt*dt)*(TMath::Exp(-alpha*dt)+
 380                      b*TMath::Exp(-dt*beta))/w ;
 381       }
 382     }
 383     fret += diff*diff ;
 384   }
 385   if(nSamples){
 386     fret/=nSamples ;
 387     if(iflag == 2){
 388       for(Int_t iparam = 0 ; iparam < 2 ; iparam++)
 389         Grad[iparam] /= nSamples ;
 390     }
 391   }
 392
 393 }
 394 //-----------------------------------------------------------------------------
 395 Double_t AliPHOSRawDecoderv1::Gamma2(Double_t dt,Double_t en,TArrayD * params){  //Function for fitting samples
 396   //parameters:
 397   //dt-time after start
 398   //en-amplutude
 399   //function parameters
 400
 401   Double_t ped=params->At(4) ;
 402   if(dt<0.)
 403     return ped ; //pedestal
 404   else
 405     return ped+en*(dt+params->At(0)*dt*dt)*(TMath::Exp(-dt*params->At(1))+params->At(2)*TMath::Exp(-dt*params->At(3))) ;
 406 }
 407