update from Per Thomas - new classes for NeuralNet and FastFit + code warning fixes
[u/mrichter/AliRoot.git] / EMCAL / AliCaloRawAnalyzer.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 <perthomas.hille@yale.edu>            *
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 // Base class for extraction 
21 // of signal amplitude and peak position
22 // From CALO Calorimeter RAW data (from the RCU)
23 // Contains some utilities for preparing / selecting
24 // Signals suitable for signal extraction
25 // By derived classes
26
27 #include "AliLog.h"
28 #include "AliCaloRawAnalyzer.h"
29 #include "AliCaloBunchInfo.h"
30 #include "AliCaloFitResults.h"
31 #include "TMath.h"
32 #include <iostream>
33 using namespace std;
34
35 ClassImp(AliCaloRawAnalyzer)  
36
37 AliCaloRawAnalyzer::AliCaloRawAnalyzer(const char *name) :  TObject(),
38                                                             fMinTimeIndex(-1),
39                                                             fMaxTimeIndex(-1),
40                                                             fFitArrayCut(5),
41                                                             fAmpCut(4),
42                                                             fNsampleCut(5),
43                                                             fIsZerosupressed( false ),
44                                                             fVerbose( false )
45 {
46   //Comment 
47   sprintf(fName, "%s", name);
48   for(int i=0; i < MAXSAMPLES; i++ )
49     {
50       fReversed[i] = 0;
51     }
52 }
53
54 AliCaloRawAnalyzer::~AliCaloRawAnalyzer()
55 {
56
57 }
58
59
60 void 
61 AliCaloRawAnalyzer::SetTimeConstraint(const int min, const int max ) 
62 {
63   //Require that the bin if the maximum ADC value is between min and max (timebin)
64   if(  ( min > max ) || min > MAXSAMPLES  || max > MAXSAMPLES  )
65     {
66       AliWarning( Form( "Attempt to set Invalid time bin range (Min , Max) = (%d, %d), Ingored",  min, max ) ); 
67     }
68   else
69     {
70       fMinTimeIndex  = min;
71       fMaxTimeIndex  = max;
72     }
73 }
74
75
76 UShort_t 
77 AliCaloRawAnalyzer::Max(const UShort_t *data, const int length ) const
78 {
79   //------------
80   UShort_t tmpmax  = data[0];
81
82   for(int i=0; i < length; i++)
83     {
84       if( tmpmax  <  data[i] )
85         {
86           tmpmax = data[i];
87         }
88     }
89   return tmpmax;
90 }
91
92
93 void 
94 AliCaloRawAnalyzer::SelectSubarray( const Double_t *fData, const int length, const short maxindex, int *const first,  int *const last ) const
95 {
96   //Selection of subset of data from one bunch that will be used for fitting or
97   //Peak finding.  Go to the left and right of index of the maximum time bin
98   //Untile the ADC value is less that fFitArrayCut
99   int tmpfirst =  maxindex;
100   int tmplast  =  maxindex;
101   
102   while((  tmpfirst  ) > 0  &&  ( fData[tmpfirst] >  fFitArrayCut   ))  
103     {
104       tmpfirst -- ;
105     }
106   
107   while(( tmplast ) <  length   && ( fData [tmplast] >  fFitArrayCut ))
108     {
109       tmplast ++;
110     }
111   
112   *first = tmpfirst +1;
113   *last =  tmplast -1;
114 }
115
116
117
118 Float_t 
119 AliCaloRawAnalyzer::ReverseAndSubtractPed( const AliCaloBunchInfo *bunch, const UInt_t /*altrocfg1*/,  const UInt_t /*altrocfg2*/, double *outarray ) const
120 {
121   //Time sample comes in reversed order, revers them back
122   //Subtract the baseline based on content of altrocfg1 and altrocfg2.
123   Int_t length =  bunch->GetLength(); 
124   const UShort_t *sig = bunch->GetData();
125
126   double ped = EvaluatePedestal( sig , length);
127
128   for( int i=0; i < length; i++ )
129     {
130       outarray[i] = sig[length -i -1] - ped;
131     }
132
133   return ped;
134 }
135
136
137
138 Float_t 
139 AliCaloRawAnalyzer::EvaluatePedestal(const UShort_t * const data, const int /*length*/ ) const
140 {
141   //  double ped = 0;
142   double tmp = 0;
143
144   if( fIsZerosupressed == false ) 
145     {
146       for(int i=0; i < 5; i++ )
147         {
148           tmp += data[i];
149         }
150     }
151
152   //  cout << __FILE__ << __LINE__ << "XXXXXXXXXXX returning " <<   tmp/5 << endl;
153
154   return  tmp/5;
155 }
156
157
158 short  
159 AliCaloRawAnalyzer::Max( const AliCaloBunchInfo *const bunch , int *const maxindex ) const
160 {
161   //comment
162   short tmpmax = -1;
163   int tmpindex = -1;
164   const UShort_t *sig = bunch->GetData();
165
166   for(int i=0; i < bunch->GetLength(); i++ )
167     {
168       if( sig[i] > tmpmax )
169         {
170           tmpmax   =  sig[i];
171           tmpindex =  i; 
172         }
173     }
174   
175   if(maxindex != 0 )
176     {
177       *maxindex =  bunch->GetLength() -1 - tmpindex + bunch->GetStartBin(); 
178     }
179   
180   return  tmpmax;
181 }
182
183
184 int 
185 AliCaloRawAnalyzer::SelectBunch( const vector<AliCaloBunchInfo> &bunchvector,short *const maxampbin, short *const maxamplitude ) const
186 {
187   //We select the bunch with the highest amplitude unless any time constraints is set
188   short max = -1;
189   short bunchindex = -1;
190   short maxall = -1;
191   int indx = -1;
192
193   for(unsigned int i=0; i < bunchvector.size(); i++ )
194     {
195       max = Max(  &bunchvector.at(i), &indx );  
196       if( IsInTimeRange( indx) )
197         {
198           if( max > maxall )
199             {
200               maxall = max;
201               bunchindex = i;
202             }
203         }
204     }
205  
206   *maxampbin     = indx;
207   *maxamplitude  = max;
208   return  bunchindex;
209 }
210
211
212 bool 
213 AliCaloRawAnalyzer::IsInTimeRange( const int maxindex  ) const
214 {
215   // Ckeck if the index of the max ADC vaue is consistent with trigger.
216   if( ( fMinTimeIndex  < 0 && fMaxTimeIndex  < 0) ||fMaxTimeIndex  < 0 )
217     {
218       return true; 
219     }
220   return ( maxindex < fMaxTimeIndex ) && ( maxindex > fMinTimeIndex  ) ? true : false;
221 }
222
223
224 void 
225 AliCaloRawAnalyzer::PrintBunches( const vector<AliCaloBunchInfo> &bvctr ) const
226 {
227   //comment
228   cout << __FILE__ << __LINE__<< "*************** Printing Bunches *******************" << endl;
229   cout << __FILE__ << __LINE__<< "*** There are " << bvctr.size() << ", bunches" << endl;
230
231   for(unsigned int i=0; i < bvctr.size() ; i++ )
232     {
233       PrintBunch( bvctr.at(i) );
234       cout << " bunch = "  <<  i  << endl;
235     }
236   cout << __FILE__ << __LINE__<< "*************** Done ... *******************" << endl;
237 }
238
239
240 void 
241 AliCaloRawAnalyzer::PrintBunch( const AliCaloBunchInfo &bunch ) const
242 {
243   //comment
244   cout << __FILE__ << ":" << __LINE__ << endl;
245   cout << __FILE__ << __LINE__   << ", startimebin = " << bunch.GetStartBin() << ", length =" <<  bunch.GetLength() << endl;
246   const UShort_t *sig =  bunch.GetData();  
247   
248   for ( Int_t j = 0; j <  bunch.GetLength();  j++) 
249     {
250       printf("%d\t", sig[j] );
251     }
252   cout << endl; 
253 }
254
255 AliCaloFitResults
256 AliCaloRawAnalyzer::Evaluate( const vector<AliCaloBunchInfo>  &/*bunchvector*/, const UInt_t /*altrocfg1*/,  const UInt_t /*altrocfg2*/)
257 { // method to do the selection of what should possibly be fitted
258   // not implemented for base class
259   return AliCaloFitResults( 0, 0, 0, 0, 0, 0, 0 );
260 }
261
262 int
263 AliCaloRawAnalyzer::PreFitEvaluateSamples( const vector<AliCaloBunchInfo>  &bunchvector, const UInt_t altrocfg1,  const UInt_t altrocfg2, Int_t & index, Float_t & maxf, short & maxamp, short & maxampindex, Float_t & ped, int & first, int & last)
264 { // method to do the selection of what should possibly be fitted
265   int nsamples  = 0;
266   index = SelectBunch( bunchvector,  &maxampindex,  &maxamp ); // select the bunch with the highest amplitude unless any time constraints is set
267
268   
269   if( index >= 0 && maxamp > fAmpCut) // something valid was found, and non-zero amplitude
270     {
271       // use more convenient numbering and possibly subtract pedestal
272       ped  = ReverseAndSubtractPed( &(bunchvector.at(index)),  altrocfg1, altrocfg2, fReversed  );
273       maxf = TMath::MaxElement( bunchvector.at(index).GetLength(),  fReversed );
274       
275       if ( maxf > fAmpCut ) // possibly significant signal
276         {
277           // select array around max to possibly be used in fit
278           maxampindex -= bunchvector.at(index).GetStartBin(); // PTH - why isn't this index also reversed for call below?
279           SelectSubarray( fReversed,  bunchvector.at(index).GetLength(),  maxampindex , &first, &last);
280
281           // sanity check: maximum should not be in first or last bin
282           // if we should do a fit
283           if (first!=maxampindex && last!=maxampindex) {
284             // calculate how many samples we have 
285             nsamples =  last - first + 1;         
286           }
287         }
288     }
289
290   return nsamples;
291 }
292