EMCAL/AliCaloRawAnalyzerLMS.cxx

   1 /**************************************************************************
   2  * This file is property of and copyright by                              *
   3  * the Relativistic Heavy Ion Group (RHIG), Yale University, US, 2009     *
   4  *                                                                        *
   5  * Primary Author: Per Thomas Hille <p.t.hille@fys.uio.no>                *
   6  *                                                                        *
   7  * Contributors are mentioned in the code where appropriate.              *
   8  * Please report bugs to p.t.hille@fys.uio.no                             *
   9  *                                                                        *
  10  * Permission to use, copy, modify and distribute this software and its   *
  11  * documentation strictly for non-commercial purposes is hereby granted   *
  12  * without fee, provided that the above copyright notice appears in all   *
  13  * copies and that both the copyright notice and this permission notice   *
  14  * appear in the supporting documentation. The authors make no claims     *
  15  * about the suitability of this software for any purpose. It is          *
  16  * provided "as is" without express or implied warranty.                  *
  17  **************************************************************************/
  18
  19
  20 // Extraction of amplitude and peak position
  21 // FRom CALO raw data using
  22 // least square fit for the
  23 // Moment assuming identical and
  24 // independent errors (equivalent with chi square)
  25 //
  26
  27 #include "AliCaloRawAnalyzerLMS.h"
  28 #include "AliCaloBunchInfo.h"
  29 #include "AliCaloFitResults.h"
  30 #include "AliLog.h"
  31 #include "TMath.h"
  32 #include <iostream>
  33 #include "TF1.h"
  34 #include "TGraph.h"
  35
  36 using namespace std;
  37
  38
  39 #define BAD 4  //CRAP PTH
  40
  41 ClassImp( AliCaloRawAnalyzerLMS )
  42
  43
  44 AliCaloRawAnalyzerLMS::AliCaloRawAnalyzerLMS() : AliCaloRawAnalyzer("Chi Square Fit"),
  45                                                  fkEulerSquared(7.389056098930650227),
  46                                                  fSig(0),
  47                                                  fTf1(0)
  48 {
  49   //comment
  50   for(int i=0; i < MAXSAMPLES; i++)
  51     {
  52       fXaxis[i] = i;
  53     }
  54
  55   fTf1 = new TF1( "myformula", "[0]*((x - [1])/2.35)^2*exp(-2*(x -[1])/2.35)",  0, 30 );
  56   //  fTf1 = new TF1( "myformula", "[0]*((x - [1])/[2])^[3]*exp(-[3]*(x -[1])/[2])",  0, 30 );
  57
  58 }
  59
  60
  61 AliCaloRawAnalyzerLMS::~AliCaloRawAnalyzerLMS()
  62 {
  63   delete fTf1;
  64 }
  65
  66
  67
  68 AliCaloFitResults
  69 AliCaloRawAnalyzerLMS::Evaluate( const vector<AliCaloBunchInfo>  &bunchvector, const UInt_t altrocfg1,  const UInt_t altrocfg2 )
  70 {
  71   // Extracting signal parameters using fitting
  72   short maxampindex; //index of maximum amplitude
  73   short maxamp; //Maximum amplitude
  74   int index = SelectBunch( bunchvector,  &maxampindex,  &maxamp );
  75
  76   if( index >= 0)
  77     {
  78       Float_t  ped  = ReverseAndSubtractPed( &(bunchvector.at(index))  ,  altrocfg1, altrocfg2, fReversed  );
  79       int first;
  80       int last;
  81       Float_t maxf = TMath::MaxElement( bunchvector.at(index).GetLength(),  fReversed );
  82
  83       if ( maxf > fAmpCut )
  84         {
  85           short maxrev = maxampindex  -  bunchvector.at(index).GetStartBin();
  86           short timebinOffset = maxampindex - (bunchvector.at(index).GetLength()-1);
  87           SelectSubarray( fReversed,  bunchvector.at(index).GetLength(),  maxrev, &first, &last);
  88           int nsamples =  last - first + 1;
  89
  90           if( ( nsamples  )  > fNsampleCut )
  91             {
  92
  93               TGraph *graph =  new TGraph(  nsamples, fXaxis,  &fReversed[first] );
  94               fTf1->SetParameter(0, maxf*fkEulerSquared );
  95               fTf1->SetParameter(1, 0.2);
  96               //              fTf1->SetParameter(2,  2);
  97               //              fTf1->SetParameter(2,  3);
  98
  99               //     return   AliCaloFitResults( -1 , -1 , -1, -1, -1, -1 , -1);
 100
 101               Short_t tmpStatus =  graph->Fit(fTf1, "Q0RW");
 102
 103               //       return   AliCaloFitResults( -1 , -1 , -1, -1, -1, -1 , -1);
 104
 105               if( fVerbose == true )
 106                 {
 107                   AliCaloRawAnalyzer::PrintBunch( bunchvector.at(index) );
 108                   PrintFitResult( fTf1 ) ;
 109                 }
 110
 111                 delete graph;
 112                 return AliCaloFitResults( maxamp, ped ,    tmpStatus,
 113                                          fTf1->GetParameter(0)/fkEulerSquared,
 114                                          fTf1->GetParameter(1) + timebinOffset,
 115                                          fTf1->GetChisquare(),
 116                                          fTf1->GetNDF());
 117
 118
 119                 //     delete graph;
 120
 121             }
 122           else
 123             {
 124               return AliCaloFitResults( maxamp, ped, -1, maxf, -1, -1, -1 );
 125             }
 126         }
 127       else
 128         {
 129           return AliCaloFitResults( maxamp , ped, -1, maxf, -1, -1, -1 );
 130         }
 131     }
 132   else
 133     {
 134       return AliCaloFitResults( -99, -99 );
 135     }
 136
 137   return AliCaloFitResults( -99, -99 );
 138
 139 }
 140
 141
 142 void
 143 AliCaloRawAnalyzerLMS::PrintFitResult(const TF1 *f) const
 144 {
 145   //comment
 146   cout << endl;
 147   cout << __FILE__ << __LINE__ << "Using this samplerange we get" << endl;
 148   cout << __FILE__ << __LINE__ << "AMPLITUDE = " << f->GetParameter(0)/fkEulerSquared << ",.. !!!!" << endl;
 149   cout << __FILE__ << __LINE__ << "TOF = " << f->GetParameter(1) << ",.. !!!!" << endl;
 150   cout << __FILE__ << __LINE__ << "NDF = " << f->GetNDF() << ",.. !!!!" << endl;
 151   //  cout << __FILE__ << __LINE__ << "STATUS = " << f->GetStatus()  << ",.. !!!!" << endl << endl;
 152   cout << endl << endl;
 153 }