|
//////////////////////////////////////////////////////////////////////////// // // BrDigitizeBB // // Digitisation Module class for Beam Beam Counters.. // The module expect to get geant hits as input and produce digitized // BB data BrBbDig. The table is labelled Raw as secondary lable. // The method for digitization is described in the Event Method. // // // //////////////////////////////////////////////////////////////////////////// // // // Revisions: // April 1 1999, fv // Included BrParametersDBmanager for better detector parameter handling // Changed SetDetectorParamsDB to be private (scheduled for deletion) // Removed at the same time the delete fParamsDB which now is done by // the database manager // November 28, 1999 fv // Fix serious error. The hits were not filled with proper geant // information but from some other hit in another module. // Use newly added method ConvertGeantModuleNo(,) for assigning // Module Numbers // // // $Id: BrDigitizeBB.cxx,v 1.3 2001/08/12 13:28:15 cholm Exp $ // $Author: cholm $ // $Date: 2001/08/12 13:28:15 $ #include "math.h" #include <stdlib.h> #ifdef UNIX #include <unistd.h> #endif #include <cassert> #include <BrIostream.h> #ifndef ROOT_TMath #include "TMath.h" #endif # include <TRandom.h> # include <TDirectory.h> # include <TH1.h> # include "BrGeometryDbManager.h" # include "BrParameterDbManager.h" # include "BrDigitizeBB.h" # include "BrDetectorParamsBB.h" # include "BrEventNode.h" # include "BrDataTable.h" #include <BrGeantHit.h> #include <BrGeantTrack.h> # include "BrBbDig.h" # include "BrVector3D.h" //__________________________________________________________________ ClassImp(BrDigitizeBB); //__________________________________________________________________ BrDigitizeBB::BrDigitizeBB() { // default // Constructor. fParamsBB_p = 0; } //____________________________________________________________________ BrDigitizeBB::BrDigitizeBB(const Text_t *Name, const Char_t *Title) :BrModule(Name,Title) { // Normal Named constructor // } //____________________________________________________________________ void BrDigitizeBB::DefineHistograms() { // define histograms to look at background contributions // // // Book histograms (called from BrModule::Book) // const Int_t nt = 200; // number of TDC channels in histogram const Axis_t tl = 250; // lower limit for TDCs const Axis_t th = tl+nt; // upper limit for TDCs const Int_t na = 200; // number of ADC channels in histogram const Axis_t al = 0; // lower limit for ADCs const Axis_t ah = al+na*30; // upper limit for ADCs if (!fParamsBB_p) { Warning("DefineHistograms","Number of tubes unknown."); return; } // Create new subdirectory for histograms TDirectory* savdir = gDirectory; // remember current directory hdir = savdir->mkdir(GetName()); if (!hdir) { Warning("DefineHistograms","Could not create histogram subdirectory"); return; } hdir->cd(); const Int_t NoLeftTubes = fParamsBB_p->GetNoLeftTubes(); const Int_t NoRightTubes = fParamsBB_p->GetNoRightTubes(); hltdc = new TH1F*[NoLeftTubes+1]; hladc = new TH1F*[NoLeftTubes+1]; hpladc = new TH1F*[NoLeftTubes+1]; hrtdc = new TH1F*[NoRightTubes+1]; hradc = new TH1F*[NoRightTubes+1]; hpradc = new TH1F*[NoRightTubes+1]; hltdc[0] = new TH1F("hltdc_0", "Time Left all", nt, tl, th); hladc[0] = new TH1F("hladc_0", "ADC Left all", 2*na, al, ah); for (Int_t TubeNo = 1; TubeNo <= NoLeftTubes; TubeNo++) { Char_t hname[20]; Char_t htitle[80]; sprintf(hname , "hltdc_%d" , TubeNo); sprintf(htitle, "Tdc Left %d", TubeNo); hltdc[TubeNo] = new TH1F(hname, htitle, nt, tl, th); sprintf(hname , "hladc_%d" , TubeNo); sprintf(htitle, "Adc Left %d", TubeNo); hladc[TubeNo] = new TH1F(hname, htitle, na, al, ah); sprintf(hname , "hpladc_%d" , TubeNo); sprintf(htitle, "primary Adc Left %d", TubeNo); hpladc[TubeNo] = new TH1F(hname, htitle, na, al, ah); } hrtdc[0] = new TH1F("hrtdc_0", "Time Right all", nt, tl, th); hradc[0] = new TH1F("hradc_0", "ADC Right all", 2*na, al, ah); for (Int_t TubeNo = 1; TubeNo <= NoRightTubes; TubeNo++) { Char_t hname[20]; Char_t htitle[80]; sprintf(hname , "hrtdc_%d" , TubeNo); sprintf(htitle, "Tdc Right %d", TubeNo); hrtdc[TubeNo] = new TH1F(hname, htitle, nt, tl, th); sprintf(hname , "hradc_%d" , TubeNo); sprintf(htitle, "Adc Right %d", TubeNo); hradc[TubeNo] = new TH1F(hname, htitle, na, al, ah); sprintf(hname , "hpradc_%d" , TubeNo); sprintf(htitle, "primary Adc Right %d", TubeNo); hpradc[TubeNo] = new TH1F(hname, htitle, na, al, ah); } hlhits = new TH1F("hlhits", "Hits Left ", 50, -0.0, 49.5); hltubes = new TH1F("hltubes", "Tubes Left ", NoLeftTubes,1 ,NoLeftTubes); hrhits = new TH1F("hrhits", "Hits Right ", 50, -0.5, 49.5); hrtubes = new TH1F("hrtubes", "Tubes Right",NoRightTubes, 1, 1+NoRightTubes); // Restore directory savdir->cd(); } //___________________________________________________________________________ BrDigitizeBB::~BrDigitizeBB() { }; //___________________________________________________________________________ void BrDigitizeBB::Init() { // Initialize entry point. BrParameterDbManager *gParamDb = BrParameterDbManager::Instance(); fParamsBB_p = (BrDetectorParamsBB*) gParamDb->GetDetectorParameters("BrDetectorParamsBB",GetName()); } //___________________________________________________________________________ void BrDigitizeBB::Event(BrEventNode* InputTable, BrEventNode* OutputTable) { // Event method to be called once per event. // This performes the actual slow simulation // of the hits on the BB. // // Description of method: // Each tube can in fact be hit by multiple hits // 1) First a loop over all hits is made to see how many // hits there are per tube. If the is only a single (charged) hit // The adc and tdc value is calculated from the position, // path lenght and flight time // 2) For multiple hits each end of the tube is dealt with // seperately. The fastest signals is found (at each end). If the // following hits are within the resolving time / gate width // (100n sec) the adc values are added to thos from the first hit. // BrDetectorParamsBB *par; Int_t NumHits; BrDataTable* GeantHits; BrGeantHit* ghit_p; BrBbDig *digbb_p; if(DebugLevel() > 4){ cout << "Inside Digitize of " << GetName() << endl; // InputTable->ListObjects(); } // // Get simulation parameters // par = GetDetectorParamsBB(); if(par == NULL){ Warning("event","parameters are not set"); return ; } // if( DebugLevel() > 2){ // cout << " Getting Digitization Parameters" << endl; // par->ListParameters(); //} // // Decode the inputtable // The Beam Beam counter hits come from 4 tables namely the // right array and the left each with two different kinds of tubes // "GeantHits PURA", PURB PULA PULB. // // An identical peice of code is executed for each of the 4 possibilities // with slightly different action (in terms of storage) // The output is a BrDataObject with two tables one for the left array hits // and one for the right array hits; // Char_t TableName[32]; int option; // int Left_hits = 0; // int Right_hits = 0; int notubes; // No of Tubes for a given kind Bool_t Right; int firstId; // Prepare the output table // BrDataTable *DigitizedBB_Right = new BrDataTable("DigBB Right"); OutputTable->AddObject(DigitizedBB_Right); BrDataTable *DigitizedBB_Left = new BrDataTable("DigBB Left"); OutputTable->AddObject(DigitizedBB_Left); // // The geant input has 4 tables for the hits in the beam counter // The output is only dicided into left and right //l int LeftId = 0; int RightId = 0; for(option=0; option < 4; option++) { switch(option){ case 0: strcpy(TableName, "GeantHits PURA0"); notubes = 5; Right = kTRUE; firstId = 1; break; case 1: strcpy(TableName,"GeantHits PURB0"); notubes = 30; Right = kTRUE; firstId=6; break; case 2: strcpy(TableName,"GeantHits PULA0"); notubes = 8; firstId=1; Right = kFALSE; break; case 3: strcpy(TableName ,"GeantHits PULB0"); notubes=36; firstId=9; Right = kFALSE; break; default: strcpy(TableName ,"**Wrong option"); } if(DebugLevel() > 4) cout << " * Inspecting " << TableName << endl; GeantHits = (BrDataTable*)InputTable->GetObject(TableName); if(GeantHits == 0 ) { Abort("Event", "no GeantHit Table for %s, Should not happen", TableName); return; } NumHits = GeantHits->Entries(); // int ip; // int tube; int subModule; int* tubetable = new int [notubes]; BrGeantHit** hittable = new BrGeantHit* [NumHits]; for(ip=0;ip<notubes;ip++) tubetable[ip]=0; if(DebugLevel() > 4) cout << " * No of Input hits = " << NumHits << endl; if(NumHits == 0) continue; // // Loop over all tubes; // Check on hits and create BBDig according to the algorithm given // in the class description. // BrGeantHit* hit_p; int nlookup, n; Float_t adc, tdc, length; // t0, Float_t a,b, time; for(int is=0; is< notubes; is++){ int hitno=0; if(DebugLevel()>3){ cout << "Check tube No. " << is << endl; } // // Find all hits for this tube // and fill the pointer table. // for(int ih=0; ih < NumHits; ih++){ ghit_p = (BrGeantHit*) GeantHits->At(ih); // // this kind of check should really be done using assert(ghit_p) // if(ghit_p == NULL){ cerr << "**error ought not to happen" << endl; break; } tube = ghit_p->Isub()-1; if(DebugLevel() > 9){ cout << " ** tubeno " << setw(4) << tube << " de/dx " << setw(8) << ghit_p->Dedx() << endl; } if (tube == is){ if(ghit_p->Dedx() > 0.0){ if(DebugLevel() > 9) cout << " ** Add tubeno " << tube << endl; tubetable[tube]++; hit_p = ghit_p; hittable[hitno++] = hit_p; } } } length = fParamsBB_p->GetUVTLength(); nlookup = tubetable[is]; // // Find the first time for the tube // as for single hits it is required that beta is above threshold // tdc = 9999.0; adc = 0.0; float adcb = 0.0; // primary adc for(n=0;n<nlookup;n++){ hit_p = hittable[n]; // // Should we accept this hit. As for single hit look at // beta etc. float ux; float uy; float uz; ux = (hit_p->LocalXPosOut()-hit_p->LocalXPosIn()); uy = (hit_p->LocalYPosOut()-hit_p->LocalYPosIn()); uz = (hit_p->LocalZPosOut()-hit_p->LocalZPosIn()); BrVector3D Vpath(ux,uy,uz); float path_length = Vpath.Norm(); float path_theta = Vpath.Theta(); if(!Right) path_theta = TMath::Pi() - path_theta; int pid = hit_p->Pid(); Double_t mass= BrGeantHit::Mass(pid); // Get beta from Pid and Pdet // Need to implement GParticle::Mass(Pid) or similar set of simple // functions. // pdet = hit_p->GetPdet(); // mass = GParticle::Mass(hit_p->GetPid()); // beta = pdet / TMath::Sqrt(pdet*pdet + mass*mass); // Double_t p = hit_p->Pdet(); Double_t beta = p/sqrt(p*p+mass*mass); // // Double_t beta_threshold = 0.67; if (DebugLevel() > 14) cout << " ** P: " << setw(8) << Float_t(p) << " MeV/c " << " Pid: " << setw(4) << pid << " Mass: " << setw(8) << Float_t(mass) << " Mev" << " Beta: " << setw(8) << Float_t(beta) << " (" << beta_threshold << ")" << endl; if (beta > beta_threshold) { float Npe; float r = beta_threshold/beta; Npe = fParamsBB_p->GetNpePercm() * (1.0 - r*r) * path_length; if (DebugLevel() > 14){ cout << " ** Pathlength: " << setw(8) << path_length << " " << setw(8) << length << " estimated PE: " << setw(8) << Npe << " T: " << setw(8) << hit_p->Tof() << " Theta path: " << setw(8) << path_theta << endl; } // // correction due to theta implemented. // if (path_theta < fParamsBB_p->GetThetaCut()) Npe = Npe *( 1.0 - path_theta / fParamsBB_p->GetThetaCut() * fParamsBB_p->GetThetaAtt()); else Npe=0.; if (DebugLevel() > 14) cout << "adjusted PE " << Npe << endl; adc += Npe; if(HistOn()){ int subModule; if(Right) subModule = fParamsBB_p->ConvertGeantModuleNo(kTRUE, is+firstId); else subModule = fParamsBB_p->ConvertGeantModuleNo(kFALSE, is+firstId); BrGeantTrack* track = hit_p->GetTrack(); if(track){ if(track->IsParent()) adcb +=Npe; } } time = hit_p->Tof(); if(time < tdc) tdc = time; } } if(HistOn()){ int subModule; if(Right) subModule = fParamsBB_p->ConvertGeantModuleNo(kTRUE, is+firstId); else subModule = fParamsBB_p->ConvertGeantModuleNo(kFALSE, is+firstId); if(Right){ hradc[subModule]->Fill(adc); hpradc[subModule]->Fill(adcb); } else{ hladc[subModule]->Fill(adc); hpladc[subModule]->Fill(adcb); } } // // Add the new BrBbDig to data. // For the time the fastest signal at each end is used. // This code does not simulate any slewing behaviour // For adc signals the total is added up for time's within 100nsec // of the fastest time. // if(DebugLevel() > 9){ cout << "Creating new data member" << is << endl; } if(adc > 0.0){ digbb_p = new BrBbDig(); if(Right){ DigitizedBB_Right->Add(digbb_p); digbb_p->SetId(++RightId); } else{ DigitizedBB_Left->Add(digbb_p); digbb_p->SetId(++LeftId); } if(Right) subModule = fParamsBB_p->ConvertGeantModuleNo(kTRUE, is+firstId); else subModule = fParamsBB_p->ConvertGeantModuleNo(kFALSE, is+firstId); digbb_p->SetTubeNo(subModule); gRandom->Rannor(a,b); tdc += a* fParamsBB_p->GetSigmaTime(); tdc = tdc/fParamsBB_p->GetTdcConv() +fParamsBB_p->GetTdcOffset(); Int_t itdc = (Int_t)tdc; digbb_p->SetTdc(itdc); Float_t Npe = adc; adc *= fParamsBB_p->GetAdcGain(); adc += fParamsBB_p->GetAdcOffset(); Float_t sigma; Bool_t BigTube = fParamsBB_p->IsModuleBigTube(Right, subModule); if(BigTube) { float eta = fParamsBB_p->GetQEBigTube(); float delta = fParamsBB_p->GetEmissionGainBigTube(); sigma = TMath::Sqrt((1.0-eta+1.0/(delta-1))/(Npe*eta)+1.0/Npe); } else { float eta = fParamsBB_p->GetQESmallTube(); float delta = fParamsBB_p->GetEmissionGainSmallTube(); sigma = TMath::Sqrt((1.0-eta+1.0/(delta-1))/(Npe*eta)+1.0/Npe); } adc *= (1.0+ b* sigma); Int_t iadc = (Int_t)adc; digbb_p->SetAdc(iadc); if(DebugLevel()>2){ digbb_p->List(); } } // adc> 0 } // is Numhits delete [] tubetable; } // option } //______________________________________________________________________ void BrDigitizeBB::ListDetectorParameters() const { // List the current value of the digitization parameters on // standard out. if( fParamsBB_p != NULL){ cout << "The name is " << GetName() <<"n"; fParamsBB_p->ListParameters(); } else cout << "No parameters set for this Digitisation module" << endl; } //____________________________________________________________________________ void BrDigitizeBB::Print(Option_t* option) const { // // Standard info printout. This should give the nme of the module, // the author, and most recent revision. This should be implemented // extracting some of the info from the cvs -id etc. // TString opt(option); opt.ToLower(); BrModule::Print(option); if (opt.Contains("d")) cout << endl << " Original author: Flemming Videbaek" << endl << " Revisted by: $Author: cholm $" << endl << " Revision date: $Date: 2001/08/12 13:28:15 $" << endl << " Revision number: $Revision: 1.3 $ " << endl << endl << "*************************************************" << endl; } ///////////////////////////////////////////////////////////////////////////////// // // $Log: BrDigitizeBB.cxx,v $ // Revision 1.3 2001/08/12 13:28:15 cholm // Compactified debug output, and increased most debuging levels. // // Revision 1.2 2001/06/22 17:40:00 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.26 2001/03/07 12:16:45 cholm // * Made the method BrModule::Info() const obsolete in favour of // BrModule::Print(Option_t* option="B") const. // // Revision 1.25 2001/01/17 01:58:35 hagel // Fixed float to int warning that made into error // // Revision 1.24 2000/12/16 01:04:22 videbaek // cosmetic changes for Digitizew // The BB rdo assigned a vertex valu of 0even if not found properly. // Constructor now sets to -999. to avoid confusion. // // Revision 1.23 2000/11/23 01:31:09 brahmlib // Cleaned up code in general. // // Revision 1.22 2000/11/22 21:14:45 videbaek // Moved stuff from constructor to Init() // // Revision 1.21 2000/11/17 22:00:38 videbaek // Move the detector param opening from constructor to Init(). // // Revision 1.20 2000/11/10 19:29:18 videbaek // More histograms to digitization // // Revision 1.19 2000/09/23 19:16:02 videbaek // Improved histograms for looking in details at digitization. // // Revision 1.18 2000/09/14 19:30:48 videbaek // Added parameters for better tube resolution estimate. // Fixed erros in module counting. // Moved to single loop over tubes rather than one. // // Revision 1.17 2000/07/16 19:52:58 videbaek // Remove volume information form Digitize BB // // Revision 1.16 2000/01/06 16:03:21 videbaek // Cosmetic changes // // Revision 1.15 1999/11/30 22:21:19 videbaek // Added method to convert geant module no to real life module no of // beam beam counters. // // Revision 1.14 1999/04/14 20:12:40 videbaek // Use DbManager methods to access parameters. // // Revision 1.13 1999/04/01 22:22:07 videbaek // Changed to use BrParameterDbManager for parameters. // // Revision 1.12 1999/01/21 23:35:42 hagel // Changed #ifndef's in .cxx files for include files to reflect changes to // .h files according to BRAT convention. It should be noted that very few // of the .cxx source files have the protections. Should this be systematically // changed. // // Revision 1.11 1998/12/21 20:28:09 videbaek // changed data constants // // Revision 1.10 1998/09/17 16:16:31 videbaek // Removed Mass() static fct from BrDigitizeBB // // Revision 1.9 1998/09/10 21:25:40 videbaek // Remove references to TSonataMath // // Revision 1.8 1998/08/28 20:04:12 videbaek // remove a simple get UVT length in list // // Revision 1.7 1998/08/21 13:33:54 videbaek // improvements to algorithm. Consider theta // // Revision 1.6 1998/08/04 20:24:26 videbaek // Changed detector parameters // // Revision 1.5 1998/08/03 20:54:08 videbaek // Fix algorithms for bb // // Revision 1.4 1998/07/31 19:32:16 videbaes // Working version of BB digitize // // Revision 1.3 1998/07/28 21:29:13 videbaek // several fixes to BB // // Revision 1.2 1998/07/27 16:27:05 videbaek // Stat on modifiing BB - not complete // // Revision 1.1 1998/07/02 19:40:36 videbaek // |
||||||
|
This page automatically generated by script docBrat by Christian Holm |
Copyright ; 2002 BRAHMS Collaboration
<brahmlib@rcf.rhic.bnl.gov>
|