diff --git a/src/THcShowerArray.cxx b/src/THcShowerArray.cxx index 57a231f2a078c8b7604317d903a3e08bfbae08ed..39944a40f4d78b3e4a67b51fd8479d22e8f6aa89 100644 --- a/src/THcShowerArray.cxx +++ b/src/THcShowerArray.cxx @@ -56,11 +56,23 @@ THcShowerArray::THcShowerArray( const char* name, THcShowerArray::~THcShowerArray() { // Destructor - delete fXPos; - delete fYPos; - delete fZPos; + + for (UInt_t i=0; i<fNRows; i++) { + delete [] fXPos[i]; + delete [] fYPos[i]; + delete [] fZPos[i]; + } + + delete [] fPedLimit; + delete [] fGain; + delete [] fPedSum; + delete [] fPedSum2; + delete [] fPedCount; + delete [] fSig; + delete [] fPed; + delete [] fThresh; - delete fADCHits; + delete fADCHits; fADCHits = NULL; delete frAdcPedRaw; frAdcPedRaw = NULL; delete frAdcErrorFlag; frAdcErrorFlag = NULL; @@ -72,11 +84,11 @@ THcShowerArray::~THcShowerArray() delete frAdcPulseInt; frAdcPulseInt = NULL; delete frAdcPulseAmp; frAdcPulseAmp = NULL; - delete [] fA; - delete [] fP; - delete [] fA_p; + // delete [] fA; + //delete [] fP; + // delete [] fA_p; - delete [] fE; + //delete [] fE; delete [] fBlock_ClusterID; } @@ -135,9 +147,13 @@ Int_t THcShowerArray::ReadDatabase( const TDatime& date ) {"cal_ped_sample_high", &fPedSampHigh, kInt, 0, 1}, {"cal_data_sample_low", &fDataSampLow, kInt, 0, 1}, {"cal_data_sample_high", &fDataSampHigh, kInt, 0, 1}, + {"cal_debug_adc", &fDebugAdc, kInt, 0, 1}, {0} }; - gHcParms->LoadParmValues((DBRequest*)&list, prefix); + + fDebugAdc = 0; // Set ADC debug parameter to false unless set in parameter file + + gHcParms->LoadParmValues((DBRequest*)&list, prefix); fADCMode=kADCDynamicPedestal; fAdcTimeWindowMin=0; fAdcTimeWindowMax=10000; @@ -307,14 +323,24 @@ Int_t THcShowerArray::ReadDatabase( const TDatime& date ) InitializePedestals(); // Event by event amplitude and pedestal - fA = new Double_t[fNelem]; - fP = new Double_t[fNelem]; - fA_p = new Double_t[fNelem]; + //fA = new Double_t[fNelem]; + //fP = new Double_t[fNelem]; + //fA_p = new Double_t[fNelem]; + + fE = vector<Double_t> (fNelem, 0.0); + fNumGoodAdcHits = vector<Int_t> (fNelem, 0.0); + fGoodAdcPulseIntRaw = vector<Double_t> (fNelem, 0.0); + fGoodAdcPed = vector<Double_t> (fNelem, 0.0); + fGoodAdcPulseInt = vector<Double_t> (fNelem, 0.0); + fGoodAdcPulseAmp = vector<Double_t> (fNelem, 0.0); + fGoodAdcPulseTime = vector<Double_t> (fNelem, 0.0); + + fBlock_ClusterID = new Int_t[fNelem]; // Energy depositions per block. - fE = new Double_t[fNelem]; + //fE = new Double_t[fNelem]; #ifdef HITPIC hitpic = new char*[fNRows]; @@ -345,35 +371,49 @@ Int_t THcShowerArray::DefineVariables( EMode mode ) fIsSetup = ( mode == kDefine ); // Register variables in global list - RVarDef vars[] = { - {"adchits", "List of ADC hits", "fADCHits.THcSignalHit.GetPaddleNumber()"}, - {"a", "Raw ADC Amplitude", "fA"}, - {"p", "Dynamic ADC Pedestal", "fP"}, - {"a_p", "Sparsified, ped-subtracted ADC Amplitudes", "fA_p"}, - { "nhits", "Number of hits", "fNhits" }, - { "nghits", "Number of good hits ( pass threshold on raw ADC)", "fNgoodhits" }, - { "nclust", "Number of clusters", "fNclust" }, - {"e", "Energy Depositions per block", "fE"}, - {"block_clusterID", "Cluster ID number", "fBlock_ClusterID"}, - {"earray", "Energy Deposition in array", "fEarray"}, - { "ntracks", "Number of shower tracks", "fNtracks" }, - - {"adcCounter", "List of ADC counter numbers.", "frAdcPulseIntRaw.THcSignalHit.GetPaddleNumber()"}, - - {"adcPedRaw", "List of raw ADC pedestals", "frAdcPedRaw.THcSignalHit.GetData()"}, - {"adcErrorFlag", "List of raw ADC pedestals", "frAdcErrorFlag.THcSignalHit.GetData()"}, - {"adcPulseIntRaw", "List of raw ADC pulse integrals.", "frAdcPulseIntRaw.THcSignalHit.GetData()"}, - {"adcPulseAmpRaw", "List of raw ADC pulse amplitudes.", "frAdcPulseAmpRaw.THcSignalHit.GetData()"}, - {"adcPulseTimeRaw", "List of raw ADC pulse times.", "frAdcPulseTimeRaw.THcSignalHit.GetData()"}, - - {"adcPed", "List of ADC pedestals", "frAdcPed.THcSignalHit.GetData()"}, - {"adcPulseInt", "List of ADC pulse integrals.", "frAdcPulseInt.THcSignalHit.GetData()"}, - {"adcPulseAmp", "List of ADC pulse amplitudes.", "frAdcPulseAmp.THcSignalHit.GetData()"}, - - { 0 } - }; + vector<RVarDef> vars; + + //{"adchits", "List of ADC hits", "fADCHits.THcSignalHit.GetPaddleNumber()"}, // appears an empty histogram in the root file + + vars.push_back(RVarDef{"adcErrorFlag", "Error Flag When FPGA Fails", "frAdcErrorFlag.THcSignalHit.GetData()"}); + + vars.push_back(RVarDef{"adcCounter", "List of ADC counter numbers.", "frAdcPulseIntRaw.THcSignalHit.GetPaddleNumber()"}); //raw occupancy + vars.push_back(RVarDef{"numGoodAdcHits", "Number of Good ADC Hits per PMT", "fNumGoodAdcHits" }); //good occupancy + + vars.push_back(RVarDef{"totNumAdcHits", "Total Number of ADC Hits", "fTotNumAdcHits" }); // raw multiplicity + vars.push_back(RVarDef{"totNumGoodAdcHits", "Total Number of Good ADC Hits", "fTotNumGoodAdcHits" }); // good multiplicity + + + vars.push_back(RVarDef{"goodAdcPulseIntRaw", "Good Raw ADC Pulse Integrals", "fGoodAdcPulseIntRaw"}); //this is defined as pulseIntRaw, NOT ADC Amplitude in FillADC_DynamicPedestal() method + vars.push_back(RVarDef{"goodAdcPed", "Good ADC Pedestals", "fGoodAdcPed"}); + vars.push_back(RVarDef{"goodAdcPulseInt", "Good ADC Pulse Integrals", "fGoodAdcPulseInt"}); //this is defined as pulseInt, which is the pedestal subtracted pulse integrals, and is defined if threshold is passed + vars.push_back(RVarDef{"goodAdcPulseAmp", "Good ADC Pulse Amplitudes", "fGoodAdcPulseAmp"}); + vars.push_back(RVarDef{"goodAdcPulseTime", "Good ADC Pulse Times", "fGoodAdcPulseTime"}); //this is defined as pulseInt, which is the pedestal subtracted pulse integrals, and is defined if threshold is passed + + + vars.push_back(RVarDef{"e", "Energy Depositions per block", "fE"}); //defined as fE = fA_p*fGain = pulseInt * Gain + vars.push_back(RVarDef{"earray", "Energy Deposition in Shower Array", "fEarray"}); //defined as a Double_t and represents a sum of the total deposited energy in the shower counter + + vars.push_back(RVarDef{"nclust", "Number of clusters", "fNclust" }); //what is the difference between nclust defined here and that in THcShower.cxx ? + vars.push_back(RVarDef{"block_clusterID", "Cluster ID number", "fBlock_ClusterID"}); // im NOT very clear about this. it is histogrammed at wither -1 or 0. + vars.push_back(RVarDef{"ntracks", "Number of shower tracks", "fNtracks" }); //number of cluster-to-track associations + + if (fDebugAdc) { + vars.push_back(RVarDef{"adcPedRaw", "List of raw ADC pedestals", "frAdcPedRaw.THcSignalHit.GetData()"}); + vars.push_back(RVarDef{"adcPulseIntRaw", "List of raw ADC pulse integrals.", "frAdcPulseIntRaw.THcSignalHit.GetData()"}); + vars.push_back(RVarDef{"adcPulseAmpRaw", "List of raw ADC pulse amplitudes.", "frAdcPulseAmpRaw.THcSignalHit.GetData()"}); + vars.push_back(RVarDef{"adcPulseTimeRaw", "List of raw ADC pulse times.", "frAdcPulseTimeRaw.THcSignalHit.GetData()"}); + + vars.push_back(RVarDef{"adcPed", "List of ADC pedestals", "frAdcPed.THcSignalHit.GetData()"}); + vars.push_back(RVarDef{"adcPulseInt", "List of ADC pulse integrals.", "frAdcPulseInt.THcSignalHit.GetData()"}); + vars.push_back(RVarDef{"adcPulseAmp", "List of ADC pulse amplitudes.", "frAdcPulseAmp.THcSignalHit.GetData()"}); + + } - return DefineVarsFromList( vars, mode ); + RVarDef end {0}; + vars.push_back(end); + + return DefineVarsFromList(vars.data(), mode ); } //_____________________________________________________________________________ @@ -382,8 +422,8 @@ void THcShowerArray::Clear( Option_t* ) // Clears the hit lists fADCHits->Clear(); - fNhits = 0; - fNgoodhits = 0; + fTotNumAdcHits = 0; + fTotNumGoodAdcHits = 0; fNclust = 0; fClustSize = 0; fNtracks = 0; @@ -408,6 +448,16 @@ void THcShowerArray::Clear( Option_t* ) frAdcPulseInt->Clear(); frAdcPulseAmp->Clear(); + for (UInt_t ielem = 0; ielem < fGoodAdcPed.size(); ielem++) { + fGoodAdcPulseIntRaw.at(ielem) = 0.0; + fGoodAdcPed.at(ielem) = 0.0; + fGoodAdcPulseInt.at(ielem) = 0.0; + fGoodAdcPulseAmp.at(ielem) = 0.0; + fGoodAdcPulseTime.at(ielem) = 0.0; + fNumGoodAdcHits.at(ielem) = 0.0; + fE.at(ielem) = 0.0; + } + } //_____________________________________________________________________________ @@ -424,7 +474,6 @@ Int_t THcShowerArray::CoarseProcess( TClonesArray& tracks ) // Fill set of unclustered shower array hits. // Reuse hit class pertained to the HMS/SOS type calorimeters. - // Save Y coordinate of the hit in Z parameter of the class. // Save energy deposition in the module as hit mean energy, do not use // positive and negative side energies. @@ -434,7 +483,7 @@ Int_t THcShowerArray::CoarseProcess( TClonesArray& tracks ) for(UInt_t j=0; j < fNColumns; j++) { for (UInt_t i=0; i<fNRows; i++) { - if (fA_p[k] > 0) { //hit + if (fGoodAdcPulseInt.at(k) > 0) { //hit THcShowerHit* hit = new THcShowerHit(i, j, fXPos[i][j], fYPos[i][j], fZPos[i][j], fE[k], 0., 0.); @@ -454,13 +503,22 @@ Int_t THcShowerArray::CoarseProcess( TClonesArray& tracks ) cout << "Debug output from THcShowerArray::CoarseProcess for " << GetName() << endl; - cout << " List of unclustered hits. Total hits: " << fNhits << endl; + cout << " List of unclustered hits. Total hits: " << fTotNumAdcHits << endl; THcShowerHitIt it = HitSet.begin(); //<set> version - for (Int_t i=0; i!=fNgoodhits; i++) { + for (Int_t i=0; i!=fTotNumGoodAdcHits; i++) { cout << " hit " << i << ": "; (*(it++))->show(); } } + + ////Sanity check. (Vardan) + + // if ((int)HitSet.size() != fTotNumGoodAdcHits) { + // cout << "***" << endl; + // cout << "*** THcShowerArray::CoarseProcess: HitSet.size = " << HitSet.size() + // << " != fTotNumGoodAdcHits = " << fTotNumGoodAdcHits << endl; + // cout << "***" << endl; + // } // Cluster hits and fill list of clusters. @@ -718,7 +776,7 @@ Int_t THcShowerArray::CoarseProcessHits() Int_t k=0; for(UInt_t j=0; j < fNColumns; j++) { for (UInt_t i=0; i<fNRows; i++) { - if(fA[k] > fThresh[k]) { + if(fGoodAdcPulseIntRaw.at(k) > fThresh[k]) { cout << " counter = " << k << " E = " << fE[k] << endl; @@ -763,19 +821,30 @@ void THcShowerArray::FillADC_DynamicPedestal() { for (Int_t ielem=0;ielem<frAdcPulseInt->GetEntries();ielem++) { Int_t npad = ((THcSignalHit*) frAdcPulseInt->ConstructedAt(ielem))->GetPaddleNumber() - 1; - Double_t pulseInt = ((THcSignalHit*) frAdcPulseInt->ConstructedAt(ielem))->GetData(); Double_t pulseIntRaw = ((THcSignalHit*) frAdcPulseIntRaw->ConstructedAt(ielem))->GetData(); + Double_t pulsePed = ((THcSignalHit*) frAdcPed->ConstructedAt(ielem))->GetData(); + Double_t pulseInt = ((THcSignalHit*) frAdcPulseInt->ConstructedAt(ielem))->GetData(); + Double_t pulseAmp = ((THcSignalHit*) frAdcPulseAmp->ConstructedAt(ielem))->GetData(); Double_t pulseTime = ((THcSignalHit*) frAdcPulseTimeRaw->ConstructedAt(ielem))->GetData(); - Double_t errorflag = ((THcSignalHit*) frAdcErrorFlag->ConstructedAt(ielem))->GetData(); + Bool_t errorflag = ((THcSignalHit*) frAdcErrorFlag->ConstructedAt(ielem))->GetData(); Bool_t pulseTimeCut = (pulseTime > fAdcTimeWindowMin) && (pulseTime < fAdcTimeWindowMax); - if (errorflag==0 && pulseTimeCut) { - fNhits++; - fA[npad] =pulseIntRaw; - if(fA[npad] > fThresh[npad]) { - fNgoodhits++; - fA_p[npad] =pulseInt ; - fE[npad] = fA_p[npad]*fGain[npad]; - fEarray += fE[npad]; + if (!errorflag && pulseTimeCut) { + fTotNumAdcHits++; + fGoodAdcPulseIntRaw.at(npad) = pulseIntRaw; + + if(fGoodAdcPulseIntRaw.at(npad) > fThresh[npad]) { + fTotNumGoodAdcHits++; + fGoodAdcPulseInt.at(npad) = pulseInt; + fE.at(npad) = fGoodAdcPulseInt.at(npad)*fGain[npad]; + fEarray += fE.at(npad); + + fGoodAdcPed.at(npad) = pulsePed; + fGoodAdcPulseAmp.at(npad) = pulseAmp; + fGoodAdcPulseTime.at(npad) = pulseTime; + + fNumGoodAdcHits.at(npad) = npad + 1; + + } } } @@ -804,9 +873,10 @@ Int_t THcShowerArray::ProcessHits(TClonesArray* rawhits, Int_t nexthit) frAdcPulseAmp->Clear(); for(Int_t i=0;i<fNelem;i++) { - fA[i] = 0; - fA_p[i] = 0; - fE[i] = 0; + //fA[i] = 0; + //fA_p[i] = 0; + //fE[i] = 0; + fBlock_ClusterID[i] = -1; } @@ -854,7 +924,7 @@ Int_t THcShowerArray::ProcessHits(TClonesArray* rawhits, Int_t nexthit) } #if 0 - if(fNgoodhits > 0) { + if(fTotNumGoodAdcHits > 0) { cout << "+"; for(Int_t column=0;column<fNColumns;column++) { cout << "-"; @@ -864,7 +934,7 @@ Int_t THcShowerArray::ProcessHits(TClonesArray* rawhits, Int_t nexthit) cout << "|"; for(Int_t column=0;column<fNColumns;column++) { Int_t counter = column*fNRows + row; - if(fA[counter]>threshold) { + if(fGoodAdcPulseIntRaw.at(counter) > threshold) { cout << "X"; } else { cout << " "; @@ -875,14 +945,14 @@ Int_t THcShowerArray::ProcessHits(TClonesArray* rawhits, Int_t nexthit) } #endif #ifdef HITPIC - if(fNgoodhits > 0) { + if(fTotNumGoodAdcHits > 0) { for(Int_t row=0;row<fNRows;row++) { if(piccolumn==0) { hitpic[row][0] = '|'; } for(Int_t column=0;column<fNColumns;column++) { Int_t counter = column*fNRows+row; - if(fA[counter] > threshold) { + if(fGoodAdcPulseIntRaw.at(counter) > threshold) { hitpic[row][piccolumn*(fNColumns+1)+column+1] = 'X'; } else { hitpic[row][piccolumn*(fNColumns+1)+column+1] = ' '; @@ -898,7 +968,7 @@ Int_t THcShowerArray::ProcessHits(TClonesArray* rawhits, Int_t nexthit) cout << "-"; } cout << "+"; - } +0 } cout << endl; for(Int_t row=0;row<fNRows;row++) { hitpic[row][(piccolumn+1)*(fNColumns+1)+1] = '\0'; diff --git a/src/THcShowerArray.h b/src/THcShowerArray.h index a7e2a0e4bcca8aa16ea80a17f5283737a05a74e9..8e45937cb1e1e15c7c38c007cc4a7f3ea289638c 100644 --- a/src/THcShowerArray.h +++ b/src/THcShowerArray.h @@ -95,11 +95,27 @@ protected: char **hitpic; Int_t piccolumn; #endif - Double_t* fA; // [fNelem] ADC amplitudes of blocks - Double_t* fP; // [fNelem] Event by event (FADC) pedestals - Double_t* fA_p; // [fNelem] sparsified, pedestal subtracted - // (FASTBUS) ADC amplitudes + + + + + //counting variables + Int_t fTotNumAdcHits; // Total number of ADC hits + Int_t fTotNumGoodAdcHits; // Total number of good ADC hits (pass threshold) + + vector<Int_t> fNumGoodAdcHits; // shower good occupancy + vector<Double_t> fGoodAdcPulseIntRaw; // [fNelem] Good Raw ADC pulse Integrals of blocks + + vector<Double_t> fGoodAdcPed; // [fNelem] Event by event (FADC) good pulse pedestals + vector<Double_t> fGoodAdcPulseInt; // [fNelem] good pedestal subtracted pulse integrals + vector<Double_t> fGoodAdcPulseAmp; + vector<Double_t> fGoodAdcPulseTime; + + vector<Double_t> fE; //[fNelem] energy deposition in shower blocks + + Int_t* fBlock_ClusterID; // [fNelem] Cluster ID of the block -1 then not in a cluster + Double_t fEarray; // Total Energy deposition in the array. TClonesArray* fADCHits; // List of ADC hits // Parameters @@ -128,7 +144,9 @@ protected: Double_t fAdcTimeWindowMin ; Double_t fAdcTimeWindowMax ; Double_t fAdcThreshold ; -Int_t fPedSampLow; // Sample range for + + Int_t fDebugAdc; + Int_t fPedSampLow; // Sample range for Int_t fPedSampHigh; // dynamic pedestal Int_t fDataSampLow; // Sample range for Int_t fDataSampHigh; // sample integration @@ -152,15 +170,7 @@ Int_t fPedSampLow; // Sample range for Float_t *fThresh; // [fNelem] ADC thresholds Double_t* fGain; // [fNelem] Gain constants from calibration - - //Energy depositions. - - Double_t* fE; // [fNelem] energy depositions in the blocks. - Int_t* fBlock_ClusterID; // [fNelem] Cluster ID of the block -1 then not in a cluster - Double_t fEarray; // Total Energy deposition in the array. - - Int_t fNhits; // Total number of hits - Int_t fNgoodhits; // Total number of good hits (pass threshold) + Int_t fNclust; // Number of hit clusters Int_t fNtracks; // Number of shower tracks, i.e. number of // cluster-to-track associations