|
//____________________________________________________________________ // // The BrTpcDeconvoluteClusterModule used to be part of // BrTpcClusterFinder // // It takes clusters sequences that come from BrTpcClusterModule and // examines the status. If removes noise clusters and tries to // deconvolute multiclusters. // The output data is (also) clusters in the input datatable // //____________________________________________________________________ // // $Id: BrTpcDeconvoluteClusterModule.cxx,v 1.8 2002/08/15 17:23:31 videbaek Exp $ // $Author: videbaek $ // $Date: 2002/08/15 17:23:31 $ // $Copyright: (C) 2001 BRAHMS Collaboration <brahmlib@rhic.bnl.gov> #ifndef BRAT_BrTpcDeconvoluteClusterModule #include "BrTpcDeconvoluteClusterModule.h" #endif #ifndef ROOT_TObjArray #include "TObjArray.h" #endif #ifndef ROOT_TH1 #include "TH1.h" #endif #ifndef ROOT_TDirectory #include "TDirectory.h" #endif #ifndef BRAT_BrParameterDbManager #include "BrParameterDbManager.h" #endif #ifndef BRAT_BrDetectorParamsTPC #include "BrDetectorParamsTPC.h" #endif #ifndef BRAT_BrTableNames #include "BrTableNames.h" #endif #ifndef WIN32 #include <iostream> #else #include <iostream.h> #endif #include <iomanip.h> //__________________________________________________________________ ClassImp(BrTpcDeconvoluteClusterModule); //__________________________________________________________________ BrTpcDeconvoluteClusterModule::BrTpcDeconvoluteClusterModule() : BrModule() { // This constructer is no good SetState(kSetup); fClusterTable = 0; fAdcTable = 0; } //__________________________________________________________________ BrTpcDeconvoluteClusterModule::BrTpcDeconvoluteClusterModule(const Char_t *name, const Char_t *title) : BrModule(name, title) { // This is the constructor that should be used // Name should be "T1", "T2", "TPM1", "TPM2" SetState(kSetup); fClusterTable = 0; fAdcTable = 0; // define default cut parameters SetHighPsigCut(); SetHighTsigCut(); SetMinPadCut(); SetMinTimeCut(); SetMinMaxADCCut(); } //_________________________________________________________________ BrTpcDeconvoluteClusterModule::~BrTpcDeconvoluteClusterModule() { if( fClusterTable ) { fClusterTable->Clear(); delete fClusterTable; } if(fAdcTable) delete fAdcTable; } //_______________________________________________________________________ Bool_t BrTpcDeconvoluteClusterModule::DeconvoluteClusters() { /************************************************************** * * * The purpose of this function is to break up any cluster * * which has the topology indicating that it might contain * * information from more than 1 track. Such clusters would * * have the feature of more than 1 local maxima. (Of course * * 2 tracks could be so close that only 1 maximum results.) * * * * This is a difficult routine to write but it is very * * important, especially for HBT pions. For now, I simply * * fill in a skeleton of how one might approach the problem. * * * * M. Justice * * 16-Nov-1999 * * It is worth noting it is important for any even modest * * track densities in the TPC as learned in other methods(fv)* * * **************************************************************/ Int_t nSubClusters = 0; Int_t nNoiseClusters = 0; Int_t nDeleted = 0; TObjArray subClusters(2); // autoexpands if more than 2 subclusters const Int_t nEntries = fClusterTable->GetEntries(); Debug(1,"DeconvoluteClusters"," Table has %d entries", nEntries); for(Int_t i = 0; i < nEntries; i++) { subClusters.Clear(); BrTpcCluster *cluster = (BrTpcCluster *)fClusterTable->At(i); // The next method checks the status of the cluster and // deconvolutes if possible It returns number of peaks // (=subclusters) and the subclusters in the TObjArray Debug(3,"FindSubClusters"," cluster # %d ", i); const Int_t nPeaks = FindSubClusters(cluster, &subClusters); const Int_t status = cluster->GetStatus(); if(status == BrTpcCluster::kNoiseHit) { Debug(5," ","is noisehit"); // remove and delete the cluster fClusterTable->DeleteObjectAt(i); nDeleted++; nNoiseClusters++; } else if (status == BrTpcCluster::kSingleHit) { Debug(5," ","is singlehit"); // do nothing } else if (status == BrTpcCluster::kMultiHit) { Debug(5," ","is multihit %d peaks",nPeaks); // remove cluster after we added good deconv. clusters Int_t nP = 0; Int_t j = 0; for (j = 0; j < nPeaks; j++) { BrTpcCluster *jcluster = (BrTpcCluster *)subClusters.At(j); // Add jcluster only if not noise if (jcluster->GetStatus() != BrTpcCluster::kNoiseHit){ // Add jcluster after the last in our loop // wait with compression // fClusterTable->AddAt(jcluster, nEntries+nSubClusters); fClusterTable->Add(jcluster); nDeleted--; nSubClusters++; nP++; } else delete jcluster; } // remove and delete the old cluster fClusterTable->DeleteObjectAt(i); nDeleted++; if(HistOn()) hNpeaks->Fill(nP); } else { // How did we ever get here ? Warning("DeconvoluteClusters", "Illegal cluster status"); } } // Important to remember to compress the table fClusterTable->Compress(); if(DebugLevel()>0) cout << "Number of final clusters = " << fClusterTable->GetEntries() << endl; if(DebugLevel()>1) { cout << "Number of noise clusters = " << nNoiseClusters << endl << "Number of sub clusters = " << nSubClusters << endl; } const Int_t n = fClusterTable->GetEntries(); if(n+nDeleted != nEntries) { Warning("deconvolute","values %f %f %f ", nEntries,nDeleted,n); } // Fill histograms if applicable if(HistOn()) { hNnoiseclusters->Fill(nNoiseClusters); hNsubclusters->Fill(nSubClusters); const Int_t n = fClusterTable->GetEntries(); hNclusters->Fill(n); Int_t i = 0; for(i = 0; i < n; i++) { BrTpcCluster *cluster = (BrTpcCluster *)fClusterTable->At(i); hNpads->Fill(cluster->GetNpads()); hNtime->Fill(cluster->GetNbins()); hMaxADC->Fill(cluster->GetMaxADC()); hSumADC->Fill(cluster->GetADCSum()); hPadvar->Fill(cluster->GetVarPad()); hTimevar->Fill(cluster->GetVarTime()); } } int nEntriesAtEnd = fClusterTable->GetEntries(); if(DebugLevel()>1) cout << "Cluster Table has " << nEntriesAtEnd << endl; return kTRUE; } //_______________________________________________________________________ void BrTpcDeconvoluteClusterModule::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* savdir = gDirectory; Char_t Dirname[32]; sprintf(Dirname,"DeconvoluteCluster_%s",GetName()); TDirectory* hdir = savdir->mkdir(Dirname); if (!hdir) { Warning("DefineHistograms","Could not create histogram subdirectory"); return; } hdir->cd(); Char_t Histname[32]; Char_t HistTitle[48]; sprintf(Histname,"hNpads_%s",GetName()); sprintf(HistTitle,"%s: # pads per cluster",GetName()); hNpads = new TH1F(Histname, HistTitle,20,0.5,20.5); hNpads->SetXTitle( "N_{pads}" ); hNpads->SetYTitle( "dN/dN_{pads}" ); sprintf(Histname,"hNtimebins_%s",GetName()); sprintf(HistTitle,"%s: # time bins per cluster",GetName()); hNtime = new TH1F(Histname, HistTitle,20,0.5,20.5); hNtime->SetXTitle( "N_{time bins}" ); hNtime->SetYTitle( "dN/dN_{time bins}" ); sprintf(Histname,"hNclusters_%s",GetName()); sprintf(HistTitle,"%s: # clusters per event",GetName()); hNclusters = new TH1F(Histname, HistTitle,800,-0.5,799.5); hNclusters->SetXTitle( "N_{clusters}" ); hNclusters->SetYTitle( "dN/dN_{clusters}" ); sprintf(Histname,"hMaxADC_%s",GetName()); sprintf(HistTitle,"%s: max ADC in clusters",GetName()); hMaxADC = new TH1F(Histname, HistTitle,1024,0.5,1024.5); hMaxADC->SetXTitle( "ADC_{max}" ); hMaxADC->SetYTitle( "dN/dADC_{max}" ); sprintf(Histname, "hSumADC_%s", GetName()); sprintf(HistTitle, "%s: integrated ADC in clusters", GetName()); hSumADC = new TH1F(Histname, HistTitle, 1024, 0, 10240); hSumADC->SetXTitle( "ADC_{sum}" ); hSumADC->SetYTitle( "dN/dADC_{sum}" ); sprintf(Histname,"hPadvar_%s",GetName()); sprintf(HistTitle,"%s: pad sigma",GetName()); hPadvar = new TH1F(Histname, HistTitle,40,0.0,10.0); hPadvar->SetXTitle( "#sigma_{pad} [pads]" ); hPadvar->SetYTitle( "dN/d#sigma_{pad}" ); sprintf(Histname,"hTimevar_%s",GetName()); sprintf(HistTitle,"%s: time sigma",GetName()); hTimevar = new TH1F(Histname, HistTitle,40,0.0,10.0); hTimevar->SetXTitle( "#sigma_{time bin} [time bin]" ); hTimevar->SetYTitle( "dN/d#sigma_{time bin}" ); sprintf(Histname,"hNnoiseclusters_%s",GetName()); sprintf(HistTitle,"%s: # noise clusters per event",GetName()); hNnoiseclusters = new TH1F(Histname, HistTitle,800,-0.5,799.5); hNnoiseclusters->SetXTitle( "N_{clusters}" ); hNnoiseclusters->SetYTitle( "dN/dN_{clusters}" ); sprintf(Histname,"hNsubclusters_%s",GetName()); sprintf(HistTitle,"%s: # good sub clusters per event",GetName()); hNsubclusters = new TH1F(Histname, HistTitle,800,-0.5,799.5); hNsubclusters->SetXTitle( "N_{subclusters}" ); hNsubclusters->SetYTitle( "dN/dN_{subclusters}" ); sprintf(Histname,"hNpeaks_%s",GetName()); sprintf(HistTitle,"%s: # good sub clusters per cluster",GetName()); hNpeaks = new TH1F(Histname, HistTitle,100,-0.5,99.5); hNpeaks->SetXTitle( "N_{subclusters}" ); hNpeaks->SetYTitle( "dN/dN_{subclusters}" ); // Restore directory savdir->cd(); } //_________________________________________________________________ void BrTpcDeconvoluteClusterModule::Event(BrEventNode *inputNode, BrEventNode *outputNode) { // Main method // Gets the TPCCluster data from the inputNode // Examines the status of the clusters and tries to deconvolute // multihits // The output is the same table as the input so no new table is // appended, i.e, the changes are made to the old table. SetState(kEvent); if(DebugLevel()>0) cout << Form("Entering BrTpcDeconvoluteClusterModule (%s) Event()", GetName()) << endl; fClusterTable = inputNode->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; } // Check that there are any clusters in the table or the id trick crashes const Int_t nEntries = fClusterTable->GetEntries(); if(nEntries==0) return; // Set the cluster id to be the id of the last of the clusters in // the table + 1 fClusterId = ((BrTpcCluster*)fClusterTable-> At((fClusterTable->GetEntries()-1)))->GetID() + 1; if(DebugLevel()>2) cout << "Deconvoluting clusters" << endl; DeconvoluteClusters(); // also low max ADC clusters are removed // The data table is still in the event node so no need to add // anything } //_______________________________________________________________________ void BrTpcDeconvoluteClusterModule::FillSubClusters(Int_t nPeaks, TObjArray *subClusters) { // // Scan the working array and create a set of subclusters // according to the fClusterIdTable identification. // The appropriate set of BrTPSequence objects that make up the // clusters are created. These are of course NOT the inputdata // but a subset. // This routine needs access to fCurrentRow, and the next available ClusterId // to be used. for(Int_t peak = 0; peak < nPeaks; peak++) { BrTpcCluster *cluster = new BrTpcCluster(); cluster->SetRow(fCurrentRow); cluster->SetID(fClusterId++); subClusters->Add(cluster); } // // scan each pad and find sequence in this pad. // const Int_t minPad = fAdcTable->GetMinPad(); const Int_t maxPad = fAdcTable->GetMaxPad() + 1; const Int_t minTime = fAdcTable->GetMinTime(); const Int_t maxTime = fAdcTable->GetMaxTime() + 1; for (Int_t pad = minPad; pad < maxPad; pad++) { Int_t timebin = minTime; while(timebin < maxTime){ // Check that pixel was used const Int_t id = fAdcTable->GetId(pad, timebin); if (id < 1 ){ timebin++; } else { const Int_t firstTimebin = timebin; // Find the length of the sequence Int_t nSeq = 0; while(timebin < maxTime && fAdcTable->GetId(pad, timebin) == id) { nSeq++; timebin++; } BrTpcSequence *seq = new BrTpcSequence(nSeq); seq->SetPad(pad); seq->SetRow(fCurrentRow); seq->SetTime(firstTimebin); for(Int_t n = 0; n < nSeq; n++) seq->SetAdc(n, fAdcTable->GetAdc(pad, firstTimebin+n)); if(DebugLevel()>3) cout << "Add seq " <<fCurrentRow<<" " << pad << " " << firstTimebin << " " << nSeq <<endl; BrTpcCluster *cluster = (BrTpcCluster *)subClusters->At(id-1); seq->SetClustNum(cluster->GetID()); // Add the sequence to the cluster (and the ownership) cluster->AddSeq(*seq, kFALSE); } } } } //____________________________________________________________ void BrTpcDeconvoluteClusterModule::FillWorkingArrays(BrTpcCluster *cluster) { // // Fill the adcTable with adc values in prepartion for // subclusterfinding. // fCurrentRow = cluster->GetRow(); // Used when table --> clusters const Int_t nSeqs = cluster->GetNseqs(); for (Int_t j = 0; j < nSeqs; j++) { BrTpcSequence *jseq = cluster->GetSeq(j); const Int_t pad = jseq->GetPad(); const Int_t nbkts = jseq->GetNseq(); const Int_t timebin = jseq->GetTime(); for (Int_t k = 0; k < nbkts; k++) { fAdcTable->SetAdc(pad, timebin + k, jseq->GetAdc(k)); fAdcTable->SetId(pad, timebin + k, 0); } } } //___________________________________________________________ Int_t BrTpcDeconvoluteClusterModule::FindHitStatus(BrTpcCluster* cluster) { // Determine the status of each cluster // The outcome can be either // noisehit (throw away) // single hit // multi hit (deconvolution needed) // Apply noise cuts if (cluster->GetNpads() <= fMinPadCut || cluster->GetNbins() <= fMinTimeCut || cluster->GetMaxADC()<= fMinMaxADCCut) { return (BrTpcCluster::kNoiseHit); } // Check that both widths are small if( (cluster->GetVarPad() < fHighPsigCut) && (cluster->GetVarTime() < fHighTsigCut) ) return (BrTpcCluster::kSingleHit); else return (BrTpcCluster::kMultiHit); } //_________________________________________________________________________ Int_t BrTpcDeconvoluteClusterModule::FindSubClusters(BrTpcCluster *primaryCluster, TObjArray *subClusters) { // primaryCluster is scanned with FindHitStatus to determine if its // noise, a good single hit cluster or a multihit candidate. The // TObjArray is empty initially. Ownership of new BrTpcClusters from // deconvolution is transferred to this array on exit. // // Deconvolution is done as follows : // - Fill working array // - Finds peaks and work towards valleys // - Identify part sequences belonging to same peak // - Create subClusters and new BrTpcSequences for these clusters // - // Int_t status = FindHitStatus(primaryCluster); if(status == BrTpcCluster::kSingleHit || status == BrTpcCluster::kNoiseHit) { primaryCluster->SetStatus(status); return 1; } primaryCluster->SetStatus(BrTpcCluster::kMultiHit); // Clear the working 2d array fAdcTable->Clear(); // Fill the ADC values into the working array FillWorkingArrays(primaryCluster); // Find Peaks const Int_t nPeaks = fAdcTable->FindPeaks(Short_t(fMinMaxADCCut)); // nPeaks should always be greater than 0 if (nPeaks < 1) { fAdcTable->Print(); Stop("DeconvoluteClusters", "Illegal number of peaks"); } // Fill the sub clusters corresponding to each peak FillSubClusters(nPeaks, subClusters); // Find and set status of subclusters for(Int_t peak = 0;peak < nPeaks; peak++){ BrTpcCluster *cluster = (BrTpcCluster*) subClusters->At(peak); status = FindHitStatus(cluster); if(status == BrTpcCluster::kNoiseHit) cluster->SetStatus(status); if(status == BrTpcCluster::kSingleHit) cluster->SetStatus(BrTpcCluster::kSingleHitDec); if(status == BrTpcCluster::kMultiHit) cluster->SetStatus(BrTpcCluster::kMultiHitDec); } return nPeaks; } //_______________________________________________________________________ void BrTpcDeconvoluteClusterModule::Init() { // Initialize module. This is the proper member fct for getting // information from the data base. SetState(kInit); BrParameterDbManager* gParamDb = BrParameterDbManager::Instance(); BrDetectorParamsTPC *params = (BrDetectorParamsTPC*) gParamDb->GetDetectorParameters("BrDetectorParamsTPC", GetName()); // set parameters const Int_t padsPerRow = params->GetPadsprow(); const Int_t noBuckets = params->GetNoBuckets(); fAdcTable = new BrTpcAdcTable(padsPerRow, noBuckets); } //____________________________________________________________________ void BrTpcDeconvoluteClusterModule::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:23:31 $" << endl << " $Revision: 1.8 $ " << endl << endl << "-------------------------------------------------" << endl; cout << " highPsigCut = " << fHighPsigCut << endl << " highTsigCut = " << fHighTsigCut << endl << " minPadCut = " << fMinPadCut << endl << " minTimeCut = " << fMinTimeCut << endl << " minMaxADCCut = " << fMinMaxADCCut << endl; } //////////////////////////////////////////////////////////////////////////// // // $Log: BrTpcDeconvoluteClusterModule.cxx,v $ // Revision 1.8 2002/08/15 17:23:31 videbaek // Add better diagnostics code. // // Revision 1.7 2002/01/03 19:53:19 cholm // Prepared to use BrTableNames class (or perhaps BrDetectorList) for table names // // Revision 1.6 2001/08/16 15:45:17 jens // // BrTpcClusterModule // BrTpcDeconvoluteClusterModule // BrTpcTrackModule // Changed some variable names for histograms and the title scheme for histogram. // Now histograms have title like e.g. "<TPC name>: max ADC in clusters". // // BrTpcClusterModule // Added the possibility to reset sequences in fSeqTable. BrTpcSequence::fClustNum // is altered in BrTpcClusterModule::FindClusters, and thus made it impossible to // reanalyze an event. I.e. previously a new instance of BrEvent had to be // allocated and filled from input. // // Revision 1.5 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.4 2001/07/17 20:30:50 ejkim // Fixed bug when number of clusters in input was 0. // // Revision 1.3 2001/07/06 10:24:05 pchristi // Changed assert to the ROOT Assert in TError.h // // Revision 1.2 2001/06/22 17:45:19 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:36 pchristi // Near final version. A little cosmetics is still needed. // // Revision 1.1 2001/05/29 14:20:26 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>
|