|
/////////////////////////////////////////////////////////////// // // BrDriverDC is a BRAHMS class for serving the real geometry // of T3, T4 and T5 detectors. This geometry directly refers // to the construction properties of these detectors. // NOTE: the units assumed in this code are in [cm]. // /////////////////////////////////////////////////////////////// // // $Id: BrDriverDC.cxx,v 1.16 2002/08/08 15:39:32 ufstasze Exp $ // $Author: ufstasze $ // $Date: 2002/08/08 15:39:32 $ #include "BrDriverDC.h" #ifndef BRAT_BrPlane3D #include "BrPlane3D.h" #endif #ifndef BRAT_BrVector3D #include "BrVector3D.h" #endif ClassImp(BrDriverDC); // //___________________________________________________________________ BrDriverDC::BrDriverDC(const Char_t *Name,const Char_t *Title) :TNamed(Name,Title) { // In this constructor some member parameter are set // to reflect the DC's construction. fNormalMode=0; // normal_mode=0 (default) when working with mapped DC data fZPosRef=1; // Assume that track parameters that go through SetTrack are defined // in the local SubModule frame. // When fZPosRef==0 then we assume that track parameters are defined // in the Detector frame. fRealSetup=kTRUE; Init(); } //___________________________________________________________________ void BrDriverDC::Init() { // Here we define some parameters that describe construction of the DCs. //Bool_t RealSetup=kFALSE; // very temporary, as long as GEANT geometry is wrong Bool_t RealSetup=kTRUE; if(!strcmp("T3",GetName())) {// case when detector = 'T3' fMaxPlaneNum=10; fWireDist=1.0; fRefDist=0.5; fActiveRangeX=41.0; fActiveRangeY=31.0; fXDCWide=fActiveRangeX*0.5; fYDCWide=fActiveRangeY*0.5; if(RealSetup){ fModuleDisp=25.5; // valid for real setup fCenterToWindowDist=-32.25; // valid for real setup }else{ fModuleDisp=35.0; // assumed in gean simul. fCenterToWindowDist=-(6.75+fModuleDisp); } fBaseWirePosXA=6.2; fBaseWirePosXB=6.45; fBaseWirePosXC=6.2; fBaseWirePosYA=-6.2; fBaseWirePosYB=-6.45; fBaseWirePosYC=-6.2; fBaseWirePosUA=7.59207; fBaseWirePosUB=7.84207; fBaseWirePosVA=-5.1378; fBaseWirePosVB=-5.3878; fXmin=-46.5; fXmax=-5.5; fYmin=-36.5; fYmax=-5.5; fDCX0=26.0; fDCY0=21.0; fView[0]=0; fView[1]=0; fView[2]=0; fView[3]=1; fView[4]=1; fView[5]=1; fView[6]=2; fView[7]=2; fView[8]=3; fView[9]=3; } if(!strcmp("T4",GetName())) { // case when detector='T4' fMaxPlaneNum=8; fWireDist=2.2; fRefDist=1.05; if(RealSetup){ fModuleDisp=22.5; // valid for real setup fCenterToWindowDist=-27.75; // valid for real setup }else{ fModuleDisp=33.0; // assumed in gean simul. fCenterToWindowDist=-(5.25+fModuleDisp); } fActiveRangeX=51.0; fActiveRangeY=36.0; fXDCWide=fActiveRangeX*0.5; fYDCWide=fActiveRangeY*0.5; fBaseWirePosXA=6.8; fBaseWirePosXB=7.35; fBaseWirePosYA=-6.8; fBaseWirePosYB=-7.35; fBaseWirePosUA=8.53618; fBaseWirePosUB=9.08618; fBaseWirePosVA=-5.67072; fBaseWirePosVB=-6.22072; fXmin=-56.5; fXmax=-5.5; fYmin=-41.5; fYmax=-5.5; fDCX0=31.0; fDCY0=23.5; fView[0]=0; fView[1]=0; fView[2]=1; fView[3]=1; fView[4]=2; fView[5]=2; fView[6]=3; fView[7]=3; } if(!strcmp("T5",GetName())) { // case when detector='T5' fMaxPlaneNum=8; fWireDist=2.2; fRefDist=1.05; if(RealSetup){ fModuleDisp=22.5; // valid for real setup fCenterToWindowDist=-27.75; // valid for real setup }else{ fModuleDisp=33.0; // assumed in gean simul. fCenterToWindowDist=-(5.25+fModuleDisp); } fActiveRangeX=51.0; fActiveRangeY=36.0; fXDCWide=fActiveRangeX*0.5; fYDCWide=fActiveRangeY*0.5; fBaseWirePosXA=6.8; fBaseWirePosXB=7.35; fBaseWirePosYA=-6.8; fBaseWirePosYB=-7.35; fBaseWirePosUA=8.53618; fBaseWirePosUB=9.08618; fBaseWirePosVA=-5.67072; fBaseWirePosVB=-6.22072; fXmin=-56.5; fXmax=-5.5; fYmin=-41.5; fYmax=-5.5; fDCX0=31.0; fDCY0=23.5; fView[0]=0; fView[1]=0; fView[2]=1; fView[3]=1; fView[4]=2; fView[5]=2; fView[6]=3; fView[7]=3; } DefineCorrections(); } //___________________________________________________________________ void BrDriverDC::ShowWord() { cout<<fSlowo[0]<<fSlowo[1]<<fSlowo[2]<<fSlowo[3]<<fSlowo[4]<<endl; } //___________________________________________________________________ Float_t BrDriverDC::ZfromDCFront( ) { Float_t z = fPlaneNumber * 1.5 + fModuleNumber*fModuleDisp; return z; } //___________________________________________________________________ void BrDriverDC::ShowPlaneSpecif() { // prints plane specifications on STDOUT cout << "Plane Specification: " << fSlowo[0] << fSlowo[1] << fSlowo[2] << fSlowo[3] << fSlowo[4] <<" "<<GetName()<< endl; } Float_t BrDriverDC::GetAngle() { Float_t Angle=1000.0; if(fSlowo[2]=='X')Angle=0.0; if(fSlowo[2]=='Y')Angle=-90.0; if(fSlowo[2]=='U')Angle=-18.0; if(fSlowo[2]=='V')Angle=+18.0; return Angle; } Float_t BrDriverDC::SubModulePlanePosition() { // This Function returns the plane Z position. // This position is in respect to the center point of // the module we are just now At. Float_t SubplanePos; SubplanePos = PlanePosition() - (fModuleNumber-1)*fModuleDisp; return SubplanePos; } Float_t BrDriverDC::PlanePosition() { // This Function returns the Z position of plane. // This position is in respect to the center point of // the whole DC detector (T3, T4 or T5). Float_t z, planePos; z = ZfromDCFront(); planePos = z + fCenterToWindowDist; return planePos; } Float_t BrDriverDC::GetPlanePosition(Int_t imod, Int_t iplane) { // This Function returns the Z position of plane. // This position is in respect to the center point of // the whole DC detector (T3, T4 or T5). Float_t z, planePos; At(imod,iplane); z = ZfromDCFront(); planePos = z + fCenterToWindowDist; return planePos; } Int_t BrDriverDC::VirtualNumber(Int_t wire) { // This function returns the wire position which is // the distance between the wire and the center of the // West-North Positioning Pin. Int_t chanOrd=(fNumberOfWires-1)-wire; return chanOrd; } Float_t BrDriverDC::WirePosition(Int_t wire) { // This function returns the wire position which is // the distance between the wire and the center of the // West-North Positioning Pin. Float_t distance=1000.0; Int_t chanOrd=VirtualNumber(wire); // X planes: if(fSlowo[2]=='X') { if(fSlowo[3]=='A')distance = fBaseWirePosXA + chanOrd*fWireDist; if(fSlowo[3]=='C')distance = fBaseWirePosXA + chanOrd*fWireDist; if(fSlowo[3]=='B')distance = fBaseWirePosXB + chanOrd*fWireDist; } // Y planes: if(fSlowo[2]=='Y') { if(fSlowo[3]=='A')distance = fBaseWirePosYA - chanOrd*fWireDist; if(fSlowo[3]=='C')distance = fBaseWirePosYA - chanOrd*fWireDist; if(fSlowo[3]=='B')distance = fBaseWirePosYB - chanOrd*fWireDist; } // U planes: if(fSlowo[2]=='U') { if(fSlowo[3]=='A')distance = fBaseWirePosUA + chanOrd*fWireDist; if(fSlowo[3]=='B')distance = fBaseWirePosUB + chanOrd*fWireDist; } // V planes: (most complicated) if(fSlowo[2]=='V') { if(fSlowo[3]=='A')distance = fBaseWirePosVA + chanOrd*fWireDist; if(fSlowo[3]=='B')distance = fBaseWirePosVB + chanOrd*fWireDist; } return distance; } Int_t BrDriverDC::WireStatus(Int_t wire) { // This function returns status of the particular channel: // 0 - when the channel 'GlobalChan' is empty // 1 - when the channel 'GlobalChan' is connected to existed wire Int_t chanOrd, status; status=0; // default 0 value chanOrd=wire; if(chanOrd >= 0) { if(!strcmp("T3",GetName())) { if(fSlowo[2]=='X') if(chanOrd<40)status=1; if(fSlowo[2]=='Y') if(chanOrd<30)status=1; if(fSlowo[2]=='U') if(chanOrd<47)status=1; if(fSlowo[2]=='V') if(chanOrd<48)status=1; } // closing T3 section if(!strcmp("T4",GetName()) || !strcmp("T5",GetName())) { if(fSlowo[2]=='X') if(chanOrd<23)status=1; if(fSlowo[2]=='Y') if(chanOrd<16)status=1; if(fSlowo[2]=='U') if(chanOrd<27)status=1; if(fSlowo[2]=='V') if(chanOrd<27)status=1; } // closing T4/T5 section } // if(chanOrd >=0) return status; } Float_t BrDriverDC::WireEndPoints(Int_t wire) { // This function calculate wire end points coordinates and // return the length of wire Float_t x1=0, x2=0, xx1, xx2, y1=0, y2=0, yy1, yy2, d, b; Float_t WireLength; Float_t z=PlanePosition(); Float_t dist=WirePosition(wire); if(fSlowo[2]=='X') // X plane { d = -dist; x1 = d; x2 = d; y1 = fYmin; y2 = fYmax; } if(fSlowo[2]=='Y') // Y plane { d = dist; x1 = fXmin; x2 = fXmax; y1 = d; y2 = d; } if(fSlowo[2]=='U') // U plane { d = -dist; b = d/0.309017; xx1 = (fYmin-b)/(-3.0776835) + fDCX0; yy1 = -3.0776835*fXmax + b + fDCY0; if(TMath::Abs(xx1)<=fXDCWide){x1=xx1-fDCX0; y1=fYmin;} if(TMath::Abs(yy1)<=fYDCWide){x1=fXmax; y1=yy1-fDCY0;} if(TMath::Abs(xx1)<=fXDCWide && TMath::Abs(yy1)<=fYDCWide) {cout<<"Logical Error (two 'min' points)"<<xx1<<" "<<yy1<<endl;} if(TMath::Abs(xx1)>fXDCWide && TMath::Abs(yy1)>fYDCWide) {cout<<"Logical Error (no 'min' point)"<<endl;} xx2 = (fYmax-b)/(-3.0776835) + fDCX0; yy2 = -3.0776835*fXmin + b + fDCY0; if(TMath::Abs(xx2)<=fXDCWide){x2=xx2-fDCX0; y2=fYmax;} if(TMath::Abs(yy2)<=fYDCWide){x2=fXmin; y2=yy2-fDCY0;} if(TMath::Abs(xx2)<=fXDCWide && TMath::Abs(yy2)<=fYDCWide) {cout<<"Logical Error (two 'max' points)"<<endl;} if(TMath::Abs(xx2)>fXDCWide && TMath::Abs(yy2)>fYDCWide) {cout<<"Logical Error (no 'max' point)"<<endl;} } if(fSlowo[2]=='V') // V plane { d = dist; b = d/0.309017; xx1 = (fYmin-b)/(3.0776835) + fDCX0; yy1 = 3.0776835*fXmin + b + fDCY0; if(TMath::Abs(xx1)<=fXDCWide){x1=xx1-fDCX0; y1=fYmin;} if(TMath::Abs(yy1)<=fYDCWide){x1=fXmin; y1=yy1-fDCY0;} if(TMath::Abs(xx1)<=fXDCWide && TMath::Abs(yy1)<=fYDCWide) {cout<<"Logical Error (two 'min' points)"<<endl;} if(TMath::Abs(xx1)>fXDCWide && TMath::Abs(yy1)>fYDCWide) {cout<<"Logical Error (no 'min' point)"<<endl;} xx2 = (fYmax-b)/(3.0776835) + fDCX0; yy2 = 3.0776835*fXmax + b + fDCY0; if(TMath::Abs(xx2)<=fXDCWide){x2=xx2-fDCX0; y2=fYmax;} if(TMath::Abs(yy2)<=fYDCWide){x2=fXmax; y2=yy2-fDCY0;} if(TMath::Abs(xx2)<=fXDCWide && TMath::Abs(yy2)<=fYDCWide) {cout<<"Logical Error (two 'max' points)"<<endl;} if(TMath::Abs(xx2)>fXDCWide && TMath::Abs(yy2)>fYDCWide) {cout<<"Logical Error (no 'max' point)"<<endl;} } // Now calculate the Wires End Points in the center of plane frame // define by x_DC = 0 // y_DC = 0 // where x_DC, y_DC are the coordinates of the middle (center) // point of each plane. x1 = x1 + fDCX0; y1 = y1 + fDCY0; x2 = x2 + fDCX0; y2 = y2 + fDCY0; // Let, copy the local variables to the public variables available in the other BRAT moduls: fX1=x1; fY1=y1; fX2=x2; fY2=y2; fZ=z; WireLength = sqrt ( (x2-x1)*(x2-x1) + (y2-y1)*(y2-y1) ); // wire length within the DC active area return WireLength; } void BrDriverDC::AtSubModule(Int_t modulNum) { // This function set the pointer to the Submodule numbered 'modulNumb' fSlowo[1]='-'; if(modulNum==0)fSlowo[0]='1'; if(modulNum==1)fSlowo[0]='2'; if(modulNum==2)fSlowo[0]='3'; fSlowo[4]='1'; fModuleNumber=modulNum; } void BrDriverDC::AtPlane(Int_t planeNum) { // This function set the pointer to the plane numbered 'planeNum'. // This function will be changed in the future !!! if(!strcmp("T3",GetName())) // T3 detector { if(planeNum==0){fSlowo[2]='X'; fSlowo[3]='A';} if(planeNum==1){fSlowo[2]='X'; fSlowo[3]='B';} if(planeNum==2){fSlowo[2]='X'; fSlowo[3]='C';} if(planeNum==3){fSlowo[2]='Y'; fSlowo[3]='A';} if(planeNum==4){fSlowo[2]='Y'; fSlowo[3]='B';} if(planeNum==5){fSlowo[2]='Y'; fSlowo[3]='C';} if(planeNum==6){fSlowo[2]='U'; fSlowo[3]='A';} if(planeNum==7){fSlowo[2]='U'; fSlowo[3]='B';} if(planeNum==8){fSlowo[2]='V'; fSlowo[3]='A';} if(planeNum==9){fSlowo[2]='V'; fSlowo[3]='B';} if(fSlowo[2]=='X')fNumberOfWires=40; if(fSlowo[2]=='Y')fNumberOfWires=30; if(fSlowo[2]=='U')fNumberOfWires=47; if(fSlowo[2]=='V')fNumberOfWires=48; } if(!strcmp("T4",GetName()) || !strcmp("T5",GetName())) // T4/T5 detector { if(planeNum==0){fSlowo[2]='X'; fSlowo[3]='A';} if(planeNum==1){fSlowo[2]='X'; fSlowo[3]='B';} if(planeNum==2){fSlowo[2]='Y'; fSlowo[3]='A';} if(planeNum==3){fSlowo[2]='Y'; fSlowo[3]='B';} if(planeNum==4){fSlowo[2]='U'; fSlowo[3]='A';} if(planeNum==5){fSlowo[2]='U'; fSlowo[3]='B';} if(planeNum==6){fSlowo[2]='V'; fSlowo[3]='A';} if(planeNum==7){fSlowo[2]='V'; fSlowo[3]='B';} if(fSlowo[2]=='X')fNumberOfWires=23; if(fSlowo[2]=='Y')fNumberOfWires=16; if(fSlowo[2]=='U')fNumberOfWires=27; if(fSlowo[2]=='V')fNumberOfWires=27; } fPlaneNumber=planeNum; } void BrDriverDC::At(Int_t modulNum, Int_t planeNum) { // This function set the pointer to the Submodule numbered 'modulNumb' // and to the plane numbered 'planeNum'. AtSubModule(modulNum); AtPlane(planeNum); } void BrDriverDC::SetMode(Int_t i) { fNormalMode=i; } Int_t BrDriverDC::GetMode( ) { return fNormalMode; } void BrDriverDC::SetDetectorTrack(Float_t* trPar) { // Use this method when your initial track parameters are defined in // local Detector frame. fAx=*trPar++; fAy=*trPar++; fX0=*trPar++; fY0=*trPar; fTrParam[0]=fAx; fTrParam[1]=fAy; fTrParam[2]=fX0; fTrParam[3]=fY0; } void BrDriverDC::SetSubModuleTrack(Float_t* trPar) { // Use this method when your initial track parameters are defined in // local SubModule frame. fAx=*trPar++; fAy=*trPar++; fX0=*trPar++; fY0=*trPar; fTrParam[0]=fAx; fTrParam[1]=fAy; fTrParam[2]=fX0; fTrParam[3]=fY0; Float_t z = - GetSubModulePosition(); // "-" sign because we trace back fTrParam[2] = fAx*z + fX0; // now fTrParam are define in center of the detector frame fTrParam[3] = fAy*z + fY0; fX0=fTrParam[2]; fY0=fTrParam[3]; } void BrDriverDC::ShowTrack() { // This method print parameters of the actual track defined // in the center of the detector frame. //if(GetZPosRef()) { // Float_t z = (1-fModuleNumber) * fModuleDisp; // Float_t x = fAx*z + fX0; // Float_t y = fAy*z + fY0; // cout<<"track1 (ax,ay,x0,y0): "<<fAx<<" "<<fAy<<" "<<x<<" "<<y<<endl; // } // else // { // cout<<"track0 (ax,ay,x0,y0): "<<fAx<<" "<<fAy<<" "<<fX0<<" "<<fY0<<endl; // } cout<<"track (ax,ay,x0,y0): "<<fAx<<" "<<fAy<<" "<<fX0<<" "<<fY0<<endl; } Float_t* BrDriverDC::GetTrack() { // This method returns pointer to the vector // that contains the actual track parameters defined // in the center of detector frame (local detector track). ///if(GetZPosRef()) { // Float_t z = (1-fModuleNumber) * fModuleDisp; // fTrParam[2] = fAx*z + fX0; // now fTrParam are define in center of the detector frame // fTrParam[3] = fAy*z + fY0; return fTrParam; } void BrDriverDC::GoThroughModule(Int_t* wires, Float_t* xdrifts) { Int_t wire=0; Float_t xdrift=0.0; for(Int_t planeNum=0; planeNum<fMaxPlaneNum; planeNum++){ AtPlane(planeNum); GoThroughPlane(wire, xdrift); if(WireStatus(wire)){ *wires++ = wire; *xdrifts++ = xdrift; }else{ if(DebugLevel()>0) { cout << " Track out of active area, (modul, plane, wire)= " <<fModuleNumber<<" "<<fPlaneNumber<<" "<<wire<<endl; } *wires++ = wire; *xdrifts++ = xdrift; } } } void BrDriverDC::GoThroughPlane(Int_t& wire, Float_t& xdrift) { // This function looks for the wire that has the smallest distance to // the track defined by 'ax', 'ay', 'x0', 'y0'. This is for module // 'ModulNum' and plane 'planeNum' (see function 'SimulMode'). // This wire number and distance between wire and track are indicated // by 'wire' and 'xdrift', respectively. Float_t wL, xd=0, FirstToHit; Float_t a, b, x, y, xp, yp, z; Int_t iw=0; wL=WireEndPoints(iw); //wire end points for wire number = 0 //if(GetZPosRef()) // z = SubModulePlanePosition(); // else z = PlanePosition(); x = fAx*z + fX0; // x and y hit position in plane 'planeNum' y = fAy*z + fY0; // cout<<"plane="<<fPlaneNumber<<" z="<<z<<endl; // cout<<x<<" "<<0.5*ActiveRangeX<<" "<<y<<" "<<0.5*ActiveRangeY<<endl; if(TMath::Abs(x)<0.5*fActiveRangeX && TMath::Abs(y)<0.5*fActiveRangeY) { if(fSlowo[2]=='X') { FirstToHit=x-fX1; // iw = (Int_t)rint(FirstToHit/fWireDist); iw = TMath::Nint(FirstToHit/fWireDist); xd = TMath::Abs(FirstToHit-iw*fWireDist); // cout << FirstToHit/WireDist<< endl; // cout <<FirstToHit<<" "<<iw<<" "<<WireDist<<" "<<iw*WireDist<<endl; } if(fSlowo[2]=='Y') { FirstToHit=y-fY1; // iw = (Int_t)rint(FirstToHit/fWireDist); iw = TMath::Nint(FirstToHit/fWireDist); xd = TMath::Abs(FirstToHit-iw*fWireDist); } if(fSlowo[2]=='U' || fSlowo[2]=='V') { a = (fY2-fY1)/(fX2-fX1); b = fY1 - a*fX1; xp=(x+a*y-a*b)/(1+a*a); yp=a*xp + b; FirstToHit=sqrt((x-xp)*(x-xp) + (y-yp)*(y-yp)); // iw = (Int_t)rint(FirstToHit/fWireDist); iw = TMath::Nint(FirstToHit/fWireDist); xd = TMath::Abs(FirstToHit-iw*fWireDist); // cout << "track0 " << x << " " << y<<" "<<a<<" "<<FirstToHit<< " " << wL<< endl; } wire = iw; xdrift=xd; // cout << "wire " << wire << " drift " << xdrift << endl; } else { wire = -1; xdrift=-1; } } Float_t BrDriverDC::HitToWireDist(Int_t iw) { // This function return the distance between the particular track // and the wire numbered 'iw'. // The track should be previously defined by calling the function // ReadTrack with the actual track parameters: fAx, fAy, fX0, fX0. // If the hit is out the plane active area it returns value '-1' Float_t wL, htwd=1000.0; Float_t a, b, x, y, xp, yp, z; if(!WireStatus(iw)) { return -1; } wL=WireEndPoints(iw); //if(GetZPosRef()) //z = SubModulePlanePosition(); //else z = PlanePosition(); x = fAx*z + fX0; // x and y hit position in plane 'planeNum' y = fAy*z + fY0; // cout<<x<<" "<<0.5*ActiveRangeX<<" "<<y<<" "<<0.5*ActiveRangeY<<endl; if(TMath::Abs(x)<0.5*fActiveRangeX && TMath::Abs(y)<0.5*fActiveRangeY)//hit is in the plane active area { if(fSlowo[2]=='X') { htwd=TMath::Abs(fX1-x); } if(fSlowo[2]=='Y') { htwd=TMath::Abs(fY1-y); } if(fSlowo[2]=='U' || fSlowo[2]=='V') { a = (fY2-fY1)/(fX2-fX1); b = fY1 - a*fX1; xp=(x+a*y-a*b)/(1+a*a); yp=a*xp + b; htwd=sqrt((x-xp)*(x-xp) + (y-yp)*(y-yp)); } } else { htwd=-1; } return htwd; } //_______________________________________________________________________ void BrDriverDC::DefineCorrections() { if(!strcmp("T3",GetName())){ Float_t PlaneCorr0[]= {0.0015,-0.0011, 0.0015, 0.0050, 0.0013, 0.0010, 0.0070, 0.0080, 0.0087, 0.0135, 0.0026,-0.0046,-0.0035,-0.0036,-0.0099,-0.0078,-0.0125,-0.0115,-0.0080,-0.0125, -0.0044,-0.0032, 0.0020, 0.0007,-0.0005, 0.0069, 0.0101, 0.0073, 0.0077, 0.0112}; Float_t PlaneCorr1[]= {-0.0062,-0.0050,-0.0035,-0.0006, 0.0006,-0.0014, 0.0006, 0.0016, 0.0040, 0.0037, 0.0007, 0.0069, 0.0061,-0.0004, 0.0015, 0.0017, 0.0048, 0.0016, 0.0001,-0.0018, -0.0072,-0.0057,-0.0048, 0.0002, 0.0011,-0.0004, 0.0005, 0.0002, 0.0048, 0.0051}; for(Int_t i=0;i<3*fMaxPlaneNum;i++)fPlaneCorr[i]=PlaneCorr0[i]+PlaneCorr1[i]; } if(!strcmp("T4",GetName())){ Float_t PlaneCorr0[]= {0.00987, 0.0076, 0.0300, 0.0413,-0.0880,-0.0838, 0.0899, 0.1000, 0.0098, 0.0076, 0.0590, 0.0502,-0.0798,-0.0868, 0.0836, 0.0902, 0.0037, 0.0045, 0.0580, 0.0583,-0.0895,-0.0917, 0.0647, 0.0785}; Float_t PlaneCorr1[]= {0.0036, 0.0006,-0.0015,-0.0002,-0.0041,-0.0075, 0.0119, 0.0052, 0.0020,-0.0014, 0.0061, 0.0115,-0.0094,-0.0047, 0.0015, 0.0076, -0.0033,-0.0021, 0.0029, 0.0037,-0.0066,-0.0066, 0.0113, 0.0099}; Float_t PlaneCorr2[]= {0.0024, 0.0009,-0.0008,-0.0006,-0.0011,-0.0006, 0.0062, 0.0000, -0.0023,-0.0011, 0.0020, 0.0052,-0.00024,-0.0005, 0.0070,-0.0012, -0.0008,-0.0038, 0.0015, 0.0027,-0.0025,-0.0023, 0.0130, 0.0028}; for(Int_t i=0;i<3*fMaxPlaneNum;i++)fPlaneCorr[i]=PlaneCorr0[i]+PlaneCorr1[i]+PlaneCorr2[i]; } if(!strcmp("T5",GetName())){ Float_t PlaneCorr0[]= {0.00968, 0.0025, 0.0263, 0.0389,-0.0724,-0.0797, 0.0788, 0.0896, -0.0057,-0.0023, 0.0254, 0.0533,-0.0826,-0.0859, 0.0801, 0.0841, 0.0008, 0.0028, 0.0407, 0.0559,-0.0700,-0.0788, 0.0625, 0.0609}; Float_t PlaneCorr1[]= {0.0027, 0.0022, 0.0093,-0.0017,-0.0165,-0.0129, 0.0119, 0.0135, 0.0080, 0.0013, 0.0187, 0.0013,-0.0136,-0.0129, 0.0131, 0.0136, -0.0026,-0.0026, 0.0142, 0.0059,-0.0168,-0.0133, 0.0150, 0.0201}; Float_t PlaneCorr2[]= {-0.0037, 0.0008, 0.0045, 0.0026,-0.0036,-0.0048, 0.0156, 0.0051, -0.0031,-0.0004, 0.0000,-0.0011, 0.0000, 0.0015, 0.0055, 0.0030, 0.0019,-0.0016, 0.0064, 0.0035,-0.0057,-0.0050, 0.0089, 0.0071}; for(Int_t i=0;i<3*fMaxPlaneNum;i++)fPlaneCorr[i]=PlaneCorr0[i]+PlaneCorr1[i]+PlaneCorr2[i]; } } //______________________________________________________________ Int_t BrDriverDC::LineAcceptanceStatus(BrLine3D line,Int_t opt) { // This method returns line acceptance status defined in the // following way: // 0 - line out of detector acceptance // 1 - line within acceptance only at detector entry or only // at detector exit // 2 - line within acceptance at entry and at exit Float_t z_in = GetPlanePosition(0,0); Float_t z_out = GetPlanePosition(2,fMaxPlaneNum-1); BrPlane3D PlaneIn(0,0,z_in, 0,1,z_in, 1,0,z_in); BrPlane3D PlaneOut(0,0,z_out, 0,1,z_out, 1,0,z_out); BrVector3D projIn(PlaneIn.GetIntersectionWithLine(line)); BrVector3D projOut(PlaneOut.GetIntersectionWithLine(line)); // define XMaxIn and YMaxIn At(0,0); Float_t Xmax = TMath::Min(TMath::Abs(GetWirePosition(0+opt)), TMath::Abs(GetWirePosition(fNumberOfWires-1-opt))); At(0,3); Float_t Ymax = TMath::Min(TMath::Abs(GetWirePosition(0+opt)), TMath::Abs(GetWirePosition(fNumberOfWires-1-opt))); Int_t LineStatus=0; if(TMath::Abs(projIn.GetX())<Xmax && TMath::Abs(projIn.GetY())<Ymax)LineStatus++; if(TMath::Abs(projOut.GetX())<Xmax && TMath::Abs(projOut.GetY())<Ymax)LineStatus++; return LineStatus; } //__________________________________________________________________ Float_t BrDriverDC::GetWirePosition(Int_t iw) { // This function return the distance between the center of the plane // and the wire numbered 'iw'. // first define correction for planes Float_t exception=0.0; if(!strcmp("T3",GetName())&&fRealSetup){ if(fModuleNumber==2 && fPlaneNumber==7 && iw==22)exception=-0.035; if(fModuleNumber==2 && fPlaneNumber==7 && iw==34)exception=-0.017; } Float_t Corr; if(fRealSetup)Corr=fPlaneCorr[fModuleNumber*fMaxPlaneNum + fPlaneNumber]; else Corr=0; Float_t wL, htwd=1000.0; Float_t a, b, x, y, xp, yp; if(!WireStatus(iw)){ cout<<" Error: there is no such a wire in this plane"<<endl; return 0; } wL=WireEndPoints(iw); x = 0.0; // x and y hit position in plane 'planeNum' (here always 0) y = 0.0; if(fSlowo[2]=='X') { // htwd=TMath::Abs(fX1-x); htwd=fX1; } if(fSlowo[2]=='Y') { // htwd=TMath::Abs(fY1-y); htwd=fY1; } if(fSlowo[2]=='U' || fSlowo[2]=='V') { a = (fY2-fY1)/(fX2-fX1); b = fY1 - a*fX1; xp=(x+a*y-a*b)/(1+a*a); yp=a*xp + b; if(fSlowo[2]=='U')htwd=sqrt((x-xp)*(x-xp) + (y-yp)*(y-yp)) * b/TMath::Abs(b); if(fSlowo[2]=='V')htwd=-sqrt((x-xp)*(x-xp) + (y-yp)*(y-yp)) * b/TMath::Abs(b); } return (htwd - Corr) - exception; } // methods used to clean the data //____________________________________________________________________ Bool_t BrDriverDC::ReadCuts() { // Here we open a file with cuts for a given run // This cuts alow to remove noisy areas in the // Time vs SignaWidth plane. if(frunNumber){ Char_t CutsFile[64]; sprintf(CutsFile,"run%d.cuts",frunNumber); ifstream IN; IN.open(CutsFile); if(IN){ cout<<"file "<<CutsFile<<" Opened"<<endl; }else{ cout<<"no CutsFile existed for this run"<<endl; cout<<"opening file for run 3128"<<endl; Int_t number = 3128; IN.close(); IN.clear(); sprintf(CutsFile,"run%d.cuts",number); IN.open(CutsFile); if(IN)cout<<"file "<<CutsFile<<" Opened"<<endl; else cout<<"couldn't open "<<CutsFile<<" check it out"<<endl; } IN.ignore(150,'n'); IN.ignore(150,'n'); IN.ignore(150,'n'); IN.ignore(150,'n'); Int_t im,ip,ncuts; for(Int_t j=0;j<30;j++){ IN>>im>>ip>>ncuts; fnCuts[im][ip]=ncuts; cout<<"im="<<im<<" "<<"ip="<<ip<<" "<<"ncuts="<<ncuts<<endl; for(Int_t i=0;i<ncuts;i++){ //cout<<"i="<<i<<endl; IN >> fTdcLowCut[im][ip][i] >> fTdcUppCut[im][ip][i]; if(fTdcLowCut[im][ip][i] > fTdcUppCut[im][ip][i])cout<<"Time: Lower Cut > Upper Cut"<<endl; } for(Int_t i=0;i<ncuts;i++){ IN >> fWidthLowCut[im][ip][i] >> fWidthUppCut[im][ip][i]; if(fWidthLowCut[im][ip][i] > fWidthUppCut[im][ip][i])cout<<"Width: Lower Cut > Upper Cut"<<endl; //cout<<fTdcLowCut[im][ip][i]<<endl; //cout<<fWidthLowCut[im][ip][i]<<endl; } } IN.close(); return kTRUE; } cout<<"BrDriverDC::ReadCuts: run number has to be set before calling this method"<<endl; return kFALSE; } //________________________________________________________________________ Bool_t BrDriverDC::HitIsOK(Int_t it,Int_t iw) { // This method check whether a given hit is from noise or from // clean Time vs HitWidth region. // It also impose common requirement on hit's width and hit's time //cout<<"In HitIsOK: "<<it<<" "<<GetMinTdc()<<" "<<GetMaxTdc()<<endl; Bool_t isOK=kTRUE; if(it<GetMinTdc()) isOK = kFALSE; if(it>GetMaxTdc()) isOK = kFALSE; if(!strcmp("T3",GetName())){ //cout<<"In HitIsOK: "<<iw<<" "<<GetMinWidth()<<endl; // this analysis makes sens only for T3 detector if(iw<GetMinWidth()) return kFALSE; } return isOK; } // methods for calibration purpose //_____________________________________________________________________ void BrDriverDC::ReadCalibData(TString CalibFile,Int_t CalibChan) { // This method read matrix from the disk Int_t im,ip,iw,it; Float_t r; Int_t CalibChanInView[4]; if(!strcmp("T3",GetName())){ CalibChanInView[0]=CalibChan; CalibChanInView[1]=CalibChan; CalibChanInView[2]=CalibChan; CalibChanInView[3]=CalibChan; } if(!strcmp("T4",GetName())){ CalibChanInView[0]=CalibChan; CalibChanInView[1]=CalibChan; CalibChanInView[2]=CalibChan; CalibChanInView[3]=CalibChan+12; if(frunNumber==5544 || frunNumber==5542 || frunNumber==5541){ fException[0][2]=10; fException[0][3]=10; fException[1][2]=-9; fException[1][3]=-9; fException[1][4]=-9; fException[1][5]=-6; } } if(!strcmp("T5",GetName())){ CalibChanInView[0]=CalibChan; CalibChanInView[1]=CalibChan; CalibChanInView[2]=CalibChan; CalibChanInView[3]=CalibChan+12; } cout<<"Calibration Channels (in X,Y,U,V): "<<CalibChanInView[0]<<" " <<CalibChanInView[1]<<" "<<CalibChanInView[2]<<" "<<CalibChanInView[3]<<endl; ifstream IN(CalibFile.Data()); cout<<"calib file: "<<CalibFile<<endl; if(IN){ cout<<"Reading Calib Data"<<endl; do { IN>>im>>ip>>it>>r; if(!(IN.eof()))fDCCalib[im][ip][it]=r; } while(!(IN.eof())); IN.close(); }else{ cout<<"can not open calib file: "<<CalibFile<<endl; } for(Int_t im=0;im<3;im++) { for(Int_t ip=0;ip<fMaxPlaneNum;ip++) { Int_t ical=CalibChanInView[fView[ip]] + fException[im][ip]; //if(im==1 && ip>2 && ip<6)ical=ical+50; if(fDCCalib[im][ip][ical]){ Float_t factor=fRefDist/fDCCalib[im][ip][ical]; for(Int_t it=0;it<1000;it++) { fDCCalib[im][ip][it]=fDCCalib[im][ip][it]*factor; } //it }else{ cout<<"can not calibrate plane "<<ip+1<<" in module"<<im+1<<endl; cout<<"check it out!!!"<<endl; } } //ip } //im } //__________________________________________________________________ void BrDriverDC::ReadOffsets(TString OffsetFile, Int_t globalTdcOffset) { // This method read matrix from the disk Int_t im,ip,iw,it,r; ifstream IN(OffsetFile.Data()); cout<<"offset file: "<<OffsetFile<<endl; if(IN){ cout<<"Reading Offsets"<<endl; do { IN>>im>>ip>>iw>>r; if(!(IN.eof()))fOffset[im][ip][iw] = r + globalTdcOffset; // 8 for AuAu above 5620 } while(!(IN.eof())); IN.close(); }else{ cout<<"can not open offset file"<<endl; } } //_____________________________________________________________________ Float_t BrDriverDC::GetDriftDist(Int_t im,Int_t ip,Int_t iw,Int_t il) { Int_t Offset=fOffset[im][ip][iw]; //if(Offset<4600 && Offset>3000) { // corrections for strange warking channels // if(il<4096 && il+512>Offset)il=il+512; //} Int_t it = Offset - il; //cout<<"Offset: "<<Offset<<" "<<it<<endl; if(it>=0 && it<800){ return fDCCalib[im][ip][it]; } else { if(it>-40)return 0.0; else return -100.0; } } //////////////// Methodes for matching global tracks (probably temporary stuff) /////////// //_________________________________________________________________________ void BrDriverDC::SetRunNumber(Int_t irun) { // This function define number of operating plane and detection efficiency // for given run. frunNumber=irun; return; // just define run number if(irun==2535){ fN=18; fEff=0.75; fMeanX=0.0; fSigmaX=1.0; fMeanY=0.71; fSigmaY=0.86; return; } if(irun==2525){ fN=30; fEff=0.30; fMeanX=0.0049; fSigmaX=0.507; fMeanY=0.08; fSigmaY=0.3611; fMeanU=-0.038; fSigmaU=0.421; return; } cout<<"Sorry, you can't follow with this run "<<irun<<endl; } //__________________________________________________________________________ void BrDriverDC::SetMatchCondition(Int_t ik, Float_t ASig) { // Here you define minimum number of planes with // hit that matches a given global track. fK=ik; fASig = ASig; } void BrDriverDC::SetMatchCondition(Float_t fraction) { if(fN){ fK=TMath::Nint(fraction*(Float_t)fN); }else{ cout<<"First you should call BrDriverDC::SetRunNumber"<<endl; } } //___________________________________________________________________________ Float_t BrDriverDC::GetMatchEfficiency() { Int_t n=fN; Double_t p=fEff; Float_t MatchEff=0.0; for(Int_t k=fK;k<=n;k++){ Double_t Cnk = Sil(n) / (Sil(k)*Sil(n-k)); //cout<<"Siln="<<Sil(n)<<" Silkkn="<<Sil(k)*Sil(n-k)<<endl; //cout<<"n="<<n<<" k="<<k<<" p="<<p<<" Cnk="<<Cnk<<endl; MatchEff += Cnk * TMath::Power(p,k) * TMath::Power(1.0-p,n-k); //cout<<TMath::Power(p,k)<<" "<<TMath::Power(1.0-p,n-k)<<" "<<MatchEff<<endl; } Float_t LostEff=0.0; for(Int_t k=0;k<fK;k++){ Double_t Cnk = Sil(n) / (Sil(k)*Sil(n-k)); //cout<<"Siln="<<Sil(n)<<" Silkkn="<<Sil(k)*Sil(n-k)<<endl; //cout<<"n="<<n<<" k="<<k<<" p="<<p<<" Cnk="<<Cnk<<endl; LostEff += Cnk * TMath::Power(p,k) * TMath::Power(1.0-p,n-k); //cout<<TMath::Power(p,k)<<" "<<TMath::Power(1.0-p,n-k)<<" "<<MatchEff<<endl; } cout<<"Match Prob. ="<<MatchEff<<endl; cout<<"Lost Prob. ="<<LostEff<<" Sum="<<MatchEff+LostEff<<endl; return MatchEff; } //____________________________________________________________________________ Double_t BrDriverDC::Sil(Int_t n) { if(n<0){ cout<<"RecProb::Sil: Improper n value"<<endl; return -1; } if(n==0)return 1; Double_t nsil=1; for(Int_t i=1;i<n+1;i++){ nsil=nsil*i; } return nsil; } //____________________________________________________________________________ void BrDriverDC::ResetMatchingStuff() { fNumbOfMatching=0; } //____________________________________________________________________________ Bool_t BrDriverDC::CheckIfMatch(Float_t dev, Int_t iplane) { Bool_t Match=kFALSE; if(iplane==0 || iplane==1 || iplane==2){ // X plane Float_t lowerLimit=fMeanX-fASig*fSigmaX; Float_t upperLimit=fMeanX+fASig*fSigmaX; if(BrMath::BTWN(lowerLimit,dev,upperLimit)) { //Coming here mean that we found a hit consistent with //this global track fNumbOfMatching =+ 1; Match=kTRUE; } } if(iplane==3 || iplane==4 || iplane==5){ // Y plane Float_t lowerLimit=fMeanY-fASig*fSigmaY; Float_t upperLimit=fMeanY+fASig*fSigmaY; if(BrMath::BTWN(lowerLimit,dev,upperLimit)) { //Coming here mean that we found a hit consistent with //this global track fNumbOfMatching += 1; Match=kTRUE; } } if(iplane>5) { // U and V plane Float_t lowerLimit=fMeanU-fASig*fSigmaU; Float_t upperLimit=fMeanU+fASig*fSigmaU; if(BrMath::BTWN(lowerLimit,dev,upperLimit)) { //Coming here mean that we found a hit consistent with //this global track fNumbOfMatching += 1; Match=kTRUE; } } return Match; } //____________________________________________________________________________ Bool_t BrDriverDC::CheckIfMatch(Float_t proj, Int_t imodule, Int_t iplane, Int_t iwire) { Bool_t Match=kFALSE; if(iplane==0 || iplane==1 || iplane==2){ // X plane At(imodule,iplane); Float_t Xpos=GetWirePosition(iwire); Float_t lowerLimit=fMeanX-1.0*fSigmaX; Float_t upperLimit=fMeanX+1.0*fSigmaX; if(BrMath::BTWN(lowerLimit,proj-Xpos,upperLimit)) { //Coming here mean that we found a hit consistent with //this global track fNumbOfMatching =+ 1; Match=kTRUE; } } if(iplane==3 || iplane==4 || iplane==5){ // Y plane At(imodule,iplane); Float_t Ypos=GetWirePosition(iwire); Float_t lowerLimit=fMeanY-1.0*fSigmaY; Float_t upperLimit=fMeanY+1.0*fSigmaY; if(BrMath::BTWN(lowerLimit,proj-Ypos,upperLimit)) { //Coming here mean that we found a hit consistent with //this global track fNumbOfMatching += 1; Match=kTRUE; } } return Match; } //____________________________________________________________________________ Bool_t BrDriverDC::GetMatchingStatus() { if(fNumbOfMatching>=fK){ return kTRUE; }else{ return kFALSE; } } // // $Log: BrDriverDC.cxx,v $ // Revision 1.16 2002/08/08 15:39:32 ufstasze // minor changes // // Revision 1.15 2002/06/10 16:36:15 ufstasze // Added member fRealSetup // // Revision 1.14 2002/05/27 13:06:05 ufstasze // The same exeption added for runs 5542 and 5541 (see ReadCalibData(...) // // Revision 1.13 2002/05/11 14:29:06 ufstasze // Added exception to support run 5544 for T4 // // Revision 1.12 2002/05/02 13:42:46 ufstasze // Added plane pos. correction for T5/II // // Revision 1.11 2002/04/23 15:42:58 videbaek // Remove ios::in and added close/clear calls when re-opening on same ifstream. // This has not been checked. // // Revision 1.10 2002/03/12 09:49:52 radziu // Change in CalibChan in V view from +10 to +12 // // Revision 1.9 2002/02/15 22:19:37 videbaek // Remove useless cout tu1 output // // Revision 1.8 2002/02/15 15:21:10 ufstasze // Modified ReadOffset method by adding globalTdcOffset to the list of arguments // // Revision 1.7 2001/11/12 10:23:12 ufstasze // Now each view have it's own calibration channel value. // // Revision 1.6 2001/08/23 11:31:24 staszel // "Second order" plane position correction for T4 and T5 added plus // minor changes. // // Revision 1.5 2001/08/13 11:48:20 staszel // fView has been added // // Revision 1.4 2001/08/03 15:39:56 staszel // some correction for planes position. // // Revision 1.3 2001/07/30 11:41:41 staszel // Changes that has to do with introduction of T4/T5 calibration // and tracking. // // Revision 1.2 2001/07/20 16:13:15 staszel // Changes to HitIsOK() and At(). // // Revision 1.1.1.1 2001/06/21 14:55:14 hagel // Initial revision of brat2 // // Revision 1.9 2001/06/20 12:17:20 staszel // *** empty log message *** // // Revision 1.8 2001/06/15 15:44:32 staszel // The last revision wasn't last. // Changes having to do with development of dc local tracking package // // Revision 1.7 2001/06/13 13:23:40 staszel // last revision of BrDriverDC // // Revision 1.6 2001/01/19 17:52:36 staszel // Changes having to do with development of BrDCTrackingModule // Few methods that are helpful in confirmation of global track // have been added. // // Revision 1.5 2001/01/04 13:04:12 staszel // Few methods added to allow check if DC hits match global tracks // from - T1,T2 and H1. // // Revision 1.4 2000/11/23 01:33:01 brahmlib // Cleaned up code in general. BrDriverDC now usess ClassImp/ClassDef. // // Revision 1.3 2000/10/25 17:58:06 hagel // Further development of BrDCTrackingModule // // Revision 1.2 2000/10/25 15:01:15 staszel // No we use TMath::Nint and TMath::Abs instead of rint and fabs. // Parameters defined in Init Method. // Minor changes in few methods according to Kris suggestions. // // Revision 1.1 2000/10/11 14:33:48 staszel // First release of BrDriverDC class. It serves geometry of // T3, T4 and T5 detectors. // // |
||||||
|
This page automatically generated by script docBrat by Christian Holm |
Copyright ; 2002 BRAHMS Collaboration
<brahmlib@rcf.rhic.bnl.gov>
|