|
// $Id: BrDigitizeTPC.cxx,v 1.5 2001/11/13 21:11:24 trulsml Exp $ // $Author: trulsml $ // $Date: 2001/11/13 21:11:24 $ // //________________________________________________________________________ // // The BRDigitize TPC is the digitization module for the // pad TPC's. The module performs a slow digitization using the // geant hit information. The output is in tables of the BrTPCsequence // i.e. a structure that very closely resembles the output from the // TPC front ends. // // The Detector geometry comes from two sources // BrDetectorVolume which is set using the BrGeometryDbManager // and from the BrDetectorParamsTPC // // The method has the standard Module entry points implemented // // Init() // Event() // DefineHistograms() (called from Book() // // Note Init() must be called explicitely (Nov 2000). It has been removed from the // constructor. // // December 1999 // The conversion between intrincis (Geant coordinates) and Pad Time // etc are now done through a set of BrDetectorParamsTPC methods. This // hides all explicit references to geometry in this code, and makes it much // more aminable to introduction of dabatbase calls. // // January 2000 // Modifications made to deal with dynamic allocation of adc value rather than fixed // size. A compile flag DYNAMIC is used to enable this feature. // // All direct references to intrinsic to pad,bucket now removed. All is // converted through calls to BrDetectorParamsTPC methods. // Note that the drift now goes towards positive Y values (up), and that // the pads are counted from the postive X towards lower values. // // August 2000 // The non dynamic code is scheduled to be removed soon // // // System includes // #include "math.h" #include <stdlib.h> #ifdef UNIX #include <unistd.h> #endif #include <BrIostream.h> // // Root includes // #ifndef ROOT_TDirectory #include "TDirectory.h" #endif #include "TRandom.h" #include "TH1.h" #include "TROOT.h" #include "TCanvas.h" #include "TView.h" // // BRAT includes // #include "BrDigitizeTPC.h" #include "BrMath.h" #include "BrGeometryDbManager.h" #include "BrParameterDbManager.h" #include "BrDetectorParamsTPC.h" #include "BrDetectorVolume.h" #include "BrEventNode.h" #include "BrDigitizeTPC.h" #include "BrDataTable.h" #include "BrGeantHit.h" #include "BrTpcSequence.h" #include "BrUnits.h" #include "TClonesArray.h" //Stuff for Event Display #include "BrDetectorTPC.h" #define DYNAMIC //________________________________________________________________________ ClassImp(BrDigitizeTPC); //________________________________________________________________________ BrDigitizeTPC::BrDigitizeTPC() { // Default Constructor. Should only be used by ROOT internal. Use // BrDigitize(Title,name) instead. fParams_p = 0; fVolumeParams_p = 0; fDetectorNode = 0; fAdcVal = 0; fUseGongTimeResp = kFALSE; fUseHyperbolic = kFALSE; fTauScale = 0.5; // Must be adjusted empirically! fTshift = 2.5; // Probably OK fN0_pad_to_tot = 1.0; // default fAbsorp = 0.0; // default fMinNbins = 2; // corresponds to original, hardcoded value SetState(kSetup); } //________________________________________________________________________ BrDigitizeTPC::BrDigitizeTPC(const Char_t *Name, const Char_t *Title) : BrModule(Name, Title) { // Constructor to call from user code. // fParams_p = 0; fVolumeParams_p = 0; fDetectorNode = 0; fAdcVal =0; fAdcCut = 4.0; fNoise = 0.0; #ifdef WIN32 fDedxTon0 = (Float_t)(1.0 / (.045*0.001)/3.0); #else fDedxTon0 = (Float_t)(1.0 / (.045*BrUnits::keV)/3.0); #endif fUseGongTimeResp = kFALSE; fUseHyperbolic = kFALSE; fTauScale = 0.5; // Must be adjusted empirically! fTshift = 2.5; // Probably OK fN0_pad_to_tot = 1.0; // default fAbsorp = 0.0; // default fMinNbins = 2; // corresponds to original, hardcoded value SetState(kSetup); } //______________________________________________________________________________ void BrDigitizeTPC::ListDetectorParameters(){ // List the current value of the digitization parameters on // standard out. if( fParams_p != NULL) fParams_p->ListParameters(); else cout << "No parameters set for this Digitisation module" << endl; if(fVolumeParams_p) fVolumeParams_p ->ListParameters(); } //______________________________________________________________________________ void BrDigitizeTPC::ListEventStatistics() { BrModule::ListEventStatistics(); // cout<<"NumGeantHits = "<<fNumGeantHits << endl; } //______________________________________________________________________________ BrDigitizeTPC::~BrDigitizeTPC(){ if(fDetectorNode) delete fDetectorNode; if(fAdcVal){ for(int i=(fNpads-1);i>=0;i--) { delete [] fAdcVal[i]; } if (DebugLevel() > 0) cout << "delete index Table" << endl; delete [] fAdcVal; } delete [] fDiffused_array; delete [] fShaped_array; delete [] fFinal_array; delete [] fShaper_lookup_array; } //______________________________________________________________________________ void BrDigitizeTPC::Init() { // The init member function/method is used. Intended to be used for setting // up internals that cannot be done at the time of module construction or // after parameters has been modified. // Job-level initialisation SetState(kInit); // Get volume data with GeometryDbManager // Note the volume is not used anywhere except in the Draw. BrGeometryDbManager *gGeomDb = BrGeometryDbManager::Instance(); fVolumeParams_p = (BrDetectorVolume*) gGeomDb->GetDetectorVolume("BrDetectorVolume", GetName()); BrParameterDbManager *gParamDb = BrParameterDbManager::Instance(); fParams_p = (BrDetectorParamsTPC*) gParamDb->GetDetectorParameters("BrDetectorParamsTPC",GetName()); if(fAdcVal){ for(int i=(fNpads-1);i>=0;i--) { delete [] fAdcVal[i]; } if (DebugLevel() > 0) cout << "delete index Table" << endl; delete [] fAdcVal; } // Note the use of fNpads and fNbuckets as private data members of the class // The are only modifed though a call to Init() to ensure intregity of fAdcVal // 2D array. fNpads = fParams_p->GetPadsprow(); fNbuckets = fParams_p->GetNoBuckets(); Float_t** val = new Float_t * [fNpads]; for(int i=0;i<fNpads;i++){ val[i] = new Float_t [fNbuckets]; } fAdcVal = val; fDiffused_array = new Float_t [fNbuckets*10]; fShaped_array = new Float_t [fNbuckets*10]; fFinal_array = new Float_t [fNbuckets*10]; fShaper_lookup_array = new Float_t [fNbuckets*10]; } //_______________________________________________________________________________ void BrDigitizeTPC::DefineHistograms() { // // Book histograms (called from BrModule::Book) // // Define histograms. They are: // <fill in here> if (GetState() != kInit) { Stop("DefineHistograms", "Must be called after Init"); return; } TDirectory* saveDir = gDirectory; TDirectory* histDir = gDirectory->mkdir(Form("%s_Digitize", GetName())); histDir->cd(); hNGeantHits = new TH1F(Form("hNGeantHits%s", GetName()), Form("Geant Hits per event%s", GetName()), 100, 0., 500.); hNSequences = new TH1F(Form("hNSequences%s", GetName()), Form("Sequences per event%s", GetName()), 500, 0., 5000.); hNTimebins = new TH1F(Form("hNTimebins%s", GetName()), Form("Timebins per sequence%s", GetName()), 25, 0., 25.); hAdcDist = new TH1F(Form("hAdcDist%s", GetName()), Form("ADC distribution%s", GetName()), 1000, 0.,1000.); hAdcSumDist = new TH1F(Form("hAdcSumDist%s", GetName()), Form("ADC sum distribution%s", GetName()), 1000, 0.,5000.); hN0Dist = new TH1F(Form("hN0Dist%s", GetName()), Form("N0 e- distribution%s", GetName()), 1000, 0.,1000.); hSigmaTDist = new TH1F(Form("hSigmaTDist%s", GetName()), Form("Transverse cluster width%s", GetName()), 50, 0.,5.); hSigmaLDist = new TH1F(Form("hSigmaLDist%s", GetName()), Form("Longitudinal cluster width%s", GetName()), 50, 0.,5.); hSigmaLTotDist = new TH1F(Form("hSigmaLTotDist%s", GetName()), Form("Longitudinal width of shaped clusters%s", GetName()), 50, 0.,5.); hSigmaTRes = new TH1F(Form("hSigmaTRes%s", GetName()), Form("Transverse residual%s", GetName()), 50, -2.,2.); hSigmaLRes = new TH1F(Form("hSigmaLRes%s", GetName()), Form("Longitudinal residual%s", GetName()), 50, -2.,2.); hPadCentroid = new TH1F(Form("hPadCentroid%s", GetName()), Form("Centroid of pad%s", GetName()), fNpads, 0.,(Float_t)fNpads); hSequenceCentroid = new TH1F(Form("hSequenceCentroid%s", GetName()), Form("Centroid of sequence%s", GetName()), 155, 0.,155); hSequenceCentroidDecimal = new TH1F(Form("hSequenceCentroidDecimal%s", GetName()), Form("Centroid of sequence - integer%s", GetName()), 200, 0.,1.0); hYmiss = new TH1F(Form("hYmiss%s", GetName()), Form("Distance between hit and sequence centroid", GetName()), 200,-3.0,3.0); hXDistToHit = new TH1F(Form("hXDistToHit%s", GetName()), Form("Distance hit - pad centroid", GetName()), 200,-3.0,3.0); hYDistToHit = new TH1F(Form("hYDistToHit%s", GetName()), Form("Distance hit - time centroid", GetName()), 200,-3.0,3.0); hPadShape = new TProfile(Form("hPadShape%s", GetName()), Form("Shape of cluster in pad%s", GetName()), 20, -10.,10.,0.0,1.0); hShapeResp = new TProfile(Form("hShapeResp%s", GetName()), Form("Shaper response%s", GetName()), 100, 0., 100.,0.0,1.0); hTimeCluster2D = new TH2F(Form("hTimeCluster2D%s", GetName()), Form("Time shape of orig cluster%s", GetName()), 30, -15., 15.,100,0.0,1.0); hTimeShape2D = new TH2F(Form("hTimeShape2D%s", GetName()), Form("Time shape of folded cluster%s", GetName()), 30, -15., 15.,100,0.0,1.0); hPadShape2D = new TH2F(Form("hPadShape2D%s", GetName()), Form("Shape of cluster in pad%s", GetName()), 20, -10.,10.,100,0.0,1.0); hSequenceShape = new TProfile(Form("hSequenceShape%s", GetName()), Form("Shape of sequence%s", GetName()), 30, -15,15.,0.0,1.0); hSequenceShapeY2D = new TH2F(Form("hSequenceShapeY2D%s", GetName()), Form("Shape of sequence vs Y%s",GetName()),20,-10,10,30,-15,15); hSequenceShapeYsl2D = new TH2F(Form("hSequenceShapeYsl2D%s", GetName()), Form("Shape of sequence vs Ysl%s",GetName()),20,-1.0,1.0,50,-25,25); hSigmaTVsY = new TProfile(Form("hSigmaTVsY%s", GetName()), Form("Pad width vs drift length%s", GetName()),20,-10,10,0.0,5.0); hSigmaTVsXsl = new TProfile(Form("hSigmaTVsXsl%s", GetName()), Form("Pad width vs X slope%s", GetName()),20,-1.0,1.0,0.0,10.0); hSigmaLVsY = new TProfile(Form("hSigmaLVsY%s", GetName()), Form("Time width vs drift length%s", GetName()),20,-10,10,0.0,5.0); hSigmaLTotVsY = new TProfile(Form("hSigmaLTotVsY%s", GetName()), Form("Tot. time width vs drift length%s", GetName()),20,-10,10,0.0,5.0); hSigmaLTotVsY2D = new TH2F(Form("hSigmaLTotVsY2D%s", GetName()), Form("Tot. time width vs drift length%s", GetName()),20,-10,10,50,0.0,5.0); hSigmaLTotVsYsl = new TProfile(Form("hSigmaLTotVsYsl%s", GetName()), Form("Tot time width vs Y slope%s", GetName()),20,-1.0,1.0,0.0,10.0); hSigmaLTotVsYsl2D = new TH2F(Form("hSigmaLTotVsYsl2D%s", GetName()), Form("Tot. time width vs Y slope%s", GetName()),20,-1.0,1.0,50,0.0,10.0); hSigmaTResVsY = new TProfile(Form("hSigmaTResVsY%s", GetName()), Form("Trans resid vs drift length%s", GetName()),20,-10,10,0.0,2.0); hSigmaLResVsY = new TProfile(Form("hSigmaLResVsY%s", GetName()), Form("Long resid vs drift length%s", GetName()),20,-10,10,0.0,2.0); hSigmaTResVsXsl = new TProfile(Form("hSigmaTResVsXsl%s", GetName()), Form("Trans resid vs Xslope%s", GetName()),20,-0.5,0.5,0.0,2.0); hSigmaLResVsYsl = new TProfile(Form("hSigmaLResVsYsl%s", GetName()), Form("Long resid vs Yslope%s", GetName()),20,-0.5,0.5,0.0,2.0); hSigmaTResVsN0 = new TProfile(Form("hSigmaTResVsN0%s", GetName()), Form("Trans resid vs charge%s", GetName()),50,0,250,0.0,2.0); hSigmaLResVsN0 = new TProfile(Form("hSigmaLResVsN0%s", GetName()), Form("Long resid vs charge%s", GetName()),50,0,250,0.0,2.0); gDirectory = saveDir; } //______________________________________________________________________________ void BrDigitizeTPC::Event(BrEventNode* InputTable, BrEventNode* OutputTable) { // Event method to be called once per event. // This performes the actual slow simulation of the hits on the TPC. // // // // // The code has been changed (Feb 99) to be able to create sequences of // TPC data. The algorithm does not complete mimic the one presently // implemented in the // readout board. The conditions to create an adc cluster is // - #consequtive bins with AdcValue above fAdcCut must be two // - once a single bins goes below the cut the sequence is generated. // This later condition is strciter than the readoutboard which requires // 2 (3?) bins below the cut. // SetState(kEvent); // const Bool_t hyperbolic = kFALSE; #ifdef WIN32 const double pi = 3.1415926; const double twopi = 2.0*pi; #else const double pi = BrUnits::pi ; const double twopi = BrUnits::twopi; #endif const Float_t sqtwopi = (Float_t)BrMath::Sqrt(twopi); const Float_t sqtwo = Float_t(sqrt(2.)); Float_t a,b,z,ax,ay,n0,xlo,xhi,dxl,dxt,xc,yc; Float_t time_bucket, Anode_Gap, sigma_0, sigma_prf,path_length; Float_t dedx_cm, factor_to_channels, xwire, ywire; Float_t sigma_T,sigma_L,diffuse_T,diffuse_L,Ldrift; Float_t sigmadel_T,sigmadel_L; // behaviour of residuals Float_t y_time,tau,sigma_time,tau_norm,tshift; Float_t y_low,y_high,tnorm,pnorm,edge_low,edge_high,qpad,qbin; Int_t i; Int_t irow,ihit,ixmin,ixmax,min_time,max_time,ipad,itime; /* // Trying numerical integration for Gong response Float_t diffused_array[fNbuckets*10]; Float_t shaped_array[fNbuckets*10]; Float_t final_array[fNbuckets*10]; Float_t shaper_lookup_array[fNbuckets*10]; */ BrDetectorParamsTPC *par; Int_t NumHits; BrGeantHit* ghit_p; if(DebugLevel()>0){ printf("Inside Digitize of %sn",GetName());; if(DebugLevel()>0) InputTable->ListObjects(); } par = GetDetectorParamsTPC(); if(0 == par){ Warning("Event","Detector Parameters not set for %s",GetName()); return ; } if( DebugLevel() > 2){ cout << " Getting Parameters" << endl; par->ListParameters(); } // Decode the inputtable // Char_t TableName[32]; #ifndef OLD_CDAT_VERSION sprintf(TableName,"GeantHits %s", GetName()); #else sprintf(TableName,"GeantHits %sPR", GetName()); #endif fGeantHits = InputTable->GetDataTable(TableName); if(fGeantHits == 0 ) { if(DebugLevel() > 0){ cout << "No Table " << TableName << " : " << GetName() << endl; } return; } NumHits = fGeantHits->GetEntries(); if(NumHits == 0 ) { if(DebugLevel() > 0){ cout << "No Hits in " << TableName << " :" << GetName() << endl; } return; } /* cout << " About to fill hist hNGeantHits! " << endl; cout << " NumHits: " << NumHits << endl; if (!hNGeantHits) cout << " Histogram does not exist! " << endl; */ if (HistOn()){ hNGeantHits->Fill(NumHits); } // cout << " Filled hist hNGeantHits! " << endl; // // ObjectList returns TObjectArray // Prepare the output table // sprintf(TableName,"TPCSequence %s",GetName()); BrDataTable *SequenceTPC = new BrDataTable(TableName,"Sim"); OutputTable->AddObject(SequenceTPC); time_bucket = par->GetTimeBucket(); Anode_Gap = par->GetAnodeGap(); sigma_0 = par->GetShapingTime() /2.36; sigma_prf = par->GetAnodeGap(); tau = fTauScale*par->GetShapingTime(); // tshift = fTshift*tau; if (DebugLevel()>0) cout << " tau = " << tau << endl; // Calculate shaper lookup array Float_t shaper_integral = 0.0; Float_t shaper_centroid = 0.0; for (Int_t ifine=0; ifine<100; ifine++){ Float_t finetime = (Float_t)ifine*0.1*par->GetTimeBucket(); // in ns Float_t arg = finetime/tau; fShaper_lookup_array[ifine] = arg*arg*BrMath::Exp(-1.*arg); shaper_integral += fShaper_lookup_array[ifine]*0.1*par->GetTimeBucket(); shaper_centroid += fShaper_lookup_array[ifine]*0.1*par->GetTimeBucket()*ifine*0.1*par->GetTimeBucket(); if (DebugLevel()>2) cout << ifine << " Time: " << finetime << " Shape: " << fShaper_lookup_array[ifine] << endl; } shaper_centroid /= shaper_integral; if (DebugLevel()>2) cout << " Shaper integral: " << shaper_integral << " Centroid: " << shaper_centroid << endl; for (Int_t ifine=0; ifine<100; ifine++){ // Renormalize to integral 1 fShaper_lookup_array[ifine] = fShaper_lookup_array[ifine]/shaper_integral; if (HistOn()) hShapeResp->Fill(ifine,fShaper_lookup_array[ifine]); } factor_to_channels = par->GetADCGain() * par->GetADCChannels() / (Float_t)1024.; for(int actrow = 1; actrow <= par->GetNumberOfActiveRows(); actrow++) { irow = par->GetActiveRowNumber(actrow); Int_t NoTests = 0; { for(int pad=0;pad<fNpads;pad++) for(int time=0;time<fNbuckets;time++) fAdcVal[pad][time]=0.0; } for(ihit=0;ihit<NumHits;ihit++) { ghit_p = (BrGeantHit*)fGeantHits->At(ihit); if( irow == ghit_p->Isub()) { ax = ( ghit_p->LocalPosOut()[0] - ghit_p->LocalPosIn()[0] ) / ( ghit_p->LocalPosOut()[2] - ghit_p->LocalPosIn()[2] ); ay = ( ghit_p->LocalPosOut()[1] - ghit_p->LocalPosIn()[1] ) / ( ghit_p->LocalPosOut()[2] - ghit_p->LocalPosIn()[2] ); path_length = (Float_t)0.0; for(i=0;i<3;i++) path_length += (Float_t)pow(ghit_p->LocalPosOut()[i]-ghit_p->LocalPosIn()[i],2.); path_length = (Float_t)sqrt(path_length); dedx_cm = ghit_p->Dedx() / path_length; z = (ghit_p->LocalPosOut()[2] + ghit_p->LocalPosIn()[2]) / (Float_t)2.0; // dxt = par->GetPadLength() * BrMath::Abs(ax) / (Float_t)2.; // dxl = par->GetPadLength() * BrMath::Abs(ay) / (Float_t)2.; // Reproduce behaviour of tracks with large slope // Division with 2 - apparent bug dxt = par->GetPadLength() * BrMath::Abs(ax); dxl = par->GetPadLength() * BrMath::Abs(ay); xc = ghit_p->LocalPosIn()[0] + ax * ( z - ghit_p->LocalPosIn()[2] ); yc = ghit_p->LocalPosIn()[1] + ay * ( z - ghit_p->LocalPosIn()[2] ); Ldrift = par->IntrinsicToDrift(yc); diffuse_T = par->GetDtrans() * BrMath::Sqrt( (Double_t)Ldrift+2. ); diffuse_L = par->GetDlong() * BrMath::Sqrt( (Double_t)Ldrift+2. ); // As a temporary solution to the problem of angled tracks, // fold in the sigma of the square response of ax and ay // Treatment of residuals according to formulae: // del_T**2 = diffuse_T**2/(padlength*n0) + dxt**2/(12.*n0_per_pad) // del_L**2 = diffuse_L**2/(padlength*n0) + dxl**2/(12.*n0_per_pad) if (fUseGongTimeResp){ sigma_L = (Float_t)sqrt( diffuse_L*diffuse_L + dxl*dxl/(Float_t)12.); // Increase sigmadel_L! sigmadel_L = diffuse_L*diffuse_L/par->GetPadLength(); // cout << " Sigmadel_L 1: " << sigmadel_L; sigmadel_L += tau*tau*par->GetDriftv()*par->GetDriftv()/par->GetPadLength(); // cout << " Sigmadel_L 2: " << sigmadel_L; sigmadel_L += dxl*dxl*fN0_pad_to_tot/(Float_t)12.; // cout << " Sigmadel_L 3: " << sigmadel_L; sigmadel_L = (Float_t)sqrt(sigmadel_L); // cout << " Sigmadel_L: " << sigmadel_L << " Sigma_L: " << sigma_L << endl; // sigmadel_L = (Float_t)sqrt( diffuse_L*diffuse_L/par->GetPadLength() + dxl*dxl*fN0_pad_to_tot/(Float_t)12.); // for residuals } else { sigma_L = (Float_t)sqrt(pow(sigma_0*par->GetDriftv(),2) + diffuse_L*diffuse_L); sigmadel_L = (Float_t)sqrt( sigma_L*sigma_L/par->GetPadLength() + dxl*dxl*fN0_pad_to_tot/(Float_t)12.); // for residuals sigma_L = (Float_t)sqrt( sigma_L*sigma_L + dxl*dxl/(Float_t)12.); } sigma_T = (Float_t)sqrt(sigma_prf*sigma_prf + diffuse_T*diffuse_T); sigmadel_T = (Float_t)sqrt( sigma_T*sigma_T /par->GetPadLength() + dxt*dxt*fN0_pad_to_tot/(Float_t)12.); // for residuals // cout << " Sigmadel_T: " << sigmadel_T << " Sigma_T: " << sigma_T << endl; sigma_T = (Float_t)sqrt( sigma_T*sigma_T + dxt*dxt/(Float_t)12.); sigma_time = sigma_L/par->GetDriftv(); tau_norm = sigma_time/tau; if (HistOn()) hSigmaTDist->Fill(sigma_T/par->GetPadDistance()); if (HistOn()) hSigmaLDist->Fill(sigma_time/par->GetTimeBucket()); if (HistOn()) hSigmaTVsY->Fill(yc,sigma_T/par->GetPadDistance()); if (HistOn()) hSigmaLVsY->Fill(yc,sigma_time/par->GetTimeBucket()); if (HistOn()) hSigmaTVsXsl->Fill(ax,sigma_T/par->GetPadDistance()); if(DebugLevel()>2) cout << par->GetDlong() << " " << Ldrift << " " << diffuse_L << endl; /* Use a 3 sigma limit on searching for pads which might * have been hit At this point we have all the cluster * variables needed ax = tan(alpha), ay=tan(lambda) yc = * Ldrift, xc which are positions before evaluating the * diffusion in drifts and response widths */ n0 = dedx_cm * fDedxTon0; n0 *= par->GetPadLength() * (Float_t)sqrt( 1.0 + ax*ax + ay*ay ); n0 *= BrMath::Exp(-1.*fAbsorp*Ldrift); if (HistOn()) hN0Dist->Fill(n0); // Evaluate the position on the Anode wires.. xwire = par->IntrinsicToDistance(xc); ywire = Ldrift; // Distribute signal according to n0 on the Pad gRandom->Rannor(a,b); // xwire += sigma_T / (Float_t)BrMath::Sqrt(n0) * a; // ywire += sigma_L / (Float_t)BrMath::Sqrt(n0) * b; xwire += sigmadel_T / (Float_t)BrMath::Sqrt(n0) * a; ywire += sigmadel_L / (Float_t)BrMath::Sqrt(n0) * b; if (HistOn()) hSigmaTRes->Fill(sigmadel_T / (Float_t)BrMath::Sqrt(n0) * a); if (HistOn()) hSigmaLRes->Fill(sigmadel_L / (Float_t)BrMath::Sqrt(n0) * b); if (HistOn()) hSigmaTResVsY->Fill(yc,BrMath::Abs(sigmadel_T / (Float_t)BrMath::Sqrt(n0) * a)); if (HistOn()) hSigmaLResVsY->Fill(yc,BrMath::Abs(sigmadel_L / (Float_t)BrMath::Sqrt(n0) * b)); if (HistOn()) hSigmaTResVsXsl->Fill(ax,BrMath::Abs(sigmadel_T / (Float_t)BrMath::Sqrt(n0) * a)); if (HistOn()) hSigmaLResVsYsl->Fill(ay,BrMath::Abs(sigmadel_L / (Float_t)BrMath::Sqrt(n0) * b)); if (HistOn()) hSigmaTResVsN0->Fill(n0,BrMath::Abs(sigmadel_T / (Float_t)BrMath::Sqrt(n0) * a)); if (HistOn()) hSigmaLResVsN0->Fill(n0,BrMath::Abs(sigmadel_L / (Float_t)BrMath::Sqrt(n0) * b)); Float_t padxwire = par->IntrinsicToPad(xwire); if (HistOn()) hPadCentroid->Fill(padxwire); // Evaluate the pad response xlo = xwire - (Float_t)3.0 * sigma_T; xhi = xwire + (Float_t)3.0 * sigma_T; if( DebugLevel() > 2) { cout << "Padrow #" << irow <<endl; cout.width(10); cout << " DT,DL " << setw(8) << diffuse_T << "," << setw(8)<<diffuse_L << endl; cout << " n0 " << setw(8)<< n0 << dedx_cm<<fDedxTon0<<endl; cout << " " << xwire<<" " << xlo << " " << xhi << endl; } // Find pads which must be looked at. ixmin = (Int_t) par->DistanceToPad(xlo); //(Int_t)( xlo / par->GetPadDistance() ); // + 1; ixmin = BrMath::Max(ixmin,0); ixmax = (Int_t) par->DistanceToPad(xhi); // (Int_t)( xhi / par->GetPadDistance() ); // + 1; ixmax = BrMath::Min(fNpads-1,ixmax); // Pad numbers are counted 0,1,2.. from lower X // values. there is the same number of Pads on each side of // the centerline longitudinal Diffusion i.e. along Drift // Convert into time buckets. Do fold the response of the // shaper into this if (fUseGongTimeResp){ y_low = ywire - (Float_t)5.0 * sigma_L; y_high = ywire + (Float_t)10.0 * sigma_L; y_time = par->DriftToTime(irow, ywire); } else { y_low = ywire - (Float_t)3.0 * sigma_L; y_high = ywire + (Float_t)3.0 * sigma_L; y_time = par->DriftToTime(irow, ywire); } min_time = (Int_t) par->DriftToTime(irow, y_low); min_time = BrMath::Max(min_time,0); max_time = (Int_t) par->DriftToTime(irow, y_high); max_time = BrMath::Min(max_time, fNbuckets-1); if( DebugLevel() > 2){ cout.width(10); cout << " " << ixmin << " "<< ixmax << endl; cout << " " << y_low << " " << y_high << " yc: " << yc << endl; cout << " " << min_time << " " << max_time << " ytime: " << y_time << endl; cout << " Diffuse_L " << diffuse_L << " Sigma_L " << sigma_L << " sigmadel_L " << sigmadel_L << endl; cout << " Sigma_time " << sigma_time << endl; } tnorm = (Float_t)1. / ( sqtwopi * sigma_L ); pnorm = (Float_t)1. / ( sqtwopi * sigma_T ); if (fUseGongTimeResp){ // Calculate time response - same for all pads if (DebugLevel()>2) cout << " ytime: " << y_time << " sigma_time: " << sigma_time << endl; for (Int_t ifinetime = 0;ifinetime<=fNbuckets*10;ifinetime++ ) { fDiffused_array[ifinetime] = 0.0; fShaped_array[ifinetime] = 0.0; fFinal_array[ifinetime] = 0.0; } // First calculate Gaussian shape for arriving cluster Float_t original_integral = 0.0; Float_t original_centroid = 0.0; for (Int_t ifinetime = min_time*10;ifinetime<=max_time*10;ifinetime++ ) { Float_t timediff = ((Float_t)ifinetime*0.1 - y_time)*par->GetTimeBucket(); Float_t timearg = timediff/sigma_time; fDiffused_array[ifinetime] = BrMath::Exp(-0.5*timearg*timearg)/ (sigma_time*sqtwo); original_integral += fDiffused_array[ifinetime]*0.1*par->GetTimeBucket(); original_centroid += fDiffused_array[ifinetime]*0.1*par->GetTimeBucket()*ifinetime*0.1*par->GetTimeBucket(); // if (NoTests<1) cout << " timediff " << timediff << " diffuse shape " << fDiffused_array[ifinetime] << endl; } original_centroid /= original_integral; // if (NoTests<10) cout << " Centroid: " << original_centroid << " Integral: " << original_integral << endl; // renormalize Gaussian for (Int_t ifinetime = min_time*10;ifinetime<=max_time*10;ifinetime++ ) { fDiffused_array[ifinetime] = fDiffused_array[ifinetime] / original_integral; Float_t bin_now = (Float_t)(ifinetime)/10.0 - y_time; if (HistOn()) hTimeCluster2D->Fill(bin_now,fDiffused_array[ifinetime]*200); // if (NoTests<1) cout << " timech " << ifinetime - y_time*10 << " diffuse shape " << fDiffused_array[ifinetime] << endl; } Float_t shaped_integral = 0.0; Float_t shaped_centroid = 0.0; // Fold cluster with shaper response function for (Int_t ifinetime = min_time*10;ifinetime<=(100+max_time*10);ifinetime++ ) { for (Int_t ishape = ifinetime; ishape<ifinetime+100;ishape++ ) { fShaped_array[ishape] += fDiffused_array[ifinetime]*fShaper_lookup_array[ishape-ifinetime]*0.1*par->GetTimeBucket(); ; } } ; for (Int_t ifinetime = min_time*10;ifinetime<=(100+max_time*10);ifinetime++ ) { shaped_centroid += fShaped_array[ifinetime]*0.1*par->GetTimeBucket()*ifinetime*0.1*par->GetTimeBucket(); shaped_integral += fShaped_array[ifinetime]*0.1*par->GetTimeBucket(); } shaped_centroid /= shaped_integral; // if (NoTests<10) cout << " Shaped integral: " << shaped_integral << " Centroid: " << shaped_centroid << " Shift: " << shaped_centroid - original_centroid << endl; Float_t shift = (shaped_centroid - original_centroid)/(0.1*par->GetTimeBucket()); Int_t ishift = (Int_t)shift; Float_t decshift = shift - (Float_t)ishift; // if (NoTests<10) cout << " Shift " << shift << " Ishift: " << ishift << " Decshift: " << decshift << endl; Float_t finalCentroid = 0.0; Float_t finalIntegral = 0.0; for (Int_t ifinetime = min_time*10;ifinetime<=max_time*10;ifinetime++ ) { // fFinal_array[ifinetime] = (1.0-decshift)*fShaped_array[ifinetime+ishift] + decshift*fShaped_array[ifinetime+ishift+1]; fFinal_array[ifinetime] = fShaped_array[ifinetime+ishift] + decshift*(fShaped_array[ifinetime+ishift+1]-fShaped_array[ifinetime+ishift]); finalIntegral += fFinal_array[ifinetime]*0.1*par->GetTimeBucket(); finalCentroid += fFinal_array[ifinetime]*0.1*par->GetTimeBucket()*ifinetime*0.1*par->GetTimeBucket(); Float_t bin_now = (Float_t)(ifinetime)/10.0 - y_time; if (HistOn()) hTimeShape2D->Fill(bin_now,fShaped_array[ifinetime]*100); //if (NoTests<1) cout << " timech " << ifinetime - y_time*10 << " Timehape " << fFinal_array[ifinetime] << endl; // if (NoTests<1) cout << " Final array: " << fFinal_array[ifinetime] << " Shaped array: " << fShaped_array[ifinetime] << endl; } // if (NoTests<10) cout << " Final integral: " << finalIntegral << " Centroid: " << finalCentroid << endl; NoTests++; } // if Gong response Float_t xqsum = 0.0; Float_t xpadsum = 0.0; for( ipad = ixmin;ipad<=ixmax;ipad++ ) { Int_t Ntests = 0; edge_low = ((Float_t)ipad + (Float_t).5)*par->GetPadDistance() - par->GetPadWidth()/(Float_t)2.; edge_high = ((Float_t)ipad + (Float_t).5)*par->GetPadDistance() + par->GetPadWidth()/(Float_t)2.; if( fUseHyperbolic ) { printf("This part has not been tested yet!!!"); qpad = (Float_t)atan((Float_t)sinh(pi*(edge_high-xwire) / (Float_t)2./Anode_Gap)) - (Float_t)atan((Float_t)sinh(pi*(edge_low -xwire) / (Float_t)2./Anode_Gap)); xqsum += qpad; xpadsum += (Float_t)ipad*qpad; if (HistOn()) hPadShape->Fill((Float_t)ipad - padxwire, qpad); if (HistOn()) hPadShape2D->Fill((Float_t)ipad - padxwire, qpad); } else { qpad = pnorm*(Float_t)0.5*(BrMath::Erf((edge_high-xwire) / (sigma_T*sqtwo)) - BrMath::Erf((edge_low-xwire) / (sigma_T*sqtwo))); xqsum += qpad; xpadsum += (Float_t)ipad*qpad; if (HistOn()) hPadShape->Fill((Float_t)ipad - padxwire, qpad); if (HistOn()) hPadShape2D->Fill((Float_t)ipad - padxwire, qpad); } Float_t centroid_test = 0.0; Float_t centroid2_test = 0.0; Float_t integral_test = 0.0; for( itime = min_time;itime<=max_time;itime++ ) { if (fUseGongTimeResp){ //Try realistic time response Float_t qsum = 0; for (Int_t ifinetime = itime*10-5; ifinetime<(itime*10+5); ifinetime++){ if (ifinetime>=min_time*10) qsum += fFinal_array[ifinetime]*0.1*par->GetTimeBucket(); } integral_test += qsum; centroid_test += qsum*itime*par->GetTimeBucket(); centroid2_test += qsum*itime*par->GetTimeBucket()*itime*par->GetTimeBucket(); // if (Ntests<1) cout << " itime: " << itime << " qsum: " << qsum << endl; Ntests++; /* Double_t lambda = ((y_time - itime)*par->GetTimeBucket()-tshift)/sigma_time + tau_norm; Double_t arg1 = tau_norm*(lambda - 0.5*tau_norm); Double_t fac1 = BrMath::Exp(arg1); Double_t gau1 = BrMath::Exp(-0.5*lambda*lambda); gau1 *= lambda/BrMath::Sqrt(2.0*pi); Double_t erf1 = BrMath::Erf(lambda/sqtwo); Double_t fac2 = (1.0 + lambda*lambda)/2.0*(1 - erf1) - gau1; Double_t timeresp = tau_norm*tau_norm*fac1*fac2; */ Double_t timeresp = qsum; if (HistOn()) hSequenceShape->Fill((Float_t)itime-y_time,timeresp); if (HistOn()) hSequenceShapeY2D->Fill(yc,(Float_t)itime-y_time); if (HistOn()) hSequenceShapeYsl2D->Fill(ay,(Float_t)itime-y_time); qbin = timeresp*qpad*n0; } else { // use Gaussian edge_low = par->TimeToDrift(itime); edge_high = par->TimeToDrift(itime+1); qbin = (Float_t)0.5 * (BrMath::Erf((edge_high - ywire) / (sigma_L*sqtwo)) - BrMath::Erf((edge_low - ywire) / (sigma_L*sqtwo))); if (HistOn()) hSequenceShape->Fill((Float_t)itime-y_time,qbin); if (HistOn()) hSequenceShapeY2D->Fill(yc,(Float_t)itime-y_time); if (HistOn()) hSequenceShapeYsl2D->Fill(ay,(Float_t)itime-y_time); qbin *= tnorm * qpad * n0; } // cout << " ipad: " << ipad << " itime: " << itime << " qbin: " << qbin << endl; fAdcVal[ipad][itime]+=qbin*factor_to_channels; } //loop over itime centroid_test /= integral_test; centroid2_test /= integral_test; Float_t time_centroid = centroid_test/par->GetTimeBucket(); Float_t time_centroid2 = centroid2_test/(par->GetTimeBucket()*par->GetTimeBucket()); Float_t timevar = BrMath::Sqrt(time_centroid2-time_centroid*time_centroid); if (HistOn()) hYmiss->Fill(time_centroid-y_time); if (HistOn()) hSigmaLTotDist->Fill(timevar); if (HistOn()) hSigmaLTotVsY2D->Fill(yc,timevar); if (HistOn()) hSigmaLTotVsY->Fill(yc,timevar); if (HistOn()) hSigmaLTotVsYsl->Fill(ay,timevar); if (HistOn()) hSigmaLTotVsYsl2D->Fill(ay,timevar); Float_t y_resulting = par->TimeToIntrinsic(irow,time_centroid); if (HistOn()) hYDistToHit->Fill(y_resulting-yc); // if (NoTests<10) cout << " Centroid: " << centroid_test << " Integral: " << integral_test << " Ytime: " << time_centroid << endl; } // loop over ipad Float_t xpadcentroid = xpadsum / xqsum; Float_t x_resulting = par->PadToIntrinsic(xpadcentroid); if (HistOn()) hXDistToHit->Fill(x_resulting-xc); } // correct row } // loop over ihit // // Add noise if fNoise non-zero // if(fNoise>0.0) { for(ipad=0;ipad<fNpads;ipad++) for(itime=0;itime<fNbuckets;itime++){ gRandom->Rannor(a,b); fAdcVal[ipad][itime]+=a*fNoise; } } if(DebugLevel()> 0) cout << "Scan for Sequences Row " << irow << endl; for(int pad=0; pad<fNpads; pad++) { if( ! par->IsPadActive(irow, pad)) continue; Bool_t insequence=kFALSE; int timebin; int nbins=0; for(int time=0; time < fNbuckets; time++){ if(fAdcVal[pad][time] > fAdcCut){ if(insequence){ nbins++; } else { nbins=1; timebin = time; insequence=kTRUE; } } else { if(insequence){ // if(nbins>1){ if(nbins>=fMinNbins){ BrTpcSequence* seq = new BrTpcSequence(nbins); SequenceTPC->Add(seq); seq->SetRow(irow); seq->SetPad(pad); seq->SetTime(timebin); Short_t* adc = seq->GetSequence(); Float_t adcsum = 0.0; Float_t adcbinsum = 0.0; for(int it=0;it<nbins;it++, timebin++){ Short_t digval = (Short_t) fAdcVal[pad][timebin]; *adc++ = Short_t ( digval > 1024 ? 1024 : digval); if (HistOn()) hAdcDist->Fill(digval); adcsum += digval; adcbinsum += digval*(Float_t)timebin; } if (HistOn()) hNTimebins->Fill(nbins); Float_t seq_centroid = adcbinsum/adcsum; Float_t seq_cent_dec = seq_centroid - (Int_t)seq_centroid; // if (DebugLevel()) cout << " Cent: " << seq_centroid << " Dec cent: " << seq_cent_dec << endl; if (HistOn()) hAdcSumDist->Fill(adcsum); if (HistOn()) hSequenceCentroid->Fill(seq_centroid); if (HistOn()) hSequenceCentroidDecimal->Fill(seq_cent_dec); // if (HistOn()) hYmiss->Fill(seq_centroid-y_time); } insequence = kFALSE; } } } // time if(insequence){ // if(nbins>1){ if(nbins>=fMinNbins){ BrTpcSequence* seq = new BrTpcSequence(nbins); SequenceTPC->Add(seq); seq->SetRow(irow); seq->SetPad(pad); seq->SetTime(timebin); Short_t* adc = seq->GetSequence(); Float_t adcsum = 0.0; Float_t adcbinsum = 0.0; for(int it=0;it<nbins;it++, timebin++){ Short_t digval = (Short_t) fAdcVal[pad][timebin]; *adc++ = Short_t ( digval > 1024 ? 1024 : digval); if (HistOn()) hAdcDist->Fill(digval); adcsum += digval; adcbinsum += digval*(Float_t)timebin; } if (HistOn()) hNTimebins->Fill(nbins); Float_t seq_centroid = adcbinsum/adcsum; Float_t seq_cent_dec = seq_centroid - (Int_t)seq_centroid; // if (DebugLevel()) cout << " Cent: " << seq_centroid << " Dec cent: " << seq_cent_dec << endl; if (HistOn()) hAdcSumDist->Fill(adcsum); if (HistOn()) hSequenceCentroid->Fill(seq_centroid); if (HistOn()) hSequenceCentroidDecimal->Fill(seq_cent_dec); // if (HistOn()) hYmiss->Fill(seq_centroid-y_time); } } } } //loop over irow if (HistOn()) hNSequences->Fill(SequenceTPC->GetEntries()); //Get stuff for statistics if(DebugLevel()>0) cout << "Total number of sequences " << SequenceTPC->GetEntries() << endl; } //_________________________________________________________ void BrDigitizeTPC::DrawGeantHits() { // // BrDetectorVolume *vol; Char_t volbuf[64]; Char_t canvasname[64]; Char_t nodename[64]; TCanvas *canvas; sprintf(canvasname,"%sCanvas",GetName()); if(!fDetectorNode) { canvas = (TCanvas*)gROOT->FindObject(canvasname); if(!canvas) canvas = new TCanvas(canvasname,canvasname,50,50,400,400); sprintf(volbuf,"TPC %s",GetName()); vol = GetDetectorVolume(); sprintf(nodename,"%sDetectorNode",GetName()); fDetectorNode = new BrDetectorTPC(nodename,nodename,GetName()); //vol); fDetectorNode->Draw(); TView *view = canvas->GetView(); // Set the view to dimensions we know for TPC's; this could be made more general // but following looks pretty good. view->SetRange(-35,-35,-35,35,35,35); Int_t irep; view->SetView(-90,90,0,irep); } fDetectorNode->Clear(); canvas = (TCanvas*)gROOT->FindObject(canvasname); //Check if canvas still exists; if not, don't draw the hits. if(canvas) { canvas->cd(); fDetectorNode->DrawGeantHits(fGeantHits); } } ///////////////////////////////////////////////////////////////////////////// // // $Log: BrDigitizeTPC.cxx,v $ // Revision 1.5 2001/11/13 21:11:24 trulsml // Put if (debuglevel) in front of annoying: cout << "delete index Table" << endl; // in the destructor and the Init method. // // Revision 1.4 2001/09/18 19:31:56 trine // Improved description of angled tracks and of residuals. // Added members for setting: // - min. no of bins in a BrTpcSequence // - ratio of total no. of electrons / electrons per pad. // // Revision 1.3 2001/08/12 13:32:58 cholm // Changed CTOR prototype to comply with BrModule CTOR (const Char_t*'s instead // of just Char_t*'s). // // Revision 1.2 2001/06/22 17:40:08 cholm // Changes t odo with data classes renaming. // // Revision 1.1.1.1 2001/06/21 14:55:06 hagel // Initial revision of brat2 // // Revision 1.40 2001/06/19 18:57:43 trine // Added possibility to use more realistic time response taken from // Gong's STAR note, instead of Gaussian default. // Moved calculation of "residuals" to after place in code where track // slopes are folded into cluster widths. // // Revision 1.38 2001/04/18 20:19:42 videbaek // initilaize fAdcVal in constructor // // Revision 1.37 2000/11/22 21:22:35 videbaek // Modified for Init handling. // Some cosmetic changes. // // Revision 1.36 2000/10/02 17:07:56 pchristi // Cleaned up headers and source files in the tpc dir // // Revision 1.35 2000/09/07 20:05:29 videbaek // Add methods for row setting, and Modified to proper row numbers 1,2... // // Revision 1.34 2000/09/04 14:58:40 videbaek // Add some redundant block brackets to make VC++6.0 happy. // // Revision 1.33 2000/08/28 18:05:54 videbaek // Move CVS log to end of file // // Revision 1.32 2000/08/28 00:18:45 videbaek // Added calibration parameter to Detectorparams, to deal better with time-distance conversion // in TPCs. Added row dependent conversion too. // // Revision 1.31 2000/04/28 21:22:13 videbaek // Add new classes for TPC. Many changes and not completely debugged. // Note that the old method BrTPCLocaltracking has been tested to work. // // Revision 1.30 2000/03/21 21:22:34 cholm // Several changes: A few hacks where needed to compile on Digital Unix, noticably in my - Christian Holm - classes for the DB and Exceptions. Proberly still need to do something to some of Konstantins stuff. Also, I removed a lot of warnings fromt the compiler, by teddying up the code. Several places, old C-style indecies in for loops were assumed. Corrected. Unused variables have been commented out. // // Revision 1.29 2000/01/18 18:39:54 videbaek // Moved covariance array form fit to calling parameter. // lized and checkout of digitization, reconstruction code. // Added (possible noise ) to digitization. Changed default value of // dedxton0. Made it parameter in DetectorParamsTPC // // Revision 1.28 2000/01/04 21:49:38 videbaek // start adding dynamic arrays rather than fixed size. // Change cluster size max to deal with noisy cosmic ray events. // // Revision 1.27 1999/12/30 19:45:00 videbaek // Code changes to use the intrinsic cocrdinates consistently. // // Revision 1.26 1999/12/18 21:37:42 videbaek // Obsolete for a long timeBrDigTTPC.cxx // // Revision 1.25 1999/08/14 17:01:59 videbaek // Updates for // a) RawData adding the Trigger TDCs' // b) Implied new data objects BrTrigBB // c) Clean up of TPC code // d) Test MyMonitor and scripts // // Revision 1.24 1999/03/09 18:38:34 videbaek // Minor fixes. Now a proper working version. Included dependeces in // Makefile. Added some member function to the tracking class, and // minor updates. // // Revision 1.23 1999/03/07 23:28:26 hagel // Slight modifications in SetDetectorParamsTPC to prevent crashing if user // uses it since BrParameterDbManager is implemented. // // Revision 1.22 1999/03/07 00:00:45 hagel // // // Revision 1.21 1999/02/28 22:14:38 hagel // Changes to make event display more robust and seamless // // Revision 1.20 1999/02/25 14:53:16 videbaek // Inserted code for using BrTPCSequnece instead of BrDigTPC for both // TPC digitization and the TPC local track reconstruction.// // Revision 1.19 1999/01/28 21:28:13 videbaek // Added BrTPSSqquences to libraries. Not yet in use, though. Expected to // replace BrDigTPC objects // Added cvs flags // Changed name of BrDigTPC data tables. // // Revision 1.18 1999/01/15 16:35:48 videbaek // Changes included using clonesarrays in digitizations. // Localtracking added cuts for type 4 clusters // // Revision 1.17 1999/01/06 23:08:04 hagel // Add support for DrawGeantTracks in Digitize // // Revision 1.16 1998/12/21 20:27:41 videbaek // Improved geometry for MTP1 and MTP2 // // Revision 1.14 1998/10/09 19:23:59 videbaek // Add some diagnostics output // // Revision 1.13 1998/09/27 17:11:02 alv // removed EventStatisticsStart and EventStatisticsEnd // (inherit them from BrModule) // In ListEventStatistics // used BrModule::ListEventStatistics to list CPU time // // Revision 1.12 1998/09/18 20:37:40 videbaek // veMove SetRowPosition from being inline // // Revision 1.11 1998/09/16 20:17:26 alv // In constructors // Added fParams_p=0 and fVolumeParams_p=0 // In SetDetectorParamsTPC // Protected against delete of par // Used copy constructor insted of new+assignment // In Event // Removed embedded x0 in format string of sprintf // Removed unused variables twopi and sqpi // // Revision 1.10 1998/08/05 19:34:27 hagel // Fix BrDigTPC memory leak // // Revision 1.9 1998/07/02 13:58:43 videbaek // comment updates mainly // // Revision 1.8 1998/05/14 21:18:15 videbaek // update makefile - do not reload .so libraries // // Revision 1.7 1998/05/13 20:30:50 hagel // Changes to compile under Solaris // // Revision 1.6 1998/05/08 20:45:56 videbaek // TPC class updates; added functionality // // Revision 1.5 1998/05/08 00:58:43 videbaek // TPC code updates // // Revision 1.4 1998/04/30 16:40:48 hagel // Changes to work with new version of GBRAHMS // // Revision 1.3 1998/04/06 21:12:00 videbaek // Clean up and additions for Win95 // // Revision 1.2 1998/04/01 22:08:36 hagel // Add Local tracking and fix bugs in Digitization // // Revision 1.1.1.1 1998/03/04 21:34:03 brahmlib // Brat tpc // // |
||||||
|
This page automatically generated by script docBrat by Christian Holm |
Copyright ; 2002 BRAHMS Collaboration
<brahmlib@rcf.rhic.bnl.gov>
|