|
//____________________________________________________________________ // // This module is devoted to the DC detector calibration //____________________________________________________________________ // // $Id: BrDcCalModule.cxx,v 1.17 2002/05/07 17:35:24 ouerdane Exp $ // $Author: ouerdane $ // $Date: 2002/05/07 17:35:24 $ // $Copyright: (C) 2001 BRAHMS Collaboration <brahmlib@rhic.bnl.gov> // /////////////////////////////////////////////////////////////////////// // // // BrDcCalModule // // // // // // // // BrDcCalModule is a calibration module for the // // brahms drift chambers // // // // It also will store data in the Calibration Database using the // // BrCalibrationManager // // // // // /////////////////////////////////////////////////////////////////////// #ifndef WIN32 #include <iostream> #else #include <iostream.h> #endif #include <iomanip.h> #include <fstream.h> // ROOT Includes #ifndef ROOT_TF1 #include "TF1.h" #endif #ifndef ROOT_TString #include "TString.h" #endif #ifndef ROOT_TMath #include "TMath.h" #endif #ifndef ROOT_TDirectory #include "TDirectory.h" #endif // Brat stuff #ifndef BRAT_BrTofPedCalModule #include "BrDcCalModule.h" #endif #ifndef BRAT_BrDataTable #include "BrDataTable.h" #endif #ifndef BRAT_BrDcDig #include "BrDcDig.h" #endif #ifndef BRAT_BrDCHit #include "BrDCHit.h" #endif #ifndef BRAT_BrParameterDbManager #include "BrParameterDbManager.h" #endif #ifndef BRAT_BrCalibrationManager #include "BrCalibrationManager.h" #endif #ifndef BRAT_BrRunInfoManager #include "BrRunInfoManager.h" #endif #ifndef BRAT_BrRunInfo #include "BrRunInfo.h" #endif #ifndef BRAT_BrBbVertex #include "BrBbVertex.h" #endif //________________________________________________________________________ ClassImp(BrDcCalModule); //________________________________________________________________________ BrDcCalModule::BrDcCalModule() : BrModule() { // // default constructor // fCommitCal = kFALSE; fCalibration = 0; fGlobalTdcOffset = BrDcCalibration::kCalException; fDriftTime = BrDcCalibration::kCalException; fCalComment[0] = ""; fCalComment[1] = ""; fLocalTrack = 0; } //________________________________________________________________________ BrDcCalModule::BrDcCalModule(Text_t *Name, Char_t *Title) : BrModule(Name, Title) { // // constructor for named object // fEvents=0; fCommitCal = kFALSE; fCalibration = 0; fGlobalTdcOffset = BrDcCalibration::kCalException; fDriftTime = BrDcCalibration::kCalException; fCalComment[0] = ""; fCalComment[1] = ""; fLocalTrack = 0; } //________________________________________________________________________ BrDcCalModule::~BrDcCalModule() { // // destructor // } //_______________________________________________________________________ void BrDcCalModule::Init() { // // Called once per session // Get detectorParameters. Inform the parameterElementManeger about // database tables to be used (i.e. filled). // DC description initialization // Finally it should go to BrDetectorParamsDC if (fCommitCal) { BrCalibrationManager* calMan = BrCalibrationManager::Instance(); fCalibration = (BrDcCalibration*)calMan->Register("BrDcCalibration", GetName()); if (!fCalibration) { Abort("Init", "couldn't get pointer to DC calibration parameter"); return; } if (fGlobalTdcOffset != BrDcCalibration::kCalException) fCalibration->Use("tdcOffset", BrCalibration::kWrite, 1); if (fDriftTime != BrDcCalibration::kCalException) fCalibration->Use("driftTime", BrCalibration::kWrite, 1); return; } fLocalTrack = new BrDCTrackingModule((Char_t*)GetName(),(Char_t*)GetName(),kFALSE); fLocalTrack -> UseMySql(); fLocalTrack->Init(); fPlaneView[0]='X'; fPlaneView[1]='Y'; fPlaneView[2]='U'; fPlaneView[3]='V'; if(!strcmp("T3",GetName())) { // case of T3 fNumPlanes=10; fWires[0]=40; fWires[1]=40; fWires[2]=40; // X view fWires[3]=30; fWires[4]=30; fWires[5]=30; // Y view fWires[6]=48; fWires[7]=48; // U view fWires[8]=48; fWires[9]=48; // V view fViewPlanes[0]=3;fViewPlanes[1]=3;fViewPlanes[2]=2;fViewPlanes[3]=2; fView[0]='X'; fView[1]='X'; fView[2]='X'; fView[3]='Y'; fView[4]='Y'; fView[5]='Y'; fView[6]='U'; fView[7]='U'; fView[8]='V'; fView[9]='V'; fViewNumber[0]=1; fViewNumber[1]=2; fViewNumber[2]=3; fViewNumber[3]=1; fViewNumber[4]=2; fViewNumber[5]=3; fViewNumber[6]=1; fViewNumber[7]=2; fViewNumber[8]=1; fViewNumber[9]=2; fRefDriftDist=0.5; cout<<"was in T3 section"<<endl; }else if(!strcmp("T4",GetName()) || !strcmp("T5",GetName())){ fNumPlanes=8; fWires[0]=23; fWires[1]=23; // X view fWires[2]=16; fWires[3]=16; // Y view fWires[4]=27; fWires[5]=27; // U view fWires[6]=27; fWires[7]=27; // V view fViewPlanes[0]=2;fViewPlanes[1]=2;fViewPlanes[2]=2;fViewPlanes[3]=2; fView[0]='X'; fView[1]='X'; fView[2]='Y'; fView[3]='Y'; fView[4]='U'; fView[5]='U'; fView[6]='V'; fView[7]='V'; fViewNumber[0]=1; fViewNumber[1]=2; fViewNumber[2]=1; fViewNumber[3]=2; fViewNumber[4]=1; fViewNumber[5]=2; fViewNumber[6]=1; fViewNumber[7]=2; fRefDriftDist=1.1; cout<<"was in T4/5 section"<<endl; }else{ cout<<"Wrong detector name: '"<<GetName()<<"'"<<endl; cout<<"should be 'T3', 'T4', or 'T5'"<<endl; } // driver settings if(!fDCDriver)InitDCDriver(); if(!strcmp("T4",GetName()) || !strcmp("T5",GetName())){ if(!ReadPromptDelayMap())cout<<"Prompt-Delay mapping file doesn't exist"<<endl; } // Job-level initialization SetState(kInit); BrModule::Init(); // for the moment we don't use it - if you do use this // form!! //fCalibration->Use("topPedestal", BrCalibration::kWrite,fParamsTof->GetNoSlats()); // NOT THIS. //fCalibration->UseForWrite("topPedestalWidth", fParamsTof->GetNoSlats()); //fCalibration->UseForWrite("botPedestal", fParamsTof->GetNoSlats()); //fCalibration->UseForWrite("botPedestalWidth", fParamsTof->GetNoSlats()); cout<<"Initialazing BadWire matrix"<<endl; for(Int_t im=0;im<3;im++){ for(Int_t ip=0;ip<10;ip++){ for(Int_t iw=0;iw<48;iw++){ fBadWire[im][ip][iw]=kFALSE; } } } cout<<"BrDcCalModule from brat2"<<endl; } //______________________________________________________________ void BrDcCalModule::InitDCDriver() { // driver settings if(!fDCDriver){ fDCDriver = new BrDriverDC(GetName(),GetName()); if(!strcmp("T3",GetName())){ fDCDriver->SetMinWidth(fMinWidth); fDCDriver->SetMinTdc(fMinTdc); fDCDriver->SetMaxTdc(fMaxTdc); fCalibChan=290; } if(!strcmp("T4",GetName()) || !strcmp("T5",GetName()) ){ fDCDriver->SetMinWidth(fMinWidth); fDCDriver->SetMinTdc(fMinTdc); fDCDriver->SetMaxTdc(fMaxTdc); fCalibChan=460; } cout<<"Cleaning settings for "<<GetName()<<" detector: " <<fMinTdc<<" "<<fMaxTdc<<" "<<fMinWidth<<endl; }else{ cout<<"fDCDriver already exist->??????"<<endl; } } //________________________________________________________________________ void BrDcCalModule::DefineHistograms() { //----------------------------------------------------------------------------// // one dim plots cout<<"booking histograms for "<<GetName()<<" detector"<<endl; Char_t HistName[64]; for(Int_t imod=0;imod<3;imod++) { // imod denotes the module number for(Int_t i=0;i<fNumPlanes;i++) { // i denotes the plane number //cout<<"imod: "<<imod<<" ip: "<<i<<" fNum: "<<fViewNumber[i]<<endl; sprintf(HistName,"%s_%d_%c%d",GetName(),imod+1,fView[i],fViewNumber[i]); fPlane[imod][i] = new TH1F(HistName,HistName,fWires[i],.5,fWires[i]+0.5); fPlane[imod][i]->SetYTitle("counts"); fPlane[imod][i]->SetXTitle("wire number"); sprintf(HistName,"%s_%d_%c%d_Time",GetName(),imod+1,fView[i],fViewNumber[i]); fTime[imod][i] = new TH1F(HistName,HistName,2000,0,6000); fTime[imod][i]->SetYTitle("dN/dt"); fTime[imod][i]->SetXTitle("time [ns]"); //sprintf(HistName,"%s_%d_%c%d_TimeVsWidth",GetName(),imod+1,fView[i],fViewNumber[i]); //fTimeVsWidth[imod][i] = new TH2F(HistName,HistName,500,4000,5500,260,0,260); } } sprintf(HistName,"%sAssoHits",GetName()); fhAssoHits = new TH1F(HistName,"Associated hits for all rec. tracks",41,-.5,40.5); fhAssoHits->SetYTitle("counts"); fhAssoHits->SetXTitle("hits number"); sprintf(HistName,"%sAssoHitsAcce",GetName()); fhAssoHitsAcce = new TH1F(HistName,"Associated hits for accepted rec. tracks",41,-.5,40.5); fhAssoHitsAcce->SetYTitle("counts"); fhAssoHitsAcce->SetXTitle("hits number"); sprintf(HistName,"%sHitsPerPlane",GetName()); fhHitsPerPlane = new TH1F(HistName,"Hits Per Plane",30,0.5,30.5); fhHitsPerPlane->SetYTitle("counts"); fhHitsPerPlane->SetXTitle("Plane Number"); sprintf(HistName,"%sAllTracks",GetName()); fhAllTracks = new TH1F(HistName,"tracks",1,0.5,1.5); } //________________________________________________________________________ void BrDcCalModule::Begin() { // Run-level initialization SetState(kBegin); if (fCommitCal) return; if(!fLocalTrack)Failure("Begin", "can not retrive parameters from BrDCTrackingModule"); fLocalTrack->Begin(); fDCDriver->SetMinTdc((fLocalTrack->GetClusterFinder())->GetMinTdc()); fDCDriver->SetMaxTdc((fLocalTrack->GetClusterFinder())->GetMaxTdc()); fPromptDelay = (fLocalTrack->GetClusterFinder())->GetPromptDelay(); fMinTdc = (fLocalTrack->GetClusterFinder())->GetMinTdc(); fMaxTdc = (fLocalTrack->GetClusterFinder())->GetMaxTdc(); if (fReadTdcMatrix) return; // check if histos are booked if (!HistBooked()) Failure("Begin", "must book histograms for calibration to work"); } //________________________________________________________________________ void BrDcCalModule::Event(BrEventNode* inNode, BrEventNode* outNode) { // // Fill the histograms that hold the Pedestal spectra // for synchronization events // // Required input e.g. : BrDataTable "DigTOF TOF1" // // SetState(kEvent); if (fCommitCal) return; if(DebugLevel() > 1) cout << "Entering Event() in BrDcCalModule for " << GetName() << endl; if (fReadTdcMatrix) return; //inNode->ListObjects(); if(fCalibLevel)DecodeRedData(inNode); else DecodeRawData(inNode); // get BB stuff BrBbVertex* bb = (BrBbVertex*)inNode->GetObject("BB Vertex"); if (!bb) { bb = (BrBbVertex*)outNode->GetObject("BB Vertex"); if (!bb) { if (DebugLevel() > 1) cout << " No beam beam stuff" << endl; } } if(!bb && fRunType=="AuAu")return; if(bb && fRunType=="AuAu"){ if (bb->GetMethod() != 2) // Small tubes only return; fT0 = bb->GetTime0(); fZ0 = bb->GetZ0(); fRightTime = bb->GetRightTime(); fLeftTime = bb->GetLeftTime(); //cout<<"fT0,fZ0: "<<fT0<<" "<<fZ0<<endl; } //cout<<"t1: "<<fNumCleanHits<<endl; fEvents++; if(fEvents%5000==0)cout<<"Events --> "<<fEvents<<" "<<flush; // for each event fill Tdc Matrix //if(fNumCleanHits)cout<<"entering FillTdcMatrix(): "<<fNumCleanHits<<endl; FillTdcMatrix(); } //________________________________________________________________________ void BrDcCalModule::Finish() { // // called once per job // // Job-level finalization SetState(kFinish); if (fCommitCal) { if (fGlobalTdcOffset != BrDcCalibration::kCalException) { Info(5, "Finish", "Committing tdc offset (%d) for %s with this comment:n %s", fGlobalTdcOffset, GetName(), fCalComment[0].Data()); fCalibration->SetTdcOffset(fGlobalTdcOffset); fCalibration->SetComment("tdcOffset", fCalComment[0].Data()); } if (fDriftTime != BrDcCalibration::kCalException) { Info(5, "Finish", "Committing drift time (%d) for %s with this comment:n %s", fDriftTime, GetName(), fCalComment[1].Data()); fCalibration->SetDriftTime(fDriftTime); fCalibration->SetComment("driftTime", fCalComment[1].Data()); } return; } if(!fReadTdcMatrix)SaveTdcMatrix(); if(!strcmp("T4",GetName()) || !strcmp("T5",GetName())){ if(!fMapExist)SavePromptDelayMap(); } if(!fReadTdcMatrix){ SaveTdcMatrix(); //return; // if tdc matrix was created in this run then it is enough to write it down. } FillWithTdcMatrix(1,3); // special treatment for this plane if(fOffsetOption){ PrepareForOffsets(); FindTdcOffsets1New(); // summing over all wires in the plane - require PrepareForOffsets(). }else{ FindTdcOffsetsNew(); } //if(!strcmp("T4",GetName()))CorrectOffsets(); //FindTdcOffsets(); // each wire treated independly //FindTdcOffsetsNew(); SaveTdcOffsets(); Int_t CreatHist=1; ConstructTimeMatrix(CreatHist); Calibrate(); SaveCalibData(); // set parameters for eventual commit into DB //fCalibration->SetBotPedestal(slat, meanFit); //fCalibration->SetBotPedestalWidth(slat, sigmaFit); //fCalibration->SetComment("topPedestal","Generated by //BrTofPedCalModule: fit with a Gaussian function"); } //____________________________________________________________________ Bool_t BrDcCalModule::ReadTdcMatrix() { // save pedestal to ascii file // but first check if fDCDriver has already been defined if(!fDCDriver)InitDCDriver(); fReadTdcMatrix=kFALSE; if (fTdcMatrixFile == "") { cout << "n You forgot to set a tdc matrix file for " << GetName() << ".n Please provide it now: " << flush; cin >> fTdcMatrixFile; cout << endl; cout<<"Good Boy!"<<endl; } Int_t im,ip,iw,it,r; ifstream IN(fTdcMatrixFile.Data(), ios::in); if(IN){ fReadTdcMatrix=kTRUE; cout<<"Reading TDC Matrix"<<endl; do { IN>>im>>ip>>iw>>it>>r; if(!(IN.eof())){ //if(fDCDriver->HitIsOK(it,fMinWidth)){ if(im==0 && ip==8 && iw>6 && iw<22)cout<<"im,ip,iw "<<im<<" "<<ip<<" "<<iw<<" not used"<<endl; else fTdcMatrix[im][ip][iw][it]+=r; //} } } while(!(IN.eof())); IN.close(); } return fReadTdcMatrix; } //_____________________________________________________ Bool_t BrDcCalModule::ReadPromptDelayMap() { // maping to ascii file fMapExist=kFALSE; Int_t im,ip,iw,id; Char_t MapFileName[64]; sprintf(MapFileName,"%sPromptDelayMap.dat",GetName()); ifstream IN(MapFileName, ios::in); if(IN){ fMapExist=kTRUE; cout<<"Reading Prompt-Delay Mapping from file: "<<MapFileName<<endl; do { IN>>im>>ip>>iw>>id; if(!(IN.eof())){ fDelayed[im][ip][iw]=id; } } while(!(IN.eof())); IN.close(); } return fMapExist; } //_____________________________________________________________________ void BrDcCalModule::FillTdcMatrix() { // This method should be called in the event loop // to fill tdc spectrum with each new event Int_t NumHits=fNumCleanHits; for(Int_t i=0;i<NumHits;i++) { Int_t im = fEveTab[i][0]; Int_t ip = fEveTab[i][1]; Int_t iw = fEveTab[i][2]; Int_t it = fEveTab[i][3]; //Int_t itBB = GetTdcFromBBRef((Double_t)it); //cout<<"itBB: "<<itBB<<endl; //if(itBB>=0 && itBB<5200)fTdcMatrix[im][ip][iw][itBB]++; fTdcMatrix[im][ip][iw][it]++; } } //_____________________________________________________________________ Int_t BrDcCalModule::GetTdcFromBBRef(Double_t it) { // it is in units of 1/2 ns //Int_t TimeOffset=109; Int_t TimeOffset=123; //Double_t TdcReference = (fT0) * 2.0; // TdcReferece is now in units of 500 ps //Double_t TdcReference = (fT0-fZ0/30.0) * 2.0; Double_t TdcReference = fLeftTime*2.0; //TdcReference=0.0; //cout<<TdcReference<<" "<<fZ0<<endl; return TMath::Nint(it + TdcReference) - TimeOffset; } //_____________________________________________________________________ void BrDcCalModule::FillTdcMatrixView() { // This method should be called in the event loop // to fill time spectrum with for each new event if(strcmp("T3",GetName())){ cout<<"This method can not be call for " <<GetName()<<" detector"<<endl; } Int_t NumHits=fNumCleanHits; for(Int_t i=0;i<NumHits;i++) { Int_t im = fEveTab[i][0]; Int_t ip = fEveTab[i][1]; Int_t iw = fEveTab[i][2]; Int_t it = fEveTab[i][3]; // comment on what follow: we multiply filling of fTdcMatrix // to gain the statistics. It is of course temporary. if(fView[ip]=='X'){ for(Int_t j=0;j<3;j++){ //cout<<ip<<" "<<it<<endl; fTdcMatrix[im][j][iw][it]++; //X view } } if(fView[ip]=='Y'){ for(Int_t j=3;j<6;j++){ fTdcMatrix[im][j][iw][it]++; //Y view } } if(fView[ip]=='U'){ for(Int_t j=6;j<8;j++){ fTdcMatrix[im][j][iw][it]++; //U view } } if(fView[ip]=='V'){ for(Int_t j=8;j<10;j++){ fTdcMatrix[im][j][iw][it]++; //V view } } } } //_____________________________________________________________________ void BrDcCalModule::SaveTdcMatrix() { // This method write TDC matrix to the disk ofstream OUT(fTdcMatrixFile.Data(), ios::out); OUT.setf(ios::showpoint); cout<<"Writing TDC Matrix"<<endl; for(Int_t im=0;im<3;im++) { for(Int_t ip=0;ip<fNumPlanes;ip++) { fDCDriver->AtPlane(ip); for(Int_t iw=0;iw<fDCDriver->GetMaxWireNum();iw++) { for(Int_t it=0;it<5200;it++) { if(fTdcMatrix[im][ip][iw][it])OUT<<im<<" "<<ip<<" "<<iw<<" "<<it<<" " <<fTdcMatrix[im][ip][iw][it]<<endl; } //it } //iw } //ip } //im OUT.close(); } //_____________________________________________________________________ void BrDcCalModule::SavePromptDelayMap() { // Char_t MapFileName[64]; sprintf(MapFileName,"%sPromptDelayMap.dat",GetName()); ofstream OUT(MapFileName, ios::out); cout<<"Writing Pprompt-Delay Mapping"<<endl; for(Int_t im=0;im<3;im++) { for(Int_t ip=0;ip<fNumPlanes;ip++) { fDCDriver->AtPlane(ip); for(Int_t iw=0;iw<fDCDriver->GetMaxWireNum();iw++) { OUT<<im<<" "<<ip<<" "<<iw<<" "<<fDelayed[im][ip][iw]<<endl; } //iw } //ip } //im OUT.close(); } //_____________________________________________________________________ void BrDcCalModule::FillWithTdcMatrix(Int_t module,Int_t plane) { cout<<"Filling for TDC spectra"<<endl; Char_t HistName[64]; fDCDriver->AtPlane(plane); for(Int_t iw=0;iw<fDCDriver->GetMaxWireNum();iw++) { sprintf(HistName,"%sTdcFor%d_%d_w%d",GetName(),module,plane,iw+1); fBTime[iw] = new TH1F(HistName,HistName,2600,2600,5200); // to monitor time fBTime[iw]->SetYTitle("dN/dt");fBTime[iw]->SetXTitle("TDC Channel"); } for(Int_t im=0;im<3;im++) { for(Int_t ip=0;ip<fNumPlanes;ip++) { fDCDriver->AtPlane(ip); for(Int_t iw=0;iw<fDCDriver->GetMaxWireNum();iw++) { for(Int_t it=fMinTdc;it<fMaxTdc;it++) { fTime[im][ip]->Fill(it,fTdcMatrix[im][ip][iw][it]); fPlane[im][ip]->Fill(iw+1,fTdcMatrix[im][ip][iw][it]); if(im==module && ip==plane)fBTime[iw]->Fill(it,fTdcMatrix[im][ip][iw][it]); } //it } //iw } //ip } //im } //_____________________________________________________________________ void BrDcCalModule::PrepareForOffsets() { // If you want to perform individual calibration for a single // plane, first sum fTdcMatrix over wires. This method should be // called befor calling FindTdcOffsets1(). cout<<"Summing TdcMatrix over wires "<<endl; for(Int_t im=0;im<3;im++){ for(Int_t ip=0;ip<fNumPlanes;ip++){ fDCDriver->AtPlane(ip); for(Int_t iw=0;iw<fDCDriver->GetMaxWireNum();iw++){ Int_t id=0; if(!strcmp("T4",GetName()) || !strcmp("T5",GetName()))id=fDelayed[im][ip][iw]; for(Int_t it=3000;it<6000;it++){ fForOffsets[im][ip][id][it] += fTdcMatrix[im][ip][iw][it]; } } } } } //_____________________________________________________________ void BrDcCalModule::FindTdcOffsets() { // This method looks for tdc channel number that refers to // zero drift distance (time). Here we treated each channel (wire) independly. cout<<"Looking for TDC offsets"<<endl; for(Int_t im=0;im<3;im++) { for(Int_t ip=0;ip<fNumPlanes;ip++) { fDCDriver->AtPlane(ip); for(Int_t iw=0;iw<fDCDriver->GetMaxWireNum();iw++) { Int_t it=5199; Int_t nstep=0; Int_t sum=0; Float_t mean,meanLow=0.0; Int_t next; Int_t current; Int_t prev0,prev1,prev2; do { nstep++; sum += fTdcMatrix[im][ip][iw][it]; mean = (Float_t)sum/(Float_t)nstep; if(it>4960)meanLow=mean; if(it<4625 && it>4450)meanLow=mean; if(it==4570)mean=0; next=fTdcMatrix[im][ip][iw][it-10]; it--; } while((mean<meanLow*3.0 || mean==0 || it>5000 ) && it>fMinTdc); do { prev2=fTdcMatrix[im][ip][iw][it-3]; prev1=fTdcMatrix[im][ip][iw][it-2]; prev0=fTdcMatrix[im][ip][iw][it-1]; current=fTdcMatrix[im][ip][iw][it]; it++; next=fTdcMatrix[im][ip][iw][it]; } while((0.2*(prev2+prev1+prev0+current+next) > mean*2.3) && it+1<fMaxTdc); fOffset[im][ip][iw]=it; // Is it realy enough? //if(im==2 && ip==6)cout<<" "<<it<<endl; //if(im==2 && ip==6){Int_t ii; cin>>ii;} } //iw } //ip } //im } //_____________________________________________________________ void BrDcCalModule::FindTdcOffsetsNew() { // This method looks for tdc channel number that refers to // zero drift distance (time). Each wire is treated independly. Bool_t WireHasProblem[48]; Float_t Param=.5; Float_t param; cout<<"FindTdcOffsetsNew::Looking for TDC offsets"<<endl; for(Int_t im=0;im<3;im++) { for(Int_t ip=0;ip<fNumPlanes;ip++) { cout<<"Looking for offsets for plane "<<ip<<" in module "<<im<<endl; //param=Param; fDCDriver->AtPlane(ip); for(Int_t iw=0;iw<fDCDriver->GetMaxWireNum();iw++) { WireHasProblem[iw]=kFALSE; param=Param; Int_t id = 0; Int_t MaxTdc = fMaxTdc; Int_t MinTdc = fMinTdc; if(!strcmp("T4",GetName()) || !strcmp("T5",GetName())){ id=fDelayed[im][ip][iw]; //if(iw && id!=fDelayed[im][ip][iw-1])param=Param; if(id){ MaxTdc=fPromptDelay; MinTdc=fMinTdc; }else{ MaxTdc=fMaxTdc; MinTdc=fPromptDelay; } } Int_t it; Int_t nstep=0; Int_t sum=0; Float_t mean,meanLow; Int_t next; Int_t current; Int_t prev0,prev1,prev2; Bool_t Continue; do { it=MaxTdc; meanLow=0.0; do { nstep++; sum += fTdcMatrix[im][ip][iw][it]; mean = (Float_t)sum/(Float_t)nstep; next=fTdcMatrix[im][ip][iw][it-1]; Int_t nextP3=fTdcMatrix[im][ip][iw][it-3]; Int_t nextP4=fTdcMatrix[im][ip][iw][it-4]; Int_t nextP5=fTdcMatrix[im][ip][iw][it-5]; Int_t nextP6=fTdcMatrix[im][ip][iw][it-6]; Continue = (next<=sum*param || nextP3<=sum*param || nextP4<=sum*param || nextP5<=sum*param || nextP6<=sum*param); it--; } while((Continue || mean==0 || it>MaxTdc-40) && it>MinTdc); if(MinTdc+10>it)param-=0.01; //cout<<param<<endl; }while (MinTdc+10>it && param>0); if(param<0){ cout<<"FindOffsets: procedure has problem for (im,ip,iw) "<<im<<" "<<ip<<" "<<iw<<" "<<param<<endl; WireHasProblem[iw]=kTRUE; } fOffset[im][ip][iw]=it+10; // Is it realy enough? //cout<<"im,ip,iw,offset "<<im<<" "<<ip<<" "<<iw<<" "<<it<<" "<<id<<endl; //if(im==2 && ip==6)cout<<" "<<it<<endl; //if(im==2 && ip==6){Int_t ii; cin>>ii;} } //iw if(Param!=param)cout<<"FindOffsets: Offset finding param corrected(im,ip): " <<im<<" "<<ip<<" "<<Param<<" "<<param<<endl; //cout<<WireHasProblem<<endl; CheckAllWiresInPlane(im,ip,WireHasProblem); } //ip } //im } //_____________________________________________________________ void BrDcCalModule::CheckAllWiresInPlane(Int_t im,Int_t ip,Bool_t* wsk) { for(Int_t iw=0;iw<fDCDriver->GetMaxWireNum();iw++){ if(*wsk){ fBadWire[im][ip][iw]=kTRUE; //take the closest offsets Bool_t OffsetCorrected=kFALSE; Int_t i=1; if(iw<fDCDriver->GetMaxWireNum()/2){ do{ if(!*(wsk+i)){ fOffset[im][ip][iw]=fOffset[im][ip][iw+i]; OffsetCorrected=kTRUE; cout<<"Offset for wire "<<iw<<" in plane "<<ip <<" and module "<<im<<" has been corrected ("<<iw+i<<")"<<endl; } i++; }while(!OffsetCorrected); }else{ do{ if(!*(wsk-i)){ fOffset[im][ip][iw]=fOffset[im][ip][iw-i]; OffsetCorrected=kTRUE; cout<<"Offset for wire "<<iw<<" in plane "<<ip <<" and module "<<im<<" has been corrected ("<<iw-i<<")"<<endl; } i++; }while(!OffsetCorrected); } } wsk++; } } //_____________________________________________________________ void BrDcCalModule::FindTdcOffsets1() { // This method looks for tdc channel number that refers to // zero drift distance (time). Here we look for offsets ofter summing // over all wires in particular plane. This approximation is useful // in the case of low statistics cout<<"Looking for TDC offsets (1)"<<endl; for(Int_t im=0;im<3;im++) { for(Int_t ip=0;ip<fNumPlanes;ip++) { fDCDriver->AtPlane(ip); for(Int_t iw=0;iw<fDCDriver->GetMaxWireNum();iw++) { Int_t id = 0; Int_t MaxTdc = fMaxTdc; Int_t MinTdc = fMinTdc; if(!strcmp("T4",GetName()) || !strcmp("T5",GetName())){ id=fDelayed[im][ip][iw]; if(id){ MaxTdc=fPromptDelay; MinTdc=fMinTdc; }else{ Int_t MaxTdc = fMaxTdc; MaxTdc=fMaxTdc; MinTdc=fPromptDelay; } } Int_t it=MaxTdc; Int_t nstep=0; Int_t sum=0; Float_t mean,meanLow=0.0; Int_t next; Int_t current; Int_t prev0,prev1,prev2; do { nstep++; sum += fForOffsets[im][ip][id][it]; mean = (Float_t)sum/(Float_t)nstep; //if(it>4960)meanLow=mean; //if(it<4625 && it>4450)meanLow=mean; //if(it==4570)mean=0; if(it>MaxTdc-70)meanLow=mean; next=fForOffsets[im][ip][id][it-10]; it--; //} while((mean<meanLow*3.0 || mean==0 || it>5000 ) && it>fMinTdc); } while((mean<meanLow*3.0 || mean==0 || it>fMaxTdc-50) && it>MinTdc); //if(im==2 && ip==6)cout<<ip<<" "<<iw<<" "<<it<<endl; do { prev2=fForOffsets[im][ip][id][it-3]; prev1=fForOffsets[im][ip][id][it-2]; prev0=fForOffsets[im][ip][id][it-1]; current=fForOffsets[im][ip][id][it]; it++; next=fForOffsets[im][ip][id][it]; //} while((0.2*(prev2+prev1+prev0+current+next) > mean*2.3) && it+1<MaxTdc); } while((0.2*(prev2+prev1+prev0+current+next) > mean*1000.) && it+1<MaxTdc); fOffset[im][ip][iw]=it; // Is it realy enough? //cout<<"im,ip,iw,offset "<<im<<" "<<ip<<" "<<iw<<" "<<it<<" "<<id<<endl; //if(im==2 && ip==6)cout<<" "<<it<<endl; //if(im==2 && ip==6){Int_t ii; cin>>ii;} } //iw } //ip } //im } //_____________________________________________________________ void BrDcCalModule::FindTdcOffsets1New() { // This method looks for tdc channel number that refers to // zero drift distance (time). Here we look for offsets ofter summing // over all wires in particular plane. This approximation is useful // in the case of low statistics Float_t Param=.5; Float_t param; cout<<"Looking for TDC offsets (1)"<<endl; for(Int_t im=0;im<3;im++) { for(Int_t ip=0;ip<fNumPlanes;ip++) { param=Param; fDCDriver->AtPlane(ip); for(Int_t iw=0;iw<fDCDriver->GetMaxWireNum();iw++) { Int_t id = 0; Int_t MaxTdc = fMaxTdc; Int_t MinTdc = fMinTdc; if(!strcmp("T4",GetName()) || !strcmp("T5",GetName())){ id=fDelayed[im][ip][iw]; if(iw && id!=fDelayed[im][ip][iw-1])param=Param; if(id){ MaxTdc=fPromptDelay; MinTdc=fMinTdc; }else{ MaxTdc=fMaxTdc; MinTdc=fPromptDelay; } } Int_t it; Int_t nstep=0; Int_t sum=0; Float_t mean,meanLow; Int_t next; Int_t current; Int_t prev0,prev1,prev2; do { it=MaxTdc; meanLow=0.0; do { nstep++; sum += fForOffsets[im][ip][id][it]; mean = (Float_t)sum/(Float_t)nstep; //if(it>4960)meanLow=mean; //if(it<4625 && it>4450)meanLow=mean; //if(it==4570)mean=0; if(it>MaxTdc-70)meanLow=mean; next=fForOffsets[im][ip][id][it-1]; it--; //} while((mean<meanLow*3.0 || mean==0 || it>5000 ) && it>fMinTdc); } while((next<sum*param || mean==0 || it>fMaxTdc-50) && it>MinTdc); if(MinTdc+10>it)param-=0.01; //cout<<param<<endl; }while (MinTdc+10>it); if(MinTdc+10>it && iw==10) cout<<"FindOffsets: procedure failed for (im,ip) "<<im<<" "<<ip<<endl; fOffset[im][ip][iw]=it+10; // Is it realy enough? //cout<<"im,ip,iw,offset "<<im<<" "<<ip<<" "<<iw<<" "<<it<<" "<<id<<endl; //if(im==2 && ip==6)cout<<" "<<it<<endl; //if(im==2 && ip==6){Int_t ii; cin>>ii;} } //iw if(Param!=param)cout<<"FindOffsets: Offset finding param corrected(im,ip): " <<im<<" "<<ip<<" "<<Param<<" "<<param<<endl; } //ip } //im } //______________________________________________________________________ void BrDcCalModule::CorrectOffsets() { if(!strcmp("T4",GetName())){ Int_t OffsetCorrD[]= {40, 35, 35, 0, 47, 0, 47, 20, 57, 0, 0, 45, 47, 60, 40, 14, 40, 50, 0, 50, 40, 24, 40, 15}; Int_t OffsetCorrP[]= { 0, 24, 0, 18, 22, 20, 20, 19, 18, 20, 20, 0, 0, 20, 0, 16, 0, 28, 23, 0, 22, 24, 15, 18}; for(Int_t im=0;im<3;im++) { for(Int_t ip=0;ip<fNumPlanes;ip++) { fDCDriver->AtPlane(ip); for(Int_t iw=0;iw<fDCDriver->GetMaxWireNum();iw++) { Int_t id=fDelayed[im][ip][iw]; if(id)fOffset[im][ip][iw] -= OffsetCorrD[im*fNumPlanes + ip]; else fOffset[im][ip][iw] -= OffsetCorrP[im*fNumPlanes + ip]; } //iw } //ip } //im } } //_____________________________________________________________________ void BrDcCalModule::SaveTdcOffsets() { // This method write TDC offsets to the disk ofstream OUT(fOffsetFile.Data(),ios::out); OUT.setf(ios::showpoint); cout<<"Writing Offsets Data"<<endl; for(Int_t im=0;im<3;im++) { for(Int_t ip=0;ip<fNumPlanes;ip++) { fDCDriver->AtPlane(ip); for(Int_t iw=0;iw<fDCDriver->GetMaxWireNum();iw++) { if(fOffset[im][ip][iw])OUT<<im<<" "<<ip<<" "<<iw<<" "<<fOffset[im][ip][iw]<<endl; } //iw } //ip } //im OUT.close(); } //__________________________________________________________________ void BrDcCalModule::ConstructTimeMatrix(Int_t opt) { // here we construct time matrix for(Int_t im=0;im<3;im++) { for(Int_t ip=0;ip<fNumPlanes;ip++) { fDCDriver->AtPlane(ip); for(Int_t iw=0;iw<fDCDriver->GetMaxWireNum();iw++) { Int_t Offset=fOffset[im][ip][iw]; Int_t MinTdc=fMinTdc; if(!strcmp("T4",GetName()) || !strcmp("T5",GetName())){ if(!fDelayed[im][ip][iw])MinTdc=fPromptDelay; } for(Int_t it=0;it<720;it++) { Int_t il = Offset-it; //if(Offset<4600 && Offset>3000) { // corrections for strange working channels // if(il<4096 && il+512>Offset)il=il+512; //} // Now we can jump from tdc to drift time: if(il>MinTdc) fTimeMatrix[im][ip][iw][it]=fTdcMatrix[im][ip][iw][il]; } //it } //iw } //ip } //im if(opt) { cout<<"Filling Calib. Histograms, it can takes a while"<<endl; Char_t HistName[64]; for(Int_t im=0;im<3;im++) { for(Int_t ip=0;ip<fNumPlanes;ip++) { // time sprintf(HistName,"%sTime_m%d_p%d",GetName(),im+1,ip+1); fTimeDist[im][ip] = new TH1F(HistName,HistName,500,0,1000); // to monitor time fTimeDist[im][ip]->SetYTitle("dN/dt");fTimeDist[im][ip]->SetXTitle("TDC Channel"); sprintf(HistName,"%sTime_m%d_p%dD",GetName(),im+1,ip+1); fTimeDistD[im][ip] = new TH1F(HistName,HistName,500,0,1000); // to monitor time fTimeDistD[im][ip]->SetYTitle("dN/dt");fTimeDist[im][ip]->SetXTitle("TDC Channel"); sprintf(HistName,"%sTime_m%d_p%dP",GetName(),im+1,ip+1); fTimeDistP[im][ip] = new TH1F(HistName,HistName,500,0,1000); // to monitor time fTimeDistP[im][ip]->SetYTitle("dN/dt");fTimeDist[im][ip]->SetXTitle("TDC Channel"); // Offset sprintf(HistName,"%sOffset_m%d_p%d",GetName(),im+1,ip+1); fhOffset[im][ip] = new TH1F(HistName,HistName,48,0.5,48.5); // to monitor time fhOffset[im][ip]->SetYTitle("TDC Offset");fhOffset[im][ip]->SetXTitle("Wire Number"); fDCDriver->AtPlane(ip); for(Int_t iw=0;iw<fDCDriver->GetMaxWireNum();iw++) { fhOffset[im][ip]->Fill(iw+1,fOffset[im][ip][iw]); //if(im==2 && ip==6)cout<<HistName<<" "<<fOffset[im][ip][iw]<<endl; Int_t id=fDelayed[im][ip][iw]; for(Int_t it=0;it<1000;it++) { fTimeDist[im][ip]->Fill(it,fTimeMatrix[im][ip][iw][it]); if(id)fTimeDistD[im][ip]->Fill(it,fTimeMatrix[im][ip][iw][it]); else fTimeDistP[im][ip]->Fill(it,fTimeMatrix[im][ip][iw][it]); } //it } //iw } //ip } //im } } //__________________________________________________________________ void BrDcCalModule::Calibrate() { Int_t CalibChan=fCalibChan; for(Int_t im=0;im<3;im++) { for(Int_t ip=0;ip<fNumPlanes;ip++) { TimeIntegral(im,ip); // here we integrate over wires in plane //TimeIntegral(im,ip,iw); } } if(CalibChan) { Char_t HistName[64]; for(Int_t im=0;im<3;im++) { for(Int_t ip=0;ip<fNumPlanes;ip++) { sprintf(HistName,"%sInteg_m%d_p%d",GetName(),im+1,ip+1); fInteg[im][ip] = new TH1F(HistName,HistName,500,0,1000); // monitor calibration fInteg[im][ip]->SetYTitle("drift distance [cm]"); fInteg[im][ip]->SetXTitle("TDC Channel"); //cout<<"im,ip "<<im+1<<" "<<ip+1<<" "<<fInteg[im][ip]<<endl; for(Int_t it=0;it<1000;it++) { fCalib[im][ip][it]=(Float_t)fIntegral[im][ip][it]/ (Float_t)fIntegral[im][ip][CalibChan] * fRefDriftDist; Float_t factor=1.0/fInteg[im][ip]->GetBinWidth(1); fInteg[im][ip]->Fill(it,factor*fCalib[im][ip][it]); } } } } } //__________________________________________________________________ void BrDcCalModule::TimeIntegral(Int_t im,Int_t ip,Int_t iw) { // This method produced integrated spectrum for // the chosen wire for(Int_t it=0;it<1000;it++) { Int_t Sum=0; for(Int_t i=0;i<it;i++) { Sum += fTimeMatrix[im][ip][iw][i]; } fIntegral[im][ip][it] = Sum; } } //__________________________________________________________________ void BrDcCalModule::TimeIntegral(Int_t im,Int_t ip) { // This method produced integrated spectrum for // the chosen wire for(Int_t it=0;it<1000;it++) { Int_t Sum=0; for(Int_t i=0;i<it;i++) { for(Int_t iw=1;iw<fWires[ip]-1;iw++){ // the most side wires are exluded if(!fBadWire[im][ip][iw])Sum += fTimeMatrix[im][ip][iw][i]; } } fIntegral[im][ip][it] = Sum; } } //_____________________________________________________________________ void BrDcCalModule::SaveCalibData() { // This method write matrix to the disk ofstream OUT(fCalibFile.Data(),ios::out); OUT.setf(ios::showpoint); cout<<"Writing Calibration Data"<<endl; for(Int_t im=0;im<3;im++) { for(Int_t ip=0;ip<fNumPlanes;ip++) { for(Int_t it=0;it<1000;it++) { if(fCalib[im][ip][it])OUT<<im<<" "<<ip<<" "<<it<<" " <<fCalib[im][ip][it]<<endl; } //it } //ip } //im OUT.close(); } void BrDcCalModule::DecodeRawData(BrEventNode* inNode) { // here we decode a raw data event and keep information in // private member matrix fEveTab BrDataTable *dctable_p; if ((dctable_p = inNode->GetDataTable(Form("DigDC %s", GetName())))) { Int_t nhit = dctable_p->GetEntries(); if(nhit>1999) return; //cout<<"Hits Numb: "<<nhit<<endl; Int_t ii=0; for(Int_t i = 0; i < nhit; i++) { BrDcDig* digdc_p = (BrDcDig*)dctable_p->At(i); Int_t im = digdc_p->GetSubModule() - 1; Int_t ip = digdc_p->GetPlane() - 1; Int_t iw = digdc_p->GetWire() - 1; Int_t it = digdc_p->GetTime(); Int_t il = digdc_p->GetWidth(); //cout<<im<<" "<<ip<<" "<<iw<<" "<<it<<" "<<il<<endl; // some cleaning: // T3 section if(!strcmp("T3",GetName())){ if(fDCDriver->HitIsOK(it,il)){ fEveTab[ii][0]=im; fEveTab[ii][1]=ip; fEveTab[ii][2]=iw; fEveTab[ii][3]=it; fEveTab[ii][4]=il; ii++; } } // T4/T5 section if(!strcmp("T4",GetName()) || !strcmp("T5",GetName())){ Int_t PromptDelay=fPromptDelay; if(fDCDriver->HitIsOK(it,il)){ if(PromptDelay){ if(iw<100){ if(it>PromptDelay){ fEveTab[ii][0]=im; fEveTab[ii][1]=ip; fEveTab[ii][2]=iw; fEveTab[ii][3]=it; fEveTab[ii][4]=il; fDelayed[im][ip][iw]=0; ii++; } }else{ if(it<=PromptDelay){ iw=iw%100; fEveTab[ii][0]=im; fEveTab[ii][1]=ip; fEveTab[ii][2]=iw; fEveTab[ii][3]=it; fEveTab[ii][4]=il; fDelayed[im][ip][iw]=1; ii++; } } //cout<<im<<" "<<ip<<" "<<iw<<" "<<it<<" "<<il<<endl; }else{ fEveTab[ii][0]=im; fEveTab[ii][1]=ip; fEveTab[ii][2]=iw%100; fEveTab[ii][3]=it; fEveTab[ii][4]=il; fDelayed[im][ip][iw]=iw/100; ii++; } } } } fNumCleanHits=ii; } } //_____________________________________________________________________________ void BrDcCalModule::DecodeRedData(BrEventNode* inNode) { // here we decode a reduced data event and keep information in // private member matrix fEveTab //inNode->ListObjects(); BrDataTable *AssoHitsTable; if ((AssoHitsTable = inNode->GetDataTable(Form("AssociatedHits %s", GetName())))) { Int_t nhit = AssoHitsTable->GetEntries(); Int_t ii=0; Int_t hitsAssoWithThisTrack=0; Int_t IdSave=0; for(Int_t i=0;i<nhit;i++){ BrDCHit* hit = (BrDCHit*)AssoHitsTable->At(i); Int_t im = hit->GetImod()/100 - 1; Int_t Id=hit->GetID(); if(Id != IdSave)fhAllTracks->Fill(1); if(Id != IdSave && i){ // we have a new track fhAssoHits->Fill(hitsAssoWithThisTrack); if(hitsAssoWithThisTrack<fMinAssoHits)ii=ii-hitsAssoWithThisTrack; else fhAssoHitsAcce->Fill(hitsAssoWithThisTrack); hitsAssoWithThisTrack=0; } for(Int_t j=0;j<hit->GetNumCombined();j++){ if(ii>1999)continue; fEveTab[ii][0]=im; fEveTab[ii][1]=hit->GetPlane()[j]; fEveTab[ii][2]=hit->GetWire()[j]; fEveTab[ii][3]=hit->GetTdc()[j]; // tdcUseful = fT0 - fEveTab[ii][3] //cout<<im<<" "<<hit->GetPlane()[j]<<" "<<hit->GetTdc()[j]<<endl; Int_t plane = im*fNumPlanes + fEveTab[ii][1]+1; fhHitsPerPlane->Fill(plane); hitsAssoWithThisTrack++; ii++; } IdSave=hit->GetID(); } fhAssoHits->Fill(hitsAssoWithThisTrack); if(hitsAssoWithThisTrack<fMinAssoHits)ii=ii-hitsAssoWithThisTrack; else fhAssoHitsAcce->Fill(hitsAssoWithThisTrack); if(ii<2000)fNumCleanHits=ii; else fNumCleanHits=0; } } //____________________________________________________________________ void BrDcCalModule::Print(Option_t* option) const { // Print module information // See BrModule::Print for options. // In addition this module defines the Option: // <fill in here> TString opt(option); opt.ToLower(); BrModule::Print(option); if (opt.Contains("d")) cout << endl << " Original author: Pawel Staszel" << endl << " Last Modifications: " << endl << " $Author: ouerdane $" << endl << " $Date: 2002/05/07 17:35:24 $" << endl << " $Revision: 1.17 $ " << endl << endl << "-------------------------------------------------" << endl; } //------------------------------------------------ // $Log: BrDcCalModule.cxx,v $ // Revision 1.17 2002/05/07 17:35:24 ouerdane // added driftTime parameter to DB and classes // // Revision 1.16 2002/05/02 14:59:41 ufstasze // Added fLocalTrack to retrieve global tdc offset from db // // Revision 1.15 2002/03/07 17:49:19 ouerdane // added BrDcCalibration and connection to DB // // Revision 1.14 2002/02/15 16:09:28 ufstasze // Changes to method that looks for Tdc offsets. // // Revision 1.13 2002/02/06 18:53:56 ufstasze // updates that have to do with pp calibration // // Revision 1.12 2001/12/07 21:22:04 videbaek // Change UserForRead to Use( ,mode,..) // // Revision 1.11 2001/11/19 17:32:48 ufstasze // added FindTdcOffsetsNew // // Revision 1.10 2001/11/02 15:29:20 staszel // Reorganized data decoding. Now decoding can be done by using // DecodeRedData or DecodeRawData depending on the value of fCalibLevel // parameter. // // Revision 1.9 2001/10/08 11:27:38 cholm // Changed to use new DB access classes // // Revision 1.8 2001/09/07 14:27:33 staszel // Added new method FindTdcOffsets1New (rather temporary). // // Revision 1.6 2001/08/23 11:28:29 staszel // Some diagnostic histograms plus CorrectOffsets method. // // Revision 1.5 2001/08/13 11:46:35 staszel // Same small changes. // // Revision 1.4 2001/07/30 11:37:47 staszel // developments regarding T4 and T5 calibration // // Revision 1.3 2001/07/20 17:25:58 krakow // Changes related to BrDriverDC // // Revision 1.2 2001/06/22 17:38:50 cholm // Minor changes // // Revision 1.1.1.1 2001/06/21 14:55:05 hagel // Initial revision of brat2 // // Revision 1.2 2001/06/14 13:03:53 staszel // first cosmetic upgrade // // Revision 1.1 2001/06/14 12:31:53 staszel // Initial release of DC calibration code // |
||||||
|
This page automatically generated by script docBrat by Christian Holm |
Copyright ; 2002 BRAHMS Collaboration
<brahmlib@rcf.rhic.bnl.gov>
|