|
//____________________________________________________________________ // // The BrTpcClusterModule is used to group the sequences into // clusters, and store the clusters in a cluster table. // This Class has been taken from what M. Justice started to put // together in Oct 99 // // It takes sequences that have been filtered with // BrTPCPreProcess as input and groups the sequences into // BrTpcCluster s // // Deconvolution has been moved to BrTpcDeconvoluteClusterModule // //____________________________________________________________________ // // $Id: BrTpcClusterModule.cxx,v 1.8 2002/06/13 18:46:33 videbaek Exp $ // $Author: videbaek $ // $Date: 2002/06/13 18:46:33 $ // $Copyright: (C) 2001 BRAHMS Collaboration <brahmlib@rhic.bnl.gov> // #include <BrTpcClusterModule.h> #ifndef WIN32 #include <iostream> #else #include <iostream.h> #endif #include <iomanip.h> #ifndef ROOT_TDirectory #include "TDirectory.h" #endif #ifndef BRAT_BrTableNames #include "BrTableNames.h" #endif //__________________________________________________________________ ClassImp(BrTpcClusterModule); //__________________________________________________________________ BrTpcClusterModule::BrTpcClusterModule() : BrModule() { // This constructer is no good SetState(kSetup); fClusterTable = 0; fSeqTable = 0; fReset = kFALSE; } //__________________________________________________________________ BrTpcClusterModule::BrTpcClusterModule(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; fSeqTable = 0; fReset = kFALSE; } //_________________________________________________________________ BrTpcClusterModule::~BrTpcClusterModule() { // empty } //_________________________________________________________________ void BrTpcClusterModule::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,"Cluster_%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}" ); // Restore directory savdir->cd(); } //_______________________________________________________________________ void BrTpcClusterModule::Init() { // Initialize module. Empty. SetState(kInit); } //_________________________________________________________________ void BrTpcClusterModule::Event(BrEventNode *inputNode, BrEventNode *outputNode) { // Main method // Gets the TPCSequence data from the inputNode // Group the sequences into clusters. // Output the clusters on outputNode SetState(kEvent); if(DebugLevel()>2) cout << "Entering " << GetName() << "BrTpcClusterModule Event()" << endl; fSeqTable = inputNode->GetDataTable(Form("%s %s", BRTABLENAMES kTPCSequence, GetName())); if(fSeqTable == 0) { if(DebugLevel() > 0 ) Warning("Event", "BrDataTable %s %s not found in this event!", BRTABLENAMES kTPCSequence, GetName()); return; } // Sort the datatable using the method // BrTpcSequence::Compare(const TObject*) const // This makes it easier to group the data fSeqTable->Sort(); // Create the output node fClusterTable = new BrDataTable(Form("%s %s", BRTABLENAMES kTPCCluster, GetName())); FindClusters(); // Append output data to output event node outputNode->AddDataTable(fClusterTable); if ( fReset ) ResetSequenceTable(); } //__________________________________________________________________ void BrTpcClusterModule::ResetSequenceTable() { // Reset BrTpcSequence::fClustNum to -1, as this variable is modified in FindClusters. if ( !fSeqTable ) return; Int_t nSeq = fSeqTable->GetEntries(); if ( !nSeq ) return; for ( Int_t iS = 0; iS < nSeq; iS++ ) ((BrTpcSequence*) fSeqTable->At( iS ))->SetClustNum( -1 ); } //__________________________________________________________________ Bool_t BrTpcClusterModule::FindClusters() { // ************************************************************** // * * // * The algorithm used here is to group together sequences * // * on adjacent pads which overlap in the temporal direction * // * as determined by the SeqOverlap function * // * * // ************************************************************** BrTpcCluster *currCluster = 0; fClusterId = 0; const Int_t nSeq = fSeqTable->Entries(); Debug(3,"Find Clusters","Entry - Total sequences %d",nSeq); for (Int_t i = 0; i < nSeq; i++) { BrTpcSequence *iseq = (BrTpcSequence*) fSeqTable->At(i); Int_t irow = iseq->GetRow(); Int_t ipad = iseq->GetPad(); Int_t iclustID = iseq->GetClustNum(); // If the cluster id is < 0 it does not belong to a cluster if (iclustID < 0) { iclustID = fClusterId++; iseq->SetClustNum( iclustID ); // make a new cluster currCluster = new BrTpcCluster(irow, *iseq); fClusterTable->Add(currCluster); currCluster->SetID(iclustID); if(DebugLevel()>5) { cout << "New Cluster for Sequence " << i << " / " << nSeq << "n Row " << irow << " - Pad " << ipad << "n ID " << iclustID << endl; } } else { currCluster = GetClusterById(iclustID); Assert( currCluster ); if(DebugLevel()>3) { cout << "Sequence " << i << " / " << nSeq << " : old cluster ID " << iclustID << endl; currCluster->Print(); } } for (Int_t j = (i+1); j < nSeq; j++) { BrTpcSequence *jseq = (BrTpcSequence*) fSeqTable->At(j); Int_t jrow = jseq->GetRow(); // if we are in the next row stop looking // (we have sorted the table) if( jrow != irow ) break; Int_t jpad = jseq->GetPad(); if( jpad != (ipad+1) ) continue; if(DebugLevel()>5) cout << " Check seq("<<j<<") : " << jseq->GetClustNum() << ", pad = " << jseq->GetPad() << " overlap" << endl; if ( SeqOverlap(iseq, jseq) ) { Int_t jclustID = jseq->GetClustNum(); if ( DebugLevel()>3 ) cout << "Cluster ID for overlapping sequence : " << jclustID << endl; // If the sequence does not belong to a cluster if (jclustID == -1) { // add sequence jseq->SetClustNum( iclustID ); currCluster->AddSeq( *jseq ); if(DebugLevel()>3) cout << " Add Sequence " << jseq->GetClustNum()<<" in pad " <<jpad << "to cluster" << endl; // // If the sequence already belongs to a cluster (but not icluster) // add whole jcluster to icluster and delete jcluster } else if ( jclustID != iclustID ) { if(DebugLevel()>3){ cout << " Add old cluster : " << setw(4) << jclustID << " to cluster : " << setw(4) << iclustID; } BrTpcCluster *clustOld = GetClusterById(jclustID); // test that cluster old exists Assert( clustOld ); TObjLink *lnk = clustOld->GetSeqList()->FirstLink(); while ( lnk ) { BrTpcSequence *seqOld = (BrTpcSequence*)lnk->GetObject(); // have to set same ID for the sequence in the BrDataTable BrTpcSequence *tabSeq = (BrTpcSequence*)fSeqTable->GetObjectList() ->FindObject( seqOld ); // test that we do get the right sequence Assert( tabSeq ); Assert( tabSeq->IsEqual(seqOld) ); tabSeq->SetClustNum( iclustID ); currCluster->AddSeq( *tabSeq ); // make a copy lnk = lnk->Next(); } // Finally we remove the old cluster fClusterTable->DeleteAndCompress(clustOld); Assert(GetClusterById(jclustID) == 0 ); } } } } if(DebugLevel()>0) cout << "Number of clusters in TPC " << GetName() << " = " << fClusterTable->GetEntries() << endl; // Fill histograms if(HistOn()) { const Int_t nEnt = fClusterTable->GetEntries(); hNclusters->Fill(nEnt); for(Int_t n = 0; n < nEnt; n++) { BrTpcCluster *cluster = (BrTpcCluster *)fClusterTable->At(n); hNpads->Fill(cluster->GetNpads()); hNtime->Fill(cluster->GetNbins()); hMaxADC->Fill(cluster->GetMaxADC()); hSumADC->Fill(cluster->GetADCSum()); hPadvar->Fill(cluster->GetVarPad()); hTimevar->Fill(cluster->GetVarTime()); } } return kTRUE; } //__________________________________________________________________________ BrTpcCluster* BrTpcClusterModule::GetClusterById(Int_t id) { // Loop over fClusterTable and return cluster with matching id // If cluster is not found return 0 const Int_t nEntries = fClusterTable->GetEntries(); // do the loop counting down since it is most likely to be close to // the end = this row for(Int_t i = nEntries-1; i >=0; i--) { BrTpcCluster *cluster = (BrTpcCluster *) fClusterTable->At(i); if(cluster->GetID() == id) return cluster; } return 0; } //____________________________________________________________________ void BrTpcClusterModule::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/06/13 18:46:33 $" << endl << " $Revision: 1.8 $ " << endl << endl << "-------------------------------------------------" << endl; } //__________________________________________________________________________ Bool_t BrTpcClusterModule::SeqOverlap(BrTpcSequence *iseq, BrTpcSequence *jseq) { // *************************************************************** // * * // * Overlap criteria: * // * * // * | | * // * | | * // * TRUE | FALSE | * // * | | | * // * | | | | * // * | | * // * | | * // * | * // * | * // * * // * This could be made more flexible by introduction * // * of a parameter controlling the amount of overlap * // * required. * // * * // *************************************************************** // // Maybe this method should be in the BrTpcSequence Int_t ifirst = iseq->GetTime(); Int_t ilast = ifirst + iseq->GetNseq()-1; Int_t jfirst = jseq->GetTime(); Int_t jlast = jfirst + jseq->GetNseq()-1; if ( ( (jfirst > ifirst) && (jfirst < ilast) ) || ( (jlast > ifirst) && (jlast < ilast) ) || ( (ifirst > jfirst) && (ifirst < jlast) ) || ( (ilast > jfirst) && (ilast < jlast) ) || ( (ifirst == jfirst) && (ilast == jlast) ) ){ return kTRUE; } else { return kFALSE; } } //////////////////////////////////////////////////////////////////////////// // // $Log: BrTpcClusterModule.cxx,v $ // Revision 1.8 2002/06/13 18:46:33 videbaek // add debug information // // Revision 1.7 2002/01/03 19:53:11 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/06 10:23:48 pchristi // Changed assert to the ROOT Assert in TError.h // // Revision 1.3 2001/06/25 14:40:58 cholm // Use BrTpcCluster::Print rather than BrTpcCluster::List // // Revision 1.2 2001/06/22 17:45:14 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:33 pchristi // Near final version. A little cosmetics is still needed. // // Revision 1.1 2001/05/29 14:20:17 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>
|