Newer
Older
///////////////////////////////////////////////////////////////////////////////////////
// //
// THcCherenkov //
// //
// Zafar Ahmed. Fourth pull request. January 30 2014. //
// Comment: No need to cahnge the Map file. //
// Comment New parameter hcer_tot_pmts = 2 is added to the hcana.param file //
// This code is for the coarse process of THcCherenkov class. //
// //
// Variable Name Description //
// //
// phototubes Nuber of Cherenkov photo tubes //
// adc Raw ADC values //
// adc_p Pedestal Subtracted ADC values //
// npe Number of Photo electrons //
// npesum Sum of Number of Photo electrons //
// ncherhit Number of Hits(Cherenkov) //
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
// //
///////////////////////////////////////////////////////////////////////////////////////
#include "THcCherenkov.h"
#include "TClonesArray.h"
#include "THcSignalHit.h"
#include "THaEvData.h"
#include "THaDetMap.h"
#include "THcDetectorMap.h"
#include "THcGlobals.h"
#include "THaCutList.h"
#include "THcParmList.h"
#include "THcHitList.h"
#include "VarDef.h"
#include "VarType.h"
#include "THaTrack.h"
#include "TClonesArray.h"
#include "TMath.h"
#include "THaTrackProj.h"
#include <cstring>
#include <cstdio>
#include <cstdlib>
#include <iostream>
using namespace std;
using std::cout;
using std::cin;
using std::endl;
#include <iomanip>
using std::setw;
using std::setprecision;
//_____________________________________________________________________________
THcCherenkov::THcCherenkov( const char* name, const char* description,
THaApparatus* apparatus ) :
THaNonTrackingDetector(name,description,apparatus)
{
// Normal constructor with name and description
cout << "fADCHits " << fADCHits << endl;
InitArrays();
cout << "fADCHits " << fADCHits << endl;
}
//_____________________________________________________________________________
THcCherenkov::THcCherenkov( ) :
THaNonTrackingDetector()
{
// Constructor
fADCHits = NULL;
InitArrays();
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
//_____________________________________________________________________________
void THcCherenkov::InitArrays()
{
fGain = NULL;
fCerWidth = NULL;
fNPMT = NULL;
fADC = NULL;
fADC_P = NULL;
fNPE = NULL;
fPedSum = NULL;
fPedSum2 = NULL;
fPedLimit = NULL;
fPedMean = NULL;
fPedCount = NULL;
fPed = NULL;
fThresh = NULL;
}
//_____________________________________________________________________________
void THcCherenkov::DeleteArrays()
{
delete [] fGain; fGain = NULL;
delete [] fCerWidth; fCerWidth = NULL;
delete [] fNPMT; fNPMT = NULL;
delete [] fADC; fADC = NULL;
delete [] fADC; fADC_P = NULL;
delete [] fNPE; fNPE = NULL;
delete [] fPedSum; fPedSum = NULL;
delete [] fPedSum2; fPedSum2 = NULL;
delete [] fPedLimit; fPedLimit = NULL;
delete [] fPedMean; fPedMean = NULL;
delete [] fPedCount; fPedCount = NULL;
delete [] fPed; fPed = NULL;
delete [] fThresh; fThresh = NULL;
}
//_____________________________________________________________________________
THcCherenkov::~THcCherenkov()
{
// Destructor
delete fADCHits; fADCHits = NULL;
DeleteArrays();
}
//_____________________________________________________________________________
THaAnalysisObject::EStatus THcCherenkov::Init( const TDatime& date )
{
static const char* const here = "Init()";
cout << "THcCherenkov::Init " << GetName() << endl;
// Should probably put this in ReadDatabase as we will know the
// maximum number of hits after setting up the detector map
Stephen A. Wood
committed
InitHitList(fDetMap, "THcCherenkovHit", 100); // 100 is max hits
EStatus status;
if( (status = THaNonTrackingDetector::Init( date )) )
return fStatus=status;
// Will need to determine which apparatus it belongs to and use the
// appropriate detector ID in the FillMap call
if( gHcDetectorMap->FillMap(fDetMap, "HCER") < 0 ) {
Error( Here(here), "Error filling detectormap for %s.",
"HCER");
return kInitError;
}
return fStatus = kOK;
}
//_____________________________________________________________________________
Int_t THcCherenkov::ReadDatabase( const TDatime& date )
{
// This function is called by THaDetectorBase::Init() once at the beginning
// of the analysis.
cout << "THcCherenkov::ReadDatabase " << GetName() << endl; // Ahmed
char prefix[2];
prefix[0]=tolower(GetApparatus()->GetName()[0]);
prefix[1]='\0';
strcpy(parname,prefix); // This is taken from
strcat(parname,"cer_tot_pmts"); // THcScintillatorPlane
fNelem = (Int_t)gHcParms->Find(parname)->GetValue(); // class.
// fNelem = 2; // Default if not defined
fNPMT = new Int_t[fNelem];
fADC = new Double_t[fNelem];
fADC_P = new Double_t[fNelem];
fNPE = new Double_t[fNelem];
fCerWidth = new Double_t[fNelem];
fGain = new Double_t[fNelem];
fPedLimit = new Int_t[fNelem];
fPedMean = new Double_t[fNelem];
{"cer_adc_to_npe", fGain, kDouble, (UInt_t) fNelem}, // Ahmed
{"cer_ped_limit", fPedLimit, kInt, (UInt_t) fNelem}, // Ahmed
{"cer_width", fCerWidth, kDouble, (UInt_t) fNelem}, // Ahmed
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
gHcParms->LoadParmValues((DBRequest*)&list,prefix);
fIsInit = true;
// Create arrays to hold pedestal results
InitializePedestals();
return kOK;
}
//_____________________________________________________________________________
Int_t THcCherenkov::DefineVariables( EMode mode )
{
// Initialize global variables for histogramming and tree
cout << "THcCherenkov::DefineVariables called " << GetName() << endl;
if( mode == kDefine && fIsSetup ) return kOK;
fIsSetup = ( mode == kDefine );
// Register variables in global list
// Do we need to put the number of pos/neg TDC/ADC hits into the variables?
// No. They show up in tree as Ndata.H.aero.postdchits for example
RVarDef vars[] = {
{"phototubes", "Nuber of Cherenkov photo tubes", "fNPMT"},
{"adc", "Raw ADC values", "fADC"},
{"adc_p", "Pedestal Subtracted ADC values", "fADC_P"},
{"npe", "Number of Photo electrons", "fNPE"},
{"npesum", "Sum of Number of Photo electrons", "fNPEsum"},
{"ncherhit", "Number of Hits(Cherenkov)", "fNCherHit"},
{ 0 }
};
return DefineVarsFromList( vars, mode );
}
//_____________________________________________________________________________
inline
void THcCherenkov::Clear(Option_t* opt)
{
// Clear the hit lists
// Clear Cherenkov variables from h_trans_cer.f
fNhits = 0; // Don't really need to do this. (Be sure this called before Decode)
fNPEsum = 0;
fNCherHit = 0;
for(Int_t itube = 0;itube < fNelem;itube++) {
fNPMT[itube] = 0;
fADC[itube] = 0;
fADC_P[itube] = 0;
fNPE[itube] = 0;
}
}
//_____________________________________________________________________________
Int_t THcCherenkov::Decode( const THaEvData& evdata )
{
// Get the Hall C style hitlist (fRawHitList) for this event
Stephen A. Wood
committed
fNhits = DecodeToHitList(evdata);
if(gHaCuts->Result("Pedestal_event")) {
AccumulatePedestals(fRawHitList);
fAnalyzePedestals = 1; // Analyze pedestals first normal events
return(0);
}
if(fAnalyzePedestals) {
CalculatePedestals();
fAnalyzePedestals = 0; // Don't analyze pedestals next event
}
Int_t ihit = 0;
while(ihit < fNhits) {
THcCherenkovHit* hit = (THcCherenkovHit *) fRawHitList->At(ihit);
THcSignalHit *sighit = (THcSignalHit*) fADCHits->ConstructedAt(nADCHits++);
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
sighit->Set(hit->fCounter, hit->fADC_pos);
}
ihit++;
}
return ihit;
}
//_____________________________________________________________________________
Int_t THcCherenkov::ApplyCorrections( void )
{
return(0);
}
//_____________________________________________________________________________
Int_t THcCherenkov::CoarseProcess( TClonesArray& ) //tracks
{
/*
------------------------------------------------------------------------------------------------------------------
h_trans_cer.f code:
hcer_num_hits = 0 <--- clear event
do tube=1,hcer_num_mirrors
hcer_npe(tube) = 0. <--- clear event
hcer_adc(tube) = 0. <--- clear event
enddo
hcer_npe_sum = 0. <--- clear event
do nhit = 1, hcer_tot_hits <--- loop over total hits. Very first line of this method
tube = hcer_tube_num(nhit) <--- tube is number of PMT on either side and it is this
line: Int_t npmt = hit->fCounter - 1
hcer_adc(tube) = hcer_raw_adc(nhit) - hcer_ped(tube) <--- This is done above:
fA_Pos_p[npmt] = hit->fADC_pos - fPosPedMean[npmt];
fA_Neg_p[npmt] = hit->fADC_neg - fNegPedMean[npmt];
if (hcer_adc(tube) .gt. hcer_width(tube)) then <--- This needs to convert in hcana
hcer_num_hits = hcer_num_hits + 1
hcer_tube_num(hcer_num_hits) = tube
hcer_npe(tube) = hcer_adc(tube) * hcer_adc_to_npe(tube)
hcer_npe_sum = hcer_npe_sum + hcer_npe(tube)
endif
enddo
------------------------------------------------------------------------------------------------------------------
*/
for(Int_t ihit=0; ihit < fNhits; ihit++) {
THcCherenkovHit* hit = (THcCherenkovHit *) fRawHitList->At(ihit); // nhit = 1, hcer_tot_hits
// Pedestal subtraction and gain adjustment
// An ADC value of less than zero occurs when that particular
// channel has been sparsified away and has not been read.
// The NPE for that tube will be assigned zero by this code.
// An ADC value of greater than 8192 occurs when the ADC overflows on
// an input that is too large. Tubes with this characteristic will
// be assigned NPE = 100.0.
Int_t npmt = hit->fCounter - 1; // tube = hcer_tube_num(nhit)
// Should probably check that npmt is in range
if ( ihit != npmt )
cout << "ihit != npmt " << endl;
fNPMT[npmt] = hit->fCounter;
fADC[npmt] = hit->fADC_pos;
fADC_P[npmt] = hit->fADC_pos - fPedMean[npmt];
if ( ( fADC_P[npmt] > fCerWidth[npmt] ) && ( hit->fADC_pos < 8000 ) ) {
fNPE[npmt] = fGain[npmt]*fADC_P[npmt];
fNCherHit ++;
} else if ( hit->fADC_pos > 8000 ) {
fNPE[npmt] = 100.0;
} else {
fNPE[npmt] = 0.0;
}
ApplyCorrections();
return 0;
}
//_____________________________________________________________________________
Int_t THcCherenkov::FineProcess( TClonesArray& tracks )
{
return 0;
}
//_____________________________________________________________________________
void THcCherenkov::InitializePedestals( )
{
fNPedestalEvents = 0;
fMinPeds = 500; // In engine, this is set in parameter file
fPedSum = new Int_t [fNelem];
fPedSum2 = new Int_t [fNelem];
fPedCount = new Int_t [fNelem];
fPed = new Double_t [fNelem];
fThresh = new Double_t [fNelem];
fPedSum[i] = 0;
fPedSum2[i] = 0;
fPedCount[i] = 0;
}
}
//_____________________________________________________________________________
void THcCherenkov::AccumulatePedestals(TClonesArray* rawhits)
{
// Extract data from the hit list, accumulating into arrays for
// calculating pedestals
Int_t nrawhits = rawhits->GetLast()+1;
Int_t ihit = 0;
while(ihit < nrawhits) {
THcCherenkovHit* hit = (THcCherenkovHit *) rawhits->At(ihit);
Int_t element = hit->fCounter - 1;
Int_t nadc = hit->fADC_pos;
if(nadc <= fPedLimit[element]) {
fPedSum[element] += nadc;
fPedSum2[element] += nadc*nadc;
fPedCount[element]++;
if(fPedCount[element] == fMinPeds/5) {
fPedLimit[element] = 100 + fPedSum[element]/fPedCount[element];
}
}
ihit++;
}
fNPedestalEvents++;
return;
}
//_____________________________________________________________________________
void THcCherenkov::CalculatePedestals( )
{
// Use the accumulated pedestal data to calculate pedestals
// Later add check to see if pedestals have drifted ("Danger Will Robinson!")
// cout << "Plane: " << fPlaneNum << endl;
for(Int_t i=0; i<fNelem;i++) {
// PMT tubes
fPed[i] = ((Double_t) fPedSum[i]) / TMath::Max(1, fPedCount[i]);
fThresh[i] = fPed[i] + 15;
// Just a copy for now, but allow the possibility that fXXXPedMean is set
// in a parameter file and only overwritten if there is a sufficient number of
// pedestal events. (So that pedestals are sensible even if the pedestal events were
// not acquired.)
if(fMinPeds > 0) {
if(fPedCount[i] > fMinPeds) {
fPedMean[i] = fPed[i];
}
}
}
// cout << " " << endl;
}
void THcCherenkov::Print( const Option_t* opt) const {
THaNonTrackingDetector::Print(opt);
// Print out the pedestals
cout << endl;
cout << "Cherenkov Pedestals" << endl;
// Ahmed
}
cout << endl;
}
ClassImp(THcCherenkov)
////////////////////////////////////////////////////////////////////////////////