|
//____________________________________________________________________ // // BrTpcHitModule takes the ouput clusters from BrTpcClusterModule and // BrTpcDeconvoluteClusterModule and converts them into BrTpcHit(s) // which are then appended to the outNode // // The method only supports weighted mean methods at the moment and // the error on the hit is estimated from a settable white noise // parameter. // Note: FV 5/7/02 I do not think the error estimate for timing bins // where a sqrt(adc) count is used. This is not the proper statistics // since the fluctuation is in total change , not on individual count. // The use of padresponse widths from a DB is better, but does depend on // total Charge in clsuter (that comes from initial number of p.e. // // The adcGain can be corrected before calculating avaraged in pad x // and time bin using database calibration constants. The action is all in the // weightedmean routines for this. //____________________________________________________________________ // // $Id: BrTpcHitModule.cxx,v 1.7 2002/08/15 17:25:09 videbaek Exp $ // $Author: videbaek $ // $Date: 2002/08/15 17:25:09 $ // $Copyright: (C) 2001 BRAHMS Collaboration <brahmlib@rhic.bnl.gov> // #ifndef BRAT_BrTpcHitModule #include "BrTpcHitModule.h" #endif #ifndef ROOT_TDirectory #include "TDirectory.h" #endif #ifndef BRAT_BrParameterDbManager #include "BrParameterDbManager.h" #endif #ifndef BRAT_BrTableNames #include "BrTableNames.h" #endif #ifndef BRAT_BrTpcHit #include "BrTpcHit.h" #endif #ifndef WIN32 #include <iostream> #else #include <iostream.h> #endif #ifndef BRAT_BrException #include "BrException.h" #endif #if !defined BRAT_BrCalibrationManager #include "BrCalibrationManager.h" #endif #if !defined BRAT_BrTpcCalibration #include "BrTpcCalibration.h" #endif //____________________________________________________________________ ClassImp(BrTpcHitModule); //____________________________________________________________________ BrTpcHitModule::BrTpcHitModule() { // Default constructor. DO NOT USE SetState(kSetup); fCalibration = 0; // set default parameters SetWhiteNoise(0); SetUseCalibratedWidths(); SetUsePadGain(); } //____________________________________________________________________ BrTpcHitModule::BrTpcHitModule(const Char_t* name, const Char_t* title) : BrModule(name, title) { // Named constructor SetState(kSetup); fCalibration = 0; // set default parameters SetWhiteNoise(); SetUseCalibratedWidths(); SetUsePadGain(); } //____________________________________________________________________ void BrTpcHitModule::DefineHistograms() { // Define the histograms used to monitor performance // This method is usually called by BrModule::Book(). // // Initialiser at job level if (GetState() != kInit) { Stop("DefineHistograms", "Must be called after Init"); return; } // Create new subdirectory for histograms TDirectory* saveDir = gDirectory; Char_t dirName[32]; sprintf(dirName,"TpcHit_%s",GetName()); TDirectory* histDir = gDirectory->mkdir(dirName); if (!histDir) { Warning("DefineHistograms","Could not create histogram subdirectory"); return; } histDir->cd(); // Make histograms here Char_t Histname[48]; Char_t HistTitle[48]; const Int_t nRows = fParams->GetNumberOfRows(); const Int_t padsPrRow = fParams->GetPadsprow(); const Int_t timeBuckets = fParams->GetNoBuckets(); const Float_t xMax = fParams->PadToIntrinsic(0); const Float_t xMin = fParams->PadToIntrinsic(padsPrRow); const Float_t yMin = fParams->TimeToIntrinsic(timeBuckets); const Float_t yMax = fParams->TimeToIntrinsic(0); sprintf(Histname,"hPadTime_%s",GetName()); sprintf(HistTitle,"Pad vs Time %s",GetName()); hPadTime = new TH2F(Histname, HistTitle, padsPrRow, 0, padsPrRow, timeBuckets, 0, timeBuckets); sprintf(Histname,"hPadRow_%s",GetName()); sprintf(HistTitle,"Pad vs Row %s",GetName()); hPadRow = new TH2F(Histname, HistTitle, padsPrRow, 0, padsPrRow, nRows, 1, nRows+1); sprintf(Histname,"hX_%s",GetName()); sprintf(HistTitle,"x position %s",GetName()); hX = new TH1F(Histname, HistTitle, 100, xMin, xMax); sprintf(Histname,"hY_%s",GetName()); sprintf(HistTitle,"y position %s",GetName()); hY = new TH1F(Histname, HistTitle, 100, yMin, yMax); sprintf(Histname,"hDx_%s",GetName()); sprintf(HistTitle,"x error %s",GetName()); hDx = new TH1F(Histname, HistTitle, 100, 0, 0.2); sprintf(Histname,"hDy_%s",GetName()); sprintf(HistTitle,"y error %s",GetName()); hDy = new TH1F(Histname, HistTitle, 100, 0, 0.2); gDirectory = saveDir; } //____________________________________________________________________ void BrTpcHitModule::Init() { // Initialiser at job level SetState(kInit); BrParameterDbManager* gParamDb = BrParameterDbManager::Instance(); fParams = (BrDetectorParamsTPC*) gParamDb->GetDetectorParameters("BrDetectorParamsTPC", GetName()); if(!fUsePadGain) return; try{ // getting calibration parameter element BrCalibrationManager* calMan = BrCalibrationManager::Instance(); if (!calMan) { Abort("Init", "could not get calibration manager"); return; } fCalibration = (BrTpcCalibration*)calMan->Register("BrTpcCalibration", GetName()); if (!fCalibration) { Abort("Init", "could not get calibration parameter element %s", GetName()); return; } BrCalibration::EAccessMode mode = BrCalibration::kRead; fCalibration->Use(BrTpcCalibration::kPadAdcGain); } catch(BrException* e) { cerr << "BrTpcDeconvolute::Init : " << *e << endl; e->Execute(); } } //____________________________________________________________________ void BrTpcHitModule::ConvertClustersToHits() { const Int_t nEntries = fClusterTable->GetEntries(); if(DebugLevel()>2) cout << "Number of clusters to convert : " << nEntries << endl; for(Int_t i = 0; i < nEntries; i++) { BrTpcCluster *cluster = (BrTpcCluster*)fClusterTable->At(i); FindHitWeightedMean(cluster); } Assert(fHitTable->GetEntries() == fClusterTable->GetEntries()); } //____________________________________________________________________ void BrTpcHitModule::Event(BrEventNode* inNode, BrEventNode* outNode) { // Main method. Process one event. // Gets the TPCCluster data from the inputNode // Extracts the hit position and the error on the hit position. // It then converts the information to the local hit position and // store the information as a BrTpcHit in the outputtable // SetState(kEvent); if(DebugLevel()>0) cout << Form("Entering BrTpcHitModule (%s) Event()", GetName()) << endl; fClusterTable = inNode->GetDataTable(Form("%s %s", BRTABLENAMES kTPCCluster, GetName())); if(fClusterTable == 0) { if(DebugLevel() > 0 ) Warning("Event", "BrDataTable %s %s not found in this event!", BRTABLENAMES kTPCCluster, GetName()); return; } //make ouput hit table fHitTable = new BrDataTable(Form("%s %s", BRTABLENAMES kTpcHit, GetName())); if(DebugLevel()>1) cout << "Converting clusters to hits" << endl; ConvertClustersToHits(); // also low max ADC clusters are removed if(DebugLevel()>1) cout << "Outputting hit table on event node" << endl; outNode->AddDataTable(fHitTable); } //_______________________________________________________ Bool_t BrTpcHitModule::FindHitWeightedMean(BrTpcCluster *cluster) { // This is the old FindHits2 method used to calculate the hit // position of the hits using the weighted mean method const Int_t nSeqs = cluster->GetNseqs(); const Int_t currentRow = cluster->GetRow(); Double_t *pcharge = new Double_t[nSeqs]; // integrated pad charge Double_t *dpcharge = new Double_t[nSeqs]; // estimated error Double_t *xpad = new Double_t[nSeqs]; // x positions of pads Int_t minTime = 1000; Int_t maxTime = 0; Int_t nHitsPadSummation = 0; Double_t totalChargePadSummation = 0.0; // Loop over sequences Int_t j; for (j = 0; j < nSeqs; j++) { BrTpcSequence *jseq = cluster->GetSeq(j); const Int_t pad = jseq->GetPad(); const Int_t nBuckets = jseq->GetNseq(); const Int_t timeBin = jseq->GetTime(); Float_t adcGain=1.0; if(fUsePadGain) adcGain = fCalibration->GetPadAdcGain(currentRow, pad); if(timeBin < minTime) minTime = timeBin; if(timeBin+nBuckets-1 > maxTime) maxTime = timeBin+nBuckets-1; xpad[j] = pad; pcharge[j] = 0.0; dpcharge[j] = 0.0; // loop over time bins in sequences for (Int_t k = 0; k < nBuckets; k++) { pcharge[j] += adcGain * (jseq->GetAdc(k)); // integrating all buckets for this pad nHitsPadSummation++; } dpcharge[j] = TMath::Sqrt(Double_t(nBuckets)) * fWhiteNoise; totalChargePadSummation += pcharge[j]; } // end loop over sequences // Fit the pad position Double_t xpos = 0.0; // fitted x position Double_t dxpos = 0.0; // fitted error Double_t xwidth = 0.0; // caculated PRF width WeightedMean (nSeqs, xpad, pcharge, dpcharge, xpos, xwidth, dxpos ); // Proceed with the time distribtuion const Int_t nTimeBuckets = maxTime - minTime + 1; if(DebugLevel()>3) cout << "[min;max] time for cluster : [" << minTime << ";" << maxTime << "] => nTimeBuckets = " << nTimeBuckets << endl; Double_t *tcharge = new Double_t[nTimeBuckets]; // integrated pad charge Double_t *dtcharge = new Double_t[nTimeBuckets]; // estimated error Double_t *ytime = new Double_t[nTimeBuckets]; // positions of timebuckets for (Int_t k = 0; k < nTimeBuckets; k++){ ytime[k] = minTime + k; tcharge[k] = 0; dtcharge[k] = 0; } Int_t nHitsTimeSummation = 0; Double_t totalChargeTimeSummation = 0; for (j = 0; j < nSeqs; j++ ) { BrTpcSequence *jseq = cluster->GetSeq(j); const Int_t nBuckets = jseq->GetNseq(); const Int_t timeBin = jseq->GetTime(); const Int_t pad = jseq->GetPad(); Int_t k = timeBin-minTime; Float_t adcGain=1.0; if(fUsePadGain) adcGain = fCalibration->GetPadAdcGain(currentRow, pad); for (Int_t n = 0; n < nBuckets; n++) { tcharge[k] += adcGain* (Float_t) jseq->GetAdc(n); // integrating all buckets for this pad dtcharge[k] += fWhiteNoise*fWhiteNoise; totalChargeTimeSummation += adcGain* jseq->GetAdc(n); k++; nHitsTimeSummation++; } // end loop over sequence timebins } // end loop over sequences for (Int_t k=0; k<nTimeBuckets; k++) dtcharge[k] = TMath::Sqrt(dtcharge[k]); Assert(nHitsPadSummation == nHitsTimeSummation); Assert(totalChargePadSummation == totalChargeTimeSummation); // Fit the time position Double_t ypos = 0.0; // fitted y position Double_t dypos = 0.0; // fitted error Double_t ywidth = 0.0; // caculated width WeightedMean(nTimeBuckets, ytime, tcharge, dtcharge, ypos, ywidth, dypos); if(HistOn()) { hPadTime->Fill(xpos, ypos); hPadRow->Fill(xpos, cluster->GetRow()); } // Add tpc hit to tpc hit table BrTpcHit *hit = new BrTpcHit(); // convert positions and scales to local positions xpos = fParams->PadToIntrinsic(xpos); dypos = fParams->TimeToIntrinsic(cluster->GetRow(), dypos+ypos); ypos = fParams->TimeToIntrinsic(cluster->GetRow(), ypos); dxpos = dxpos*fParams->GetPadDistance(); dypos = TMath::Abs(ypos - dypos); Double_t zpos = fParams->GetRowPosition(cluster->GetRow()); // set error on positions if asked for if(fUseCalibratedWidths) { dxpos = fParams->GetHitWidthX(); dypos = fParams->GetHitWidthY(); if(dxpos <= 0 || dypos <= 0) Abort("FindHitWeightedMean", "dxpos or dypos negative"); } if(HistOn()) { hX->Fill(xpos); hY->Fill(ypos); hDx->Fill(dxpos); hDy->Fill(dypos); } // set positions hit->SetX(xpos); hit->SetY(ypos); hit->SetZ(zpos); hit->SetXError(dxpos); hit->SetYError(dypos); // maybe z error could be set to 1/12 of the pad length // Now set the TPC specific details hit->SetPadWidth(xwidth); hit->SetTimeWidth(ywidth); hit->SetClusterId(cluster->GetID()); hit->SetRow(cluster->GetRow()); hit->SetStatus(cluster->GetStatus()); hit->SetAdcMax(cluster->GetMaxADC()); hit->SetAdcSum(cluster->GetADCSum()); fHitTable->Add(hit); // Clean up and exit delete [] tcharge; delete [] dtcharge; delete [] ytime; delete [] pcharge; delete [] dpcharge; delete [] xpad; return kTRUE; } //____________________________________________________________________ void BrTpcHitModule::Print(Option_t* option) const { // Module Information method // // Options: (see also BrModule::Print) // BrModule::Print(option); TString opt(option); opt.ToLower(); if (opt.Contains("d")) cout << endl << " Original author: Peter H.L. Christiansen" << endl << " Last Modifications: " << endl << " $Author: videbaek $" << endl << " $Date: 2002/08/15 17:25:09 $" << endl << " $Revision: 1.7 $ " << endl << endl << "-------------------------------------------------" << endl; cout << " WhiteNoise = " << fWhiteNoise << endl; } //_______________________________________________________ void BrTpcHitModule::WeightedMean(const Int_t npts, const Double_t *x, const Double_t *y, const Double_t *dy, Double_t ¢roid, Double_t &width, Double_t &uncertainty) { // This is the routine that calculates the best estimate for the // weighted mean and the error and the width of the distribtuion Double_t mom1 = 0.0; Double_t mom2 = 0.0; Double_t ysum = 0.0; for( Int_t i=0; i<npts; i++ ) { ysum += y[i]; mom1 += x[i] * y[i]; mom2 += x[i] * x[i] * y[i]; } mom1 = mom1 / ysum; mom2 = mom2 / ysum - ( mom1 * mom1 ); if ( mom2 > 0.0 ) mom2 = TMath::Sqrt( mom2 ); else mom2 = 0; // Assuming the noise is uncorrelated uncertainty = 0.0; Int_t i; for (i = 0; i < npts; i++) { Double_t temp = (x[i] - mom1) * dy[i]; uncertainty += temp * temp; } centroid = mom1; width = mom2; uncertainty = TMath::Sqrt( uncertainty ) / ysum ; if (DebugLevel()>3){ cout << " Centroid: " << centroid << " Width: " << width << " Unc: " << uncertainty << endl; } return; } //____________________________________________________________________ // // $Log: BrTpcHitModule.cxx,v $ // Revision 1.7 2002/08/15 17:25:09 videbaek // Add code to correct for adcgain if in database. // he default is not to use this. // // Revision 1.6 2002/01/03 19:53:26 cholm // Prepared to use BrTableNames class (or perhaps BrDetectorList) for table names // // Revision 1.5 2001/10/12 10:56:00 pchristi // Added option to use the calibrated hit widths from the database. // You have to do SetUseCalibratedWidths(kTRUE) to and have them // read in. // // Revision 1.4 2001/08/14 19:29:33 pchristi // Changed the Print methods to pass the option to BrModule. // Also removed debug statement in BrTpcTrackFollowModule. // // Revision 1.3 2001/07/06 10:24:18 pchristi // Changed assert to the ROOT Assert in TError.h // // Revision 1.2 2001/06/22 17:45:29 cholm // Change names so that every data class has the same format, so that // we will not have to worry about that later on. The affected classes // are: // // BrTPCSequence -> BrTpcSequence // BrTPCCluster -> BrTpcCluster // BrTPCClusterTable -> BrTpcClusterTable // // These changes has ofcourse been propegated to the modules as well, // giving the changes // // BrTPCClusterFinder -> BrTpcClusterFinder // BrTPCSequenceAdder -> BrTpcSequenceAdder // // Revision 1.1.1.1 2001/06/21 14:55:13 hagel // Initial revision of brat2 // // Revision 1.2 2001/06/17 17:29:45 pchristi // Near final version. A little cosmetics is still needed. // // Revision 1.1 2001/05/29 14:20:39 pchristi // Initial import of new Tpc classes // // |
||||||
|
This page automatically generated by script docBrat by Christian Holm |
Copyright ; 2002 BRAHMS Collaboration
<brahmlib@rcf.rhic.bnl.gov>
|