452 #include "StMagUtilities.h"
457 #include "TGraphErrors.h"
461 #include "StTpcDb/StTpcDb.h"
462 #include "tables/St_MagFactor_Table.h"
463 #include "StDetectorDbMaker/St_tpcHVPlanesC.h"
464 #include "StDetectorDbMaker/St_tpcAnodeHVavgC.h"
465 #include "StDetectorDbMaker/St_tpcHighVoltagesC.h"
466 #include "StDetectorDbMaker/St_tpcFieldCageShortC.h"
467 #include "StDetectorDbMaker/St_tpcOmegaTauC.h"
468 #include "StDetectorDbMaker/St_tpcGridLeakC.h"
469 #include "StDetectorDbMaker/St_spaceChargeCorC.h"
470 #include "StDetectorDbMaker/StTpcSurveyC.h"
471 #include "StDetectorDbMaker/St_trigDetSumsC.h"
472 #include "StDetectorDbMaker/St_tpcPadConfigC.h"
473 #include "StDetectorDbMaker/St_tpcCalibResolutionsC.h"
474 #include "StDetectorDbMaker/St_tpcChargeEventC.h"
475 #include "StDetectorDbMaker/St_tpcSCGLC.h"
476 #include "StDbUtilities/StTpcCoordinateTransform.hh"
478 static Float_t gFactor = 1.0 ;
480 static const Float_t PiOver12 = TMath::Pi()/12. ;
481 static const Float_t PiOver6 = TMath::Pi()/6. ;
482 TNtuple *StMagUtilities::fgDoDistortion = 0;
483 TNtuple *StMagUtilities::fgUnDoDistortion = 0;
484 static const size_t threeFloats = 3 *
sizeof(Float_t);
497 const Float_t GL_charge_y_lo[4] {52.04,121.80-0.85,121.80 ,191.49} ;
498 const Float_t GL_charge_y_hi[4] {52.85,121.80 ,121.80+0.99,192.53} ;
499 Float_t GL_q_inner_of_innerSec(Float_t voltage=1100.0) {
return -42.15 + TMath::Exp(0.005345 * voltage); }
500 Float_t GL_q_outer_of_innerSec(Float_t voltage=1100.0) {
return -21.76 + TMath::Exp(0.005123 * voltage); }
501 Float_t GL_q_inner_of_outerSec(Float_t voltage=1390.0) {
return -16.68 + TMath::Exp(0.004731 * voltage); }
502 Float_t GL_q_outer_of_outerSec(Float_t voltage=1390.0) {
return -23.38 + TMath::Exp(0.004751 * voltage); }
503 Float_t GL_rho_inner_of_innerSec(Float_t voltage=1100.0)
504 {
return GL_q_inner_of_innerSec(voltage) / (GL_charge_y_hi[0]- GL_charge_y_lo[0]) ; }
505 Float_t GL_rho_outer_of_innerSec(Float_t voltage=1100.0)
506 {
return GL_q_outer_of_innerSec(voltage) / (GL_charge_y_hi[1]- GL_charge_y_lo[1]) ; }
507 Float_t GL_rho_inner_of_outerSec(Float_t voltage=1390.0)
508 {
return GL_q_inner_of_outerSec(voltage) / (GL_charge_y_hi[2]- GL_charge_y_lo[2]) ; }
509 Float_t GL_rho_outer_of_outerSec(Float_t voltage=1390.0)
510 {
return GL_q_outer_of_outerSec(voltage) / (GL_charge_y_hi[3]- GL_charge_y_lo[3]) ; }
512 double SpaceChargeRadialDependence(
double Radius) {
513 return ( 3191./(Radius*Radius) + 122.5/Radius - 0.395 ) / 15823. ;
515 double SpaceChargeFXTRadialDependenceEast(
double Radius) {
517 return ( 2427.3/Radius + 42.2021 ) / 1164374.8 ;
519 double SpaceChargeFXTRadialDependenceWest(
double Radius) {
521 return ( 5706.44/Radius + 16.2186 ) / 1173715.5 ;
529 Float_t sector, xL, yL, zL, xLC, yLC, zLC;
532 static const Char_t *Dnames = {
"sector:xL:yL:zL:xLC:yLC:zLC"};
536 void StMagUtilities::SetDoDistortionT (TFile *f) {
539 fgDoDistortion =
new TNtuple(
"DoDist",
"Result of DoDistrotion in TPC CS",Dnames);
542 void StMagUtilities::SetUnDoDistortionT(TFile *f) {
545 fgUnDoDistortion =
new TNtuple(
"UnDoDist",
"Result of UnDoDistrotion in TPC CS",Dnames);
552 cout <<
"ReInstate StMagUtilities. Be sure that this is want you want !" << endl;
553 SafeDelete(fgInstance);
556 GetDistoSmearing(mode);
559 GetTPCVoltages( mode );
565 GetGridLeak( mode ) ;
566 GetAbortGapCharge() ;
567 CommonStart( mode ) ;
568 UseManualSCForPredict(kFALSE) ;
576 cout <<
"ReInstate StMagUtilities. Be sure that this is want you want !" << endl;
577 SafeDelete(fgInstance);
580 GetDistoSmearing(0) ;
585 ManualSpaceCharge(0) ;
586 ManualSpaceChargeR2(0,1) ;
590 CommonStart( mode ) ;
591 UseManualSCForPredict(kFALSE) ;
597 void StMagUtilities::GetDistoSmearing (Int_t mode)
599 fCalibResolutions = ((mode & kDistoSmearing) > 0 ?St_tpcCalibResolutionsC::instance() : 0);
600 mRandom = (fCalibResolutions ?
new TRandom(time(NULL)) : 0);
606 void StMagUtilities::GetMagFactor ()
608 gFactor = StarMagField::Instance()->GetFactor();
610 void StMagUtilities::GetTPCParams ()
616 if (! StTpcDb::IsOldScheme()) {
621 mDistortionMode = kDisableTwistClock;
624 XTWIST = 1e3*glob->TpcEFieldRotationY() ;
625 YTWIST = -1e3*glob->TpcEFieldRotationX() ;
626 EASTCLOCKERROR = 1e3*cages->EastClockError();
627 WESTCLOCKERROR = 1e3*cages->WestClockError();
630 StarDriftV = 1e-6*StTpcDb::instance()->DriftVelocity() ;
631 TPC_Z0 = dims->gatingGridZ() ;
632 IFCShift = cages->InnerFieldCageShift();
633 for (Int_t sec = 1; sec <= 24; sec++) {
634 INNER[sec-1] = pads->innerPadRows(sec);
635 TPCROWS[sec-1] = pads->padRows(sec);
636 for ( Int_t row = 1 ; row <= TPCROWS[sec-1] ; row++ )
637 TPCROWR[sec-1][row-1] = pads->radialDistanceAtRow(sec,row);
640 OFCRadius = dims->senseGasOuterRadius();
641 INNERGGFirst = wires->firstInnerSectorGatingGridWire();
642 OUTERGGFirst = wires->firstOuterSectorGatingGridWire();
643 INNERGGLast = INNERGGFirst +
644 wires->gatingGridWirePitch() * (wires->numInnerSectorGatingGridWires() - 1);
645 OUTERGGLast = OUTERGGFirst +
646 wires->gatingGridWirePitch() * (wires->numOuterSectorGatingGridWires() - 1);
647 GAPRADIUS = 0.5 * (INNERGGLast + OUTERGGFirst);
650 WIREGAP = OUTERGGFirst - INNERGGLast;
653 void StMagUtilities::GetE()
656 Float_t R_0 = 2.130 ;
657 Float_t RStep = 2.000 ;
658 Float_t R_182 = 0.310 ;
659 Rtot = R_0 + 181*RStep + R_182 ;
660 Rfrac = TPC_Z0*RStep/(Rtot*RPitch) ;
664 GGeffectiveness = 0.981;
665 deltaGGeffectiveness = 0.964;
666 GGideal = CathodeV * (1.0 - Rfrac) / GGeffectiveness ;
667 StarMagE = TMath::Abs((CathodeV-GG)/TPC_Z0) ;
668 if (TMath::Abs(StarMagE) < 1e-6) {
669 cout <<
"StMagUtilities ERROR **** Calculations fail with extremely small or zero primary E field:" << endl;
670 cout <<
"StMagUtilities StarMagE = (CathodeV - GG) / TPC_Z0 = (" << CathodeV
671 <<
" - " << GG <<
") / " << TPC_Z0 <<
" = " << StarMagE <<
" V/cm" << endl;
676 void StMagUtilities::GetTPCVoltages (Int_t mode)
679 fTpcVolts = St_tpcHighVoltagesC::instance() ;
680 CathodeV = fTpcVolts->getCathodeVoltage() * 1000 ;
681 GG = fTpcVolts->getGGVoltage() ;
683 if (mode & kPadrow40) {
684 for (Int_t i = 0 ; i < 24; i++ ) {
685 Inner_GLW_Voltage[i] = fTpcVolts->getGridLeakWallTip(i+1) ;
686 Outer_GLW_Voltage[i] = fTpcVolts->getGridLeakWallSide(i+1) ;
688 }
else if (mode & kPadrow13) {
689 for (Int_t i = 0 ; i < 24; i++ ) {
690 Inner_GLW_Voltage[i] = 222 ;
691 Outer_GLW_Voltage[i] = 222 ;
696 if (mode & k3DGridLeak) {
699 double maxInner = 1170;
700 double maxOuter = 1390;
701 double stepsInner = 35;
702 double stepsOuter = 45;
703 TH1I innerVs(
"innerVs",
"innerVs",5,maxInner-3.5*stepsInner,maxInner+1.5*stepsInner);
704 TH1I outerVs(
"outerVs",
"outerVs",5,maxOuter-3.5*stepsOuter,maxOuter+1.5*stepsOuter);
705 for (Int_t i = 1 ; i < 25; i++ ) {
706 int inner = St_tpcPadConfigC::instance()->innerPadRows(i);
707 innerVs.Fill(anodeVolts->voltagePadrow(i,inner));
708 outerVs.Fill(anodeVolts->voltagePadrow(i,inner+1));
710 double cmnInner = innerVs.GetBinCenter(innerVs.GetMaximumBin());
711 double cmnOuter = outerVs.GetBinCenter(outerVs.GetMaximumBin());
712 cout <<
"StMagUtilities assigning common anode voltages as " << cmnInner <<
" , " << cmnOuter << endl;
713 for (Int_t i = 1 ; i < 25; i++ ) {
714 GLWeights[i] = ( ( TMath::Abs(anodeVolts->voltagePadrow(i,INNER[i-1] ) - cmnInner) < stepsInner/2. ) &&
715 ( TMath::Abs(anodeVolts->voltagePadrow(i,INNER[i-1]+1) - cmnOuter) < stepsOuter/2. ) ? 1 : -1 );
717 }
else if (mode & kFullGridLeak) {
721 float norm = ( 122.595*122.595 - 121.0*121.0 ) /
722 ( GL_rho_outer_of_innerSec() *
723 (GL_charge_y_hi[1]*GL_charge_y_hi[1] - GL_charge_y_lo[1]*GL_charge_y_lo[1]) +
724 GL_rho_inner_of_outerSec() *
725 (GL_charge_y_hi[2]*GL_charge_y_hi[2] - GL_charge_y_lo[2]*GL_charge_y_lo[2]) ) ;
727 for (Int_t i = 0 ; i < 24; i++ ) {
728 GLWeights[i ] = GL_rho_inner_of_innerSec(anodeVolts->voltagePadrow(i+1, 1)) * norm ;
729 GLWeights[i+24] = GL_rho_outer_of_innerSec(anodeVolts->voltagePadrow(i+1,INNER[i] )) * norm ;
730 GLWeights[i+48] = GL_rho_inner_of_outerSec(anodeVolts->voltagePadrow(i+1,INNER[i]+1)) * norm ;
731 GLWeights[i+72] = GL_rho_outer_of_outerSec(anodeVolts->voltagePadrow(i+1,TPCROWS[i])) * norm ;
737 Bool_t StMagUtilities::UpdateTPCHighVoltages ()
739 static tpcHighVoltages_st* voltagesTable = 0;
742 tpcHighVoltages_st* new_voltagesTable = voltagesChair->Struct();
743 Bool_t update = (new_voltagesTable != voltagesTable) || DoOnce;
745 if (update) { GetTPCVoltages(kPadrow40); voltagesTable = new_voltagesTable; }
749 void StMagUtilities::GetSpaceCharge ()
751 static spaceChargeCor_st* spaceTable = 0;
755 spaceChargeCor_st* new_spaceTable = spaceChair->Struct();
757 if (new_spaceTable == spaceTable && new_scalers == scalers)
return;
758 fSpaceCharge = spaceChair;
759 spaceTable = new_spaceTable;
760 scalers = new_scalers;
762 SpaceCharge = fSpaceCharge->getSpaceChargeCoulombs((
double)gFactor) ;
765 void StMagUtilities::GetSpaceChargeR2 ()
767 static spaceChargeCor_st* spaceTable = 0;
771 spaceChargeCor_st* new_spaceTable = spaceChair->Struct();
773 if (new_spaceTable == spaceTable && new_scalers == scalers)
return;
774 fSpaceChargeR2 = spaceChair;
775 spaceTable = new_spaceTable;
776 scalers = new_scalers;
778 SpaceChargeR2 = fSpaceChargeR2->getSpaceChargeCoulombs((
double)gFactor) ;
779 SmearCoefSC = (fCalibResolutions ? mRandom->Gaus(1,fCalibResolutions->SpaceCharge()) : 1.0);
780 SpaceChargeEWRatio = fSpaceChargeR2->getEWRatio() ;
783 void StMagUtilities::GetShortedRing ()
786 ShortTableRows = (Int_t) shortedRingsChair->GetNRows() ;
787 for ( Int_t i = 0 ; i < ShortTableRows ; i++)
789 if ( i >= 10 ) break ;
790 Side[i] = (Int_t) shortedRingsChair->side(i) ;
791 Cage[i] = (Int_t) shortedRingsChair->cage(i) ;
792 Ring[i] = (Float_t) shortedRingsChair->ring(i) ;
793 MissingResistance[i] = (Float_t) shortedRingsChair->MissingResistance(i) ;
794 Resistor[i] = (Float_t) shortedRingsChair->resistor(i) ;
798 Bool_t StMagUtilities::UpdateShortedRing ()
800 static tpcFieldCageShort_st* shortsTable = 0;
803 tpcFieldCageShort_st* new_shortsTable = shortedRingsChair->Struct();
804 Bool_t update = (new_shortsTable != shortsTable) || DoOnce;
805 if (update) { GetShortedRing(); shortsTable = new_shortsTable; }
815 Float_t RingNumber, Float_t MissingRValue, Float_t ExtraRValue )
820 Cage[0] = InnerOuter ;
821 Ring[0] = RingNumber ;
822 MissingResistance[0] = MissingRValue ;
823 Resistor[0] = ExtraRValue ;
824 if ( ( Side[0] + Cage[0] + Ring[0] + TMath::Abs(MissingResistance[0]) + TMath::Abs(Resistor[0]) ) == 0.0 ) ShortTableRows = 0 ;
833 void StMagUtilities::GetOmegaTau ()
835 fOmegaTau = St_tpcOmegaTauC::instance();
836 TensorV1 = fOmegaTau->getOmegaTauTensorV1();
837 TensorV2 = fOmegaTau->getOmegaTauTensorV2();
838 mCorrectionsMode = fOmegaTau->distortionCorrectionsMode();
859 if (mDistortionMode & kSpaceCharge) {
860 if (fSpaceCharge)
return 10;
863 if (mDistortionMode & kSpaceChargeR2) {
864 if (fSpaceChargeR2)
return 20;
867 if (mDistortionMode & kSpaceChargeFXT) {
868 if (fSpaceChargeR2)
return 30;
874 void StMagUtilities::GetGridLeak ( Int_t mode )
876 fGridLeak = St_tpcGridLeakC::instance() ;
877 InnerGridLeakStrength = fGridLeak -> getGridLeakStrength ( kGLinner ) ;
878 InnerGridLeakRadius = fGridLeak -> getGridLeakRadius ( kGLinner ) ;
879 InnerGridLeakWidth = fGridLeak -> getGridLeakWidth ( kGLinner ) ;
880 MiddlGridLeakStrength = fGridLeak -> getGridLeakStrength ( kGLmiddl ) ;
881 MiddlGridLeakRadius = fGridLeak -> getGridLeakRadius ( kGLmiddl ) ;
882 MiddlGridLeakWidth = fGridLeak -> getGridLeakWidth ( kGLmiddl ) ;
883 OuterGridLeakStrength = fGridLeak -> getGridLeakStrength ( kGLouter ) ;
884 OuterGridLeakRadius = fGridLeak -> getGridLeakRadius ( kGLouter ) ;
885 OuterGridLeakWidth = fGridLeak -> getGridLeakWidth ( kGLouter ) ;
886 if (mode & kFullGridLeak) {
887 if (InnerGridLeakWidth <= 0) memset( GLWeights ,0,24*
sizeof(Float_t));
888 if (MiddlGridLeakWidth <= 0) memset(&(GLWeights[24]),0,48*
sizeof(Float_t));
889 if (OuterGridLeakWidth <= 0) memset(&(GLWeights[72]),0,24*
sizeof(Float_t));
891 SmearCoefGL = (fCalibResolutions && fCalibResolutions->GridLeak() > 0 ?
892 mRandom->Gaus(1,fCalibResolutions->GridLeak()) : 1.0);
895 void StMagUtilities::ManualGridLeakStrength (Double_t inner, Double_t middle, Double_t outer)
897 InnerGridLeakStrength = inner ;
898 MiddlGridLeakStrength = middle ;
899 OuterGridLeakStrength = outer ;
904 void StMagUtilities::ManualGridLeakRadius (Double_t inner, Double_t middle, Double_t outer)
906 InnerGridLeakRadius = inner ;
907 MiddlGridLeakRadius = middle ;
908 OuterGridLeakRadius = outer ;
911 void StMagUtilities::ManualGridLeakWidth (Double_t inner, Double_t middle, Double_t outer)
913 InnerGridLeakWidth = inner ;
914 MiddlGridLeakWidth = middle ;
915 OuterGridLeakWidth = outer ;
918 void StMagUtilities::GetHVPlanes ()
921 fHVPlanes = St_tpcHVPlanesC::instance() ;
922 tpcHVPlanes_st* HVplanes = fHVPlanes -> Struct();
923 Float_t deltaVGGCathode = GG - GGideal;
924 deltaVGGEast = ( HVplanes -> GGE_shift_z * StarMagE) + deltaVGGCathode;
925 deltaVGGWest = (- HVplanes -> GGW_shift_z * StarMagE) + deltaVGGCathode;
928 void StMagUtilities::ManualGGVoltError (Double_t east, Double_t west)
934 void StMagUtilities::GetAbortGapCharge() {
935 fAbortGapCharge = St_tpcChargeEventC::instance() ;
936 AbortGapCharges = fAbortGapCharge->getCharges() ;
937 AbortGapTimes = fAbortGapCharge->getTimes() ;
938 AbortGapChargeCoef = (St_tpcSCGLC::instance()->SC())[0] ;
939 IonDriftVel = (St_tpcSCGLC::instance()->SC())[1] ;
948 Float_t StMagUtilities::eRList[EMap_nR] = { 48.0, 49.0,
949 50.0, 52.0, 54.0, 56.0, 58.0, 60.0, 62.0, 64.0, 66.0, 68.0,
950 70.0, 72.0, 74.0, 76.0, 78.0, 80.0, 82.0, 84.0, 86.0, 88.0,
951 90.0, 92.0, 94.0, 96.0, 98.0, 100.0, 102.0, 104.0, 106.0, 108.0,
952 110.0, 112.0, 114.0, 116.0, 118.0, 120.0, 122.0, 124.0, 126.0, 128.0,
953 130.0, 132.0, 134.0, 136.0, 138.0, 140.0, 142.0, 144.0, 146.0, 148.0,
954 150.0, 152.0, 154.0, 156.0, 158.0, 160.0, 162.0, 164.0, 166.0, 168.0,
955 170.0, 172.0, 174.0, 176.0, 178.0, 180.0, 182.0, 184.0, 186.0, 188.0,
956 190.0, 192.0, 193.0, 194.0, 195.0, 196.0, 197.0, 198.0, 199.0, 199.5 } ;
958 Float_t StMagUtilities::ePhiList[EMap_nPhi] = { 0.0000, 0.5236, 1.0472, 1.5708, 2.0944, 2.6180, 3.1416,
959 3.6652, 4.1888, 4.7124, 5.2360, 5.7596, 6.2832 } ;
961 Float_t StMagUtilities::eZList[EMap_nZ] = { -208.5, -208.0, -207.0, -206.0, -205.0, -204.0, -202.0,
962 -200.0, -198.0, -196.0, -194.0, -192.0, -190.0, -188.0, -186.0, -184.0, -182.0,
963 -180.0, -178.0, -176.0, -174.0, -172.0, -170.0, -168.0, -166.0, -164.0, -162.0,
964 -160.0, -158.0, -156.0, -154.0, -152.0, -150.0, -148.0, -146.0, -144.0, -142.0,
965 -140.0, -138.0, -136.0, -134.0, -132.0, -130.0, -128.0, -126.0, -124.0, -122.0,
966 -120.0, -118.0, -116.0, -114.0, -112.0, -110.0, -108.0, -106.0, -104.0, -102.0,
967 -100.0, -98.0, -96.0, -94.0, -92.0, -90.0, -88.0, -86.0, -84.0, -82.0,
968 -80.0, -78.0, -76.0, -74.0, -72.0, -70.0, -68.0, -66.0, -64.0, -62.0,
969 -60.0, -58.0, -56.0, -54.0, -52.0, -50.0, -48.0, -46.0, -44.0, -42.0,
970 -40.0, -38.0, -36.0, -34.0, -32.0, -30.0, -28.0, -26.0, -24.0, -22.0,
971 -20.0, -18.0, -16.0, -14.0, -12.0, -10.0, -8.0, -6.0, -4.0, -2.0,
972 -1.0, -0.5, -0.2, -0.1, -0.05, 0.05, 0.1, 0.2, 0.5, 1.0,
973 2.0, 4.0, 6.0, 8.0, 10.0, 12.0, 14.0, 16.0, 18.0, 20.0,
974 22.0, 24.0, 26.0, 28.0, 30.0, 32.0, 34.0, 36.0, 38.0, 40.0,
975 42.0, 44.0, 46.0, 48.0, 50.0, 52.0, 54.0, 56.0, 58.0, 60.0,
976 62.0, 64.0, 66.0, 68.0, 70.0, 72.0, 74.0, 76.0, 78.0, 80.0,
977 82.0, 84.0, 86.0, 88.0, 90.0, 92.0, 94.0, 96.0, 98.0, 100.0,
978 102.0, 104.0, 106.0, 108.0, 110.0, 112.0, 114.0, 116.0, 118.0, 120.0,
979 122.0, 124.0, 126.0, 128.0, 130.0, 132.0, 134.0, 136.0, 138.0, 140.0,
980 142.0, 144.0, 146.0, 148.0, 150.0, 152.0, 154.0, 156.0, 158.0, 160.0,
981 162.0, 164.0, 166.0, 168.0, 170.0, 172.0, 174.0, 176.0, 178.0, 180.0,
982 182.0, 184.0, 186.0, 188.0, 190.0, 192.0, 194.0, 196.0, 198.0, 200.0,
983 202.0, 204.0, 205.0, 206.0, 207.0, 208.0, 208.5 } ;
991 void StMagUtilities::CommonStart ( Int_t mode )
993 cout <<
"StMagUtilities::CommonSta Magnetic Field scale factor is " << gFactor << endl ;
994 if ( StTpcDb::instance() == 0 )
996 cout <<
"StMagUtilities::CommonSta ***NO TPC DB, Using default TPC parameters. You sure it is OK??? ***" << endl ;
997 cout <<
"StMagUtilities::CommonSta ***NO TPC DB, Using default TPC parameters. You sure it is OK??? ***" << endl ;
998 cout <<
"StMagUtilities::CommonSta ***NO TPC DB, Using default TPC parameters. You sure it is OK??? ***" << endl ;
999 cout <<
"StMagUtilities::CommonSta ***NO TPC DB, Using default TPC parameters. You sure it is OK??? ***" << endl ;
1002 #ifndef __NO_TWIST__
1006 IFCShift = 0.0080 ; // Shift of the IFC towards the West Endcap (cm) (2/1/2002)
1007 #ifndef __NO_CLOCK__
1008 EASTCLOCKERROR = 0.0 ;
1009 WESTCLOCKERROR = -0.43 ;
1013 INNERGGFirst = 53.0 ;
1014 INNERGGLast = 121.0 ;
1015 OUTERGGFirst = 122.595 ;
1016 OUTERGGLast = 191.395 ;
1019 for ( Int_t sec = 0; sec < 24; sec++)
1023 for ( Int_t row = 0 ; row < TPCROWS[sec] ; row++ )
1026 TPCROWR[sec][row] = 60.0 + row*4.8 ;
1027 else if ( row < INNER[sec] )
1028 TPCROWR[sec][row] = 93.6 + (row-8+1)*5.2 ;
1030 TPCROWR[sec][row] = 127.195 + (row-INNER[sec])*2.0 ;
1033 mCorrectionsMode = 0;
1035 cout <<
"StMagUtilities::CommonSta WARNING -- Using hard-wired TPC parameters. " << endl ;
1037 else cout <<
"StMagUtilities::CommonSta Using TPC parameters from DataBase. " << endl ;
1039 if ( fTpcVolts == 0 )
1041 CathodeV = -27950.0 ;
1044 for (Int_t i = 0 ; i < 25; i++ ) GLWeights[i] = 1;
1045 for (Int_t i = 0 ; i < 24; i++ ) { Inner_GLW_Voltage[i] = 999; Outer_GLW_Voltage[i] = 999; }
1046 cout <<
"StMagUtilities::CommonSta WARNING -- Using manually selected TpcVoltages setting. " << endl ;
1048 else cout <<
"StMagUtilities::CommonSta Using TPC voltages from the DB." << endl ;
1050 if ( fOmegaTau == 0 )
1054 cout <<
"StMagUtilities::CommonSta WARNING -- Using manually selected OmegaTau parameters. " << endl ;
1056 else cout <<
"StMagUtilities::CommonSta Using OmegaTau parameters from the DB." << endl ;
1058 if (fSpaceCharge) cout <<
"StMagUtilities::CommonSta Using SpaceCharge values from the DB." << endl ;
1059 else cout <<
"StMagUtilities::CommonSta WARNING -- Using manually selected SpaceCharge settings. " << endl ;
1061 if (fSpaceChargeR2) cout <<
"StMagUtilities::CommonSta Using SpaceChargeR2 values from the DB." << endl ;
1062 else cout <<
"StMagUtilities::CommonSta WARNING -- Using manually selected SpaceChargeR2 settings. " << endl ;
1064 if ( fGridLeak == 0 )
1066 ManualGridLeakStrength(0.0,15.0,0.0);
1067 ManualGridLeakRadius(53.0,GAPRADIUS,195.0);
1068 ManualGridLeakWidth(0.0,3.0,0.0);
1069 cout <<
"StMagUtilities::CommonSta WARNING -- Using manually selected GridLeak parameters. " << endl ;
1071 else cout <<
"StMagUtilities::CommonSta Using GridLeak parameters from the DB." << endl ;
1073 if ( fHVPlanes == 0 )
1075 ManualGGVoltError(0.0,0.0);
1076 cout <<
"StMagUtilities::CommonSta WARNING -- Using manually selected HV planes parameters. " << endl ;
1078 else cout <<
"StMagUtilities::CommonSta Using HV planes parameters from the DB." << endl ;
1080 if ( fAbortGapCharge == 0 )
1082 IonDriftVel = 181.67;
1083 cout <<
"StMagUtilities::CommonSta WARNING -- Using default Ion Drift Velocity. " << endl ;
1090 mDistortionMode |= mode;
1091 if (mDistortionMode & kDisableTwistClock) {
1092 mDistortionMode &= ~(mDistortionMode & kTwist);
1093 mDistortionMode &= ~(mDistortionMode & kClock);
1094 mDistortionMode &= ~(mDistortionMode & kFast2DBMap);
1096 if ( !( mode & ( kBMap | kFast2DBMap | kPadrow13 | kPadrow40 | kTwist | kClock | kMembrane | kEndcap |
1097 kIFCShift | kSpaceCharge | kSpaceChargeR2 | kSpaceChargeFXT | kAbortGap | kShortedRing |
1098 kGridLeak | k3DGridLeak | kFullGridLeak | kGGVoltError | kSectorAlign)))
1100 mDistortionMode |= kPadrow40 ;
1101 mDistortionMode |= kAbortGap ;
1102 if (! (mDistortionMode & kDisableTwistClock)) {
1103 mDistortionMode |= kFast2DBMap ;
1104 mDistortionMode |= kTwist ;
1105 mDistortionMode |= kClock ;
1107 mDistortionMode |= kIFCShift ;
1108 printf(
"StMagUtilities::CommonSta Default mode selection\n");
1110 else printf(
"StMagUtilities::CommonSta Using mode option 0x%X\n",mode);
1112 printf(
"StMagUtilities::CommonSta Using correction mode 0x%X\n",mCorrectionsMode);
1113 iterateDistortion = mCorrectionsMode & kIterateUndo;
1114 iterationFailCounter = -1;
1115 printf(
"StMagUtilities::CommonSta Version ");
1116 if ( mDistortionMode & kBMap ) printf (
"3D Mag Field Distortions") ;
1117 if ( mDistortionMode & kFast2DBMap ) printf (
"2D Mag Field Distortions") ;
1118 if ( mDistortionMode & kPadrow13 ) printf (
" + Padrow 13") ;
1119 if ( mDistortionMode & kPadrow40 ) printf (
" + Padrow 40") ;
1120 if ( mDistortionMode & kTwist ) printf (
" + Twist") ;
1121 if ( mDistortionMode & kClock ) printf (
" + Clock") ;
1122 if ( mDistortionMode & kIFCShift ) printf (
" + IFCShift") ;
1123 if ( mDistortionMode & kSpaceCharge ) printf (
" + SpaceCharge") ;
1124 if ( mDistortionMode & kSpaceChargeR2 ) printf (
" + SpaceChargeR2") ;
1125 if ( mDistortionMode & kSpaceChargeFXT) printf (
" + SpaceChargeFXT") ;
1126 if ( mDistortionMode & kAbortGap ) printf (
" + AbortGapSpaceR2") ;
1127 if ( mDistortionMode & kMembrane ) printf (
" + Central Membrane") ;
1128 if ( mDistortionMode & kEndcap ) printf (
" + Endcap") ;
1129 if ( mDistortionMode & kShortedRing ) printf (
" + ShortedRing") ;
1130 if ( mDistortionMode & kGridLeak ) printf (
" + GridLeak") ;
1131 if ( mDistortionMode & k3DGridLeak ) printf (
" + 3DGridLeak") ;
1132 if ( mDistortionMode & kSectorAlign ) printf (
" + SectorAlign") ;
1133 if ( mDistortionMode & kFullGridLeak ) printf (
" + FullGridLeak") ;
1134 if ( mDistortionMode & kDistoSmearing ) printf (
" + DistoSmearing") ;
1135 if ( ! StTpcDb::IsOldScheme()) printf (
" + New TPC Alignment schema") ;
1136 usingCartesian = kTRUE;
1140 Float_t B[3], X[3] = { 0, 0, 0 } ;
1158 OmegaTau = -10.0 * B[2] * StarDriftV / StarMagE ;
1160 Const_0 = 1. / ( 1. + TensorV2*TensorV2*OmegaTau*OmegaTau ) ;
1161 Const_1 = TensorV1*OmegaTau / ( 1. + TensorV1*TensorV1*OmegaTau*OmegaTau ) ;
1162 Const_2 = TensorV2*TensorV2*OmegaTau*OmegaTau / ( 1. + TensorV2*TensorV2*OmegaTau*OmegaTau ) ;
1164 cout <<
"StMagUtilities::BField = " << B[2] <<
" kGauss at (0,0,0)" << endl ;
1165 cout <<
"StMagUtilities::DriftVel = " << StarDriftV <<
" cm/microsec" << endl ;
1166 cout <<
"StMagUtilities::TPC_Z0 = " << TPC_Z0 <<
" cm" << endl ;
1167 cout <<
"StMagUtilities::TensorV1+V2 = " << TensorV1 <<
" " << TensorV2 << endl ;
1168 cout <<
"StMagUtilities::OmegaTau1+2 = " << OmegaTau * TensorV1 <<
" " << OmegaTau * TensorV2 << endl ;
1169 cout <<
"StMagUtilities::XTWIST = " << XTWIST <<
" mrad" << endl ;
1170 cout <<
"StMagUtilities::YTWIST = " << YTWIST <<
" mrad" << endl ;
1171 cout <<
"StMagUtilities::SpaceCharge = " << SpaceCharge <<
" Coulombs/epsilon-nought" << endl ;
1172 cout <<
"StMagUtilities::SpaceChargeR2 = " << SpaceChargeR2 <<
" Coulombs/epsilon-nought" <<
" EWRatio = "
1173 << SpaceChargeEWRatio << endl ;
1174 if (mDistortionMode & kDistoSmearing) {
1175 cout <<
"StMagUtilities::SmearCoefSC = " << SmearCoefSC << endl;
1176 cout <<
"StMagUtilities::SmearCoefGL = " << SmearCoefGL << endl;
1178 if (mDistortionMode & kAbortGap) {
1179 cout <<
"StMagUtilities::AbortGapCoef = " << AbortGapChargeCoef << endl;
1181 cout <<
"StMagUtilities::IFCShift = " << IFCShift <<
" cm" << endl ;
1182 cout <<
"StMagUtilities::CathodeV = " << CathodeV <<
" volts" << endl ;
1183 cout <<
"StMagUtilities::GG = " << GG <<
" volts" << endl ;
1184 if (mDistortionMode & kPadrow40) {
1185 cout <<
"StMagUtilities::Inner_GLW_V = " ;
1186 for ( Int_t i = 0 ; i < 24 ; i++ ) cout << Inner_GLW_Voltage[i] <<
" " ; cout <<
"volts" << endl ;
1187 cout <<
"StMagUtilities::Outer_GLW_V = " ;
1188 for ( Int_t i = 0 ; i < 24 ; i++ ) cout << Outer_GLW_Voltage[i] <<
" " ; cout <<
"volts" << endl ;
1190 cout <<
"StMagUtilities::EastClock = " << EASTCLOCKERROR <<
" mrad" << endl;
1191 cout <<
"StMagUtilities::WestClock = " << WESTCLOCKERROR <<
" mrad" << endl;
1192 cout <<
"StMagUtilities::Side = " ;
1193 for ( Int_t i = 0 ; i < ShortTableRows ; i++ ) cout << Side[i] <<
" " ; cout <<
"Location of Short E=0 / W=1 " << endl;
1194 cout <<
"StMagUtilities::Cage = " ;
1195 for ( Int_t i = 0 ; i < ShortTableRows ; i++ ) cout << Cage[i] <<
" " ; cout <<
"Location of Short IFC = 0 / OFC = 1" << endl;
1196 cout <<
"StMagUtilities::Ring = " ;
1197 for ( Int_t i = 0 ; i < ShortTableRows ; i++ ) cout << Ring[i] <<
" " ; cout <<
"Rings - Location of Short counting from the CM" << endl;
1198 cout <<
"StMagUtilities::MissingOhms = " ;
1199 for ( Int_t i = 0 ; i < ShortTableRows ; i++ ) cout << MissingResistance[i] <<
" " ; cout <<
"MOhms Missing Resistance" << endl;
1200 cout <<
"StMagUtilities::CompResistor = " ;
1201 for ( Int_t i = 0 ; i < ShortTableRows ; i++ ) cout << Resistor[i] <<
" " ; cout <<
"MOhm Compensating Resistor Value" << endl;
1202 cout <<
"StMagUtilities::InnerGridLeak = " << InnerGridLeakStrength <<
" " << InnerGridLeakRadius <<
" " << InnerGridLeakWidth << endl;
1203 cout <<
"StMagUtilities::MiddlGridLeak = " << MiddlGridLeakStrength <<
" " << MiddlGridLeakRadius <<
" " << MiddlGridLeakWidth << endl;
1204 cout <<
"StMagUtilities::OuterGridLeak = " << OuterGridLeakStrength <<
" " << OuterGridLeakRadius <<
" " << OuterGridLeakWidth << endl;
1205 cout <<
"StMagUtilities::deltaVGG = " << deltaVGGEast <<
" V (east) : " << deltaVGGWest <<
" V (west)" << endl;
1206 cout <<
"StMagUtilities::GLWeights = " ;
1207 if (mDistortionMode & k3DGridLeak) {
1208 for ( Int_t i = 1 ; i < 25 ; i++ ) cout << GLWeights[i] <<
" " ;
1210 }
else if (mDistortionMode & kFullGridLeak) {
1212 for ( Int_t i = 0 ; i < 24 ; i++ ) {
1213 for ( Int_t j = 0 ; j < 96 ; j+=24 ) cout << std::fixed <<
" " << std::setprecision(2) << GLWeights[i+j];
1216 }
else cout <<
"N/A" << endl;
1218 doingDistortion = kFALSE;
1230 void StMagUtilities::UndoDistortion(
const Float_t x[], Float_t Xprime[] , Int_t Sector )
1239 const Float_t Xtemp1[3] = {100.,0.,100.};
1241 Bool_t tempIterDist = iterateDistortion;
1242 iterateDistortion = kFALSE;
1243 UndoDistortion(Xtemp1,Xtemp2,3);
1244 iterateDistortion = tempIterDist;
1248 Float_t Xprime1[3], Xprime2[3] ;
1250 SectorNumber( Sector, x ) ;
1252 if (iterateDistortion) {
1254 iterateDistortion = kFALSE;
1256 memcpy(Xprime1,x,threeFloats);
1257 memcpy(Xprime,x,threeFloats);
1258 const Float_t MINDIST_SQ = 1e-8;
1259 Float_t dist_sq = 1e10;
1261 while (dist_sq > MINDIST_SQ) {
1262 UndoDistortion ( Xprime1, Xprime2, Sector ) ;
1263 Xprime1[0] = x[0] - (Xprime1[0] - Xprime2[0]) ;
1264 Xprime1[1] = x[1] - (Xprime1[1] - Xprime2[1]) ;
1265 Xprime1[2] = x[2] - (Xprime1[2] - Xprime2[2]) ;
1266 dist_sq = (Xprime1[0]-Xprime[0])*(Xprime1[0]-Xprime[0])
1267 + (Xprime1[1]-Xprime[1])*(Xprime1[1]-Xprime[1])
1268 + (Xprime1[2]-Xprime[2])*(Xprime1[2]-Xprime[2]) ;
1269 if (++iter > 20 && dist_sq > MINDIST_SQ) {
1271 Xprime[0] = 0.5*(Xprime[0] + Xprime1[0]);
1272 Xprime[1] = 0.5*(Xprime[1] + Xprime1[1]);
1273 Xprime[2] = 0.5*(Xprime[2] + Xprime1[2]);
1275 if (iterationFailCounter >=0) iterationFailCounter++;
1277 memcpy(Xprime,Xprime1,threeFloats);
1281 iterateDistortion = kTRUE;
1298 if ( TMath::Abs(x[2]) >= TPC_Z0 )
1299 { memcpy(Xprime,x,threeFloats); return ; }
1300 Cart2Polar(x,Xprime1[0],Xprime1[1]);
1301 if (Xprime1[0] > OFCRadius || Xprime1[0] < IFCRadius )
1302 { memcpy(Xprime,x,threeFloats); return ; }
1304 usingCartesian = kFALSE;
1306 if (mDistortionMode & kBMap) {
1308 memcpy(Xprime1,Xprime2,threeFloats);
1311 if (mDistortionMode & kFast2DBMap) {
1313 memcpy(Xprime1,Xprime2,threeFloats);
1316 if ((mDistortionMode & kBMap) && (mDistortionMode & kFast2DBMap)) {
1317 cout <<
"StMagUtilities ERROR **** Do not use kBMap and kFast2DBMap at the same time" << endl ;
1318 cout <<
"StMagUtilities ERROR **** These routines have duplicate functionality so don't do both." << endl ;
1322 if ((mDistortionMode & kPadrow13) && (mDistortionMode & kPadrow40)) {
1323 cout <<
"StMagUtilities ERROR **** Do not use kPadrow13 and kPadrow40 at the same time" << endl ;
1324 cout <<
"StMagUtilities ERROR **** These routines have duplicate functionality so don't do both." << endl ;
1328 if (mDistortionMode & kPadrow13) {
1330 memcpy(Xprime1,Xprime2,threeFloats);
1333 if (mDistortionMode & kPadrow40) {
1335 memcpy(Xprime1,Xprime2,threeFloats);
1338 if (mDistortionMode & kTwist) {
1340 memcpy(Xprime1,Xprime2,threeFloats);
1343 if (mDistortionMode & kClock) {
1345 memcpy(Xprime1,Xprime2,threeFloats);
1348 if (mDistortionMode & kMembrane) {
1350 memcpy(Xprime1,Xprime2,threeFloats);
1353 if (mDistortionMode & kEndcap) {
1355 memcpy(Xprime1,Xprime2,threeFloats);
1358 if (mDistortionMode & kIFCShift) {
1360 memcpy(Xprime1,Xprime2,threeFloats);
1363 if (mDistortionMode & (kSpaceCharge | kSpaceChargeR2 | kSpaceChargeFXT)) {
1365 memcpy(Xprime1,Xprime2,threeFloats);
1368 if (mDistortionMode & kAbortGap) {
1370 memcpy(Xprime1,Xprime2,threeFloats);
1373 if (mDistortionMode & kShortedRing) {
1375 memcpy(Xprime1,Xprime2,threeFloats);
1378 if (mDistortionMode & (kGridLeak | k3DGridLeak | kFullGridLeak)) {
1380 memcpy(Xprime1,Xprime2,threeFloats);
1383 if (mDistortionMode & kGGVoltError) {
1385 memcpy(Xprime1,Xprime2,threeFloats);
1388 if (mDistortionMode & kSectorAlign) {
1390 memcpy(Xprime1,Xprime2,threeFloats);
1397 Polar2Cart(Xprime1[0],Xprime1[1],Xprime);
1398 Xprime[2] = Xprime1[2];
1399 usingCartesian = kTRUE;
1402 if (! iterateDistortion && fgUnDoDistortion) {
1410 fgUnDoDistortion->Fill(&D.sector);
1423 Bool_t tempIterDist = iterateDistortion;
1424 Bool_t tempDoingDist = doingDistortion;
1425 iterateDistortion = kFALSE;
1426 doingDistortion = kTRUE;
1428 UndoDistortion ( x, Xprime, Sector ) ;
1430 Xprime[0] = 2*x[0] - Xprime[0] ;
1431 Xprime[1] = 2*x[1] - Xprime[1] ;
1432 Xprime[2] = 2*x[2] - Xprime[2] ;
1434 iterateDistortion = tempIterDist;
1436 if (fgDoDistortion) {
1444 fgDoDistortion->Fill(&D.sector);
1446 doingDistortion = tempDoingDist;
1466 Int_t sign, index = 1 , NSTEPS ;
1468 sign = SectorSide( Sector, x ) ;
1472 Xprime[2] = sign * TPC_Z0 ;
1474 for ( NSTEPS = 5 ; NSTEPS < 1000 ; NSTEPS += 2 )
1476 ah = ( x[2] - sign * TPC_Z0 ) / ( NSTEPS - 1 ) ;
1477 if ( TMath::Abs(ah) < 1.0 )
break ;
1480 for ( Int_t i = 1; i <= NSTEPS; ++i )
1482 if ( i == NSTEPS ) index = 1 ;
1483 Xprime[2] += index*(ah/3) ;
1484 B3DFieldTpc( Xprime, B , Sector) ;
1485 if ( TMath::Abs(B[2]) > 0.001 )
1487 Xprime[0] += index*(ah/3)*( Const_2*B[0] - Const_1*B[1] ) / B[2] ;
1488 Xprime[1] += index*(ah/3)*( Const_2*B[1] + Const_1*B[0] ) / B[2] ;
1490 if ( index != 4 ) index = 4;
else index = 2 ;
1508 Int_t sign, index = 1 , NSTEPS ;
1510 sign = SectorSide( Sector, x ) ;
1514 Xprime[2] = sign * TPC_Z0 ;
1516 for ( NSTEPS = 5 ; NSTEPS < 1000 ; NSTEPS += 2 )
1518 ah = ( x[2] - sign * TPC_Z0 ) / ( NSTEPS - 1 ) ;
1519 if ( TMath::Abs(ah) < 1.0 )
break ;
1522 for ( Int_t i = 1; i <= NSTEPS; ++i )
1524 if ( i == NSTEPS ) index = 1 ;
1525 Xprime[2] += index*(ah/3) ;
1526 BFieldTpc( Xprime, B, Sector);
1527 if ( TMath::Abs(B[2]) > 0.001 )
1529 Xprime[0] += index*(ah/3)*( Const_2*B[0] - Const_1*B[1] ) / B[2] ;
1530 Xprime[1] += index*(ah/3)*( Const_2*B[1] + Const_1*B[0] ) / B[2] ;
1532 if ( index != 4 ) index = 4;
else index = 2 ;
1548 static Float_t dx3D[EMap_nPhi][EMap_nR][EMap_nZ], dy3D[EMap_nPhi][EMap_nR][EMap_nZ] ;
1549 static Int_t ilow = 0, jlow = 0, klow = 0 ;
1550 const Int_t PHIORDER = 1 ;
1551 const Int_t ORDER = 1 ;
1555 Float_t save_x[ORDER+1], saved_x[ORDER+1] ;
1556 Float_t save_y[ORDER+1], saved_y[ORDER+1] ;
1559 if (usingCartesian) Cart2Polar(x,r,phi);
1560 else { r = x[0]; phi = x[1]; }
1561 if ( phi < 0 ) phi += TMath::TwoPi() ;
1562 Float_t z = LimitZ( Sector, x ) ;
1566 cout <<
"StMagUtilities::FastUndoD Please wait for the tables to fill ... ~90 seconds" << endl ;
1567 for ( k = 0 ; k < EMap_nPhi ; k++ )
1569 for ( i = 0 ; i < EMap_nR ; i++ )
1571 xx[0] = eRList[i] * TMath::Cos(ePhiList[k]) ;
1572 xx[1] = eRList[i] * TMath::Sin(ePhiList[k]) ;
1573 for ( j = 0 ; j < EMap_nZ ; j++ )
1577 dx3D[k][i][j] = Xprime[0] - xx[0] ;
1578 dy3D[k][i][j] = Xprime[1] - xx[1] ;
1584 Search( EMap_nPhi, ePhiList, phi, klow ) ;
1585 Search( EMap_nR, eRList, r, ilow ) ;
1586 Search( EMap_nZ, eZList, z, jlow ) ;
1587 if ( klow < 0 ) klow = 0 ;
1588 if ( ilow < 0 ) ilow = 0 ;
1589 if ( jlow < 0 ) jlow = 0 ;
1590 if ( klow + ORDER >= EMap_nPhi-1 ) klow = EMap_nPhi - 1 - ORDER ;
1591 if ( ilow + ORDER >= EMap_nR-1 ) ilow = EMap_nR - 1 - ORDER ;
1592 if ( jlow + ORDER >= EMap_nZ-1 ) jlow = EMap_nZ - 1 - ORDER ;
1594 for ( k = klow ; k < klow + ORDER + 1 ; k++ )
1596 for ( i = ilow ; i < ilow + ORDER + 1 ; i++ )
1598 save_x[i-ilow] = Interpolate( &eZList[jlow], &dx3D[k][i][jlow], ORDER, z ) ;
1599 save_y[i-ilow] = Interpolate( &eZList[jlow], &dy3D[k][i][jlow], ORDER, z ) ;
1601 saved_x[k-klow] = Interpolate( &eRList[ilow], save_x, ORDER, r ) ;
1602 saved_y[k-klow] = Interpolate( &eRList[ilow], save_y, ORDER, r ) ;
1605 if (usingCartesian) {
1606 Xprime[0] = Interpolate( &ePhiList[klow], saved_x, PHIORDER, phi ) + x[0] ;
1607 Xprime[1] = Interpolate( &ePhiList[klow], saved_y, PHIORDER, phi ) + x[1];
1609 Polar2Cart(x[0],x[1],xx);
1610 xx[0] += Interpolate( &ePhiList[klow], saved_x, PHIORDER, phi ) ;
1611 xx[1] += Interpolate( &ePhiList[klow], saved_y, PHIORDER, phi ) ;
1612 Cart2Polar(xx,Xprime[0],Xprime[1]);
1632 static Float_t dR[EMap_nR][EMap_nZ], dRPhi[EMap_nR][EMap_nZ] ;
1633 static Int_t ilow = 0, jlow = 0 ;
1634 const Int_t ORDER = 1 ;
1638 Float_t save_dR[ORDER+1], saved_dR ;
1639 Float_t save_dRPhi[ORDER+1], saved_dRPhi ;
1642 if (usingCartesian) Cart2Polar(x,r,phi);
1643 else { r = x[0]; phi = x[1]; }
1644 if ( phi < 0 ) phi += TMath::TwoPi() ;
1645 Float_t z = LimitZ( Sector, x ) ;
1649 cout <<
"StMagUtilities::FastUndo2 Please wait for the tables to fill ... ~5 seconds" << endl ;
1650 for ( i = 0 ; i < EMap_nR ; i++ )
1654 for ( j = 0 ; j < EMap_nZ ; j++ )
1658 dR[i][j] = Xprime[0] ;
1659 dRPhi[i][j] = Xprime[1] ;
1664 Search( EMap_nR, eRList, r, ilow ) ;
1665 Search( EMap_nZ, eZList, z, jlow ) ;
1666 if ( ilow < 0 ) ilow = 0 ;
1667 if ( jlow < 0 ) jlow = 0 ;
1668 if ( ilow + ORDER >= EMap_nR-1 ) ilow = EMap_nR - 1 - ORDER ;
1669 if ( jlow + ORDER >= EMap_nZ-1 ) jlow = EMap_nZ - 1 - ORDER ;
1671 for ( i = ilow ; i < ilow + ORDER + 1 ; i++ )
1673 save_dR[i-ilow] = Interpolate( &eZList[jlow], &dR[i][jlow], ORDER, z ) ;
1674 save_dRPhi[i-ilow] = Interpolate( &eZList[jlow], &dRPhi[i][jlow], ORDER, z ) ;
1677 saved_dR = Interpolate( &eRList[ilow], save_dR, ORDER, r ) ;
1678 saved_dRPhi = Interpolate( &eRList[ilow], save_dRPhi, ORDER, r ) ;
1684 phi = phi + saved_dRPhi / r ;
1685 if ( phi < 0 ) phi += 2*TMath::Pi() ;
1688 if (usingCartesian) Polar2Cart(r,phi,Xprime);
1689 else { Xprime[0] = r; Xprime[1] = phi; }
1713 Float_t z = LimitZ( Sector, x ) ;
1714 sign = SectorSide( Sector, x ) ;
1716 Zdrift = sign * ( TPC_Z0 - TMath::Abs(z) ) ;
1717 if (usingCartesian) {
1718 Xprime[0] = x[0] - ( Const_1 * YTWIST - Const_2 * XTWIST ) * Zdrift/1000 ;
1719 Xprime[1] = x[1] - ( -1* Const_1 * XTWIST - Const_2 * YTWIST ) * Zdrift/1000 ;
1722 Polar2Cart(x[0],x[1],xx);
1723 xx[0] -= ( Const_1 * YTWIST - Const_2 * XTWIST ) * Zdrift/1000 ;
1724 xx[1] -= ( -1* Const_1 * XTWIST - Const_2 * YTWIST ) * Zdrift/1000 ;
1725 Cart2Polar(xx,Xprime[0],Xprime[1]);
1745 const Int_t ORDER = 2 ;
1746 const Int_t NZDRIFT = 19 ;
1747 const Int_t NYARRAY = 37 ;
1748 const Int_t TERMS = 400 ;
1749 const Float_t SCALE = 0.192 ;
1750 const Float_t BOX = 200.0 - GAPRADIUS ;
1751 const Float_t PI = TMath::Pi() ;
1757 static Float_t ZDriftArray[NZDRIFT] = {0,1,2,3,4,5,7.5,10,12.5,15,17.5,20,25,30,50,75,100,210,220} ;
1759 static Float_t YArray[NYARRAY] = { 50.0, 75.0, 100.0,
1760 103.5, 104.0, 104.5,
1761 108.7, 109.2, 109.7,
1762 113.9, 114.4, 114.9,
1763 118.9, 119.6, 119.8,
1764 120.0, 120.25, 120.5, 120.75,
1765 121.0, 121.5, 122.1, 122.6,
1774 static Double_t C[TERMS] ;
1775 static Float_t SumArray[NZDRIFT][NYARRAY] ;
1776 static Int_t ilow = 0, jlow = 0 ;
1778 Float_t y, z, Zdrift, save_sum[3] ;
1779 Double_t r, phi, phi0, sum = 0.0 ;
1783 cout <<
"StMagUtilities::PadRow13 Please wait for the tables to fill ... ~5 seconds" << endl ;
1784 C[0] = WIREGAP * GG * SCALE / ( 2 * BOX ) ;
1785 for ( Int_t i = 1 ; i < TERMS ; i++ )
1786 C[i] = 2 * GG * SCALE * TMath::Sin( WIREGAP*i*PI/( 2*BOX ) ) / ( i * PI ) ;
1787 for ( Int_t i = 0; i < NZDRIFT ; i++ )
1789 Zdrift = ZDriftArray[i] ;
1790 for ( Int_t j = 0; j < NYARRAY ; j++ )
1794 for ( Int_t k = 1 ; k < TERMS ; k++ )
1796 sum += ( C[k] / StarMagE ) * ( 1. - TMath::Exp(-1*k*PI*Zdrift/BOX) )
1797 * TMath::Sin(k*PI*(y-GAPRADIUS)/BOX) ;
1799 SumArray[i][j] = sum ;
1804 if (usingCartesian) Cart2Polar(x,r,phi);
1805 else { r = x[0]; phi = x[1]; }
1806 phi0 = ( (Int_t)((TMath::Abs(phi)+PiOver12)/PiOver6 + 6.0 ) - 6.0 ) * PiOver6 ;
1807 if ( phi < 0 ) phi0 *= -1. ;
1808 y = r * TMath::Cos( phi0 - phi ) ;
1809 z = LimitZ( Sector, x ) ;
1810 Zdrift = TPC_Z0 - TMath::Abs(z) ;
1812 Search ( NZDRIFT, ZDriftArray, Zdrift, ilow ) ;
1813 Search ( NYARRAY, YArray, y, jlow ) ;
1815 if ( ilow < 0 ) ilow = 0 ;
1816 if ( jlow < 0 ) jlow = 0 ;
1817 if ( ilow + ORDER >= NZDRIFT - 1 ) ilow = NZDRIFT - 1 - ORDER ;
1818 if ( jlow + ORDER >= NYARRAY - 1 ) jlow = NYARRAY - 1 - ORDER ;
1820 for ( Int_t i = ilow ; i < ilow + ORDER + 1 ; i++ )
1822 save_sum[i-ilow] = Interpolate( &YArray[jlow], &SumArray[i][jlow], ORDER, y ) ;
1825 sum = Interpolate( &ZDriftArray[ilow], save_sum, ORDER, Zdrift ) ;
1829 phi = phi - ( Const_1*(-1*sum)*TMath::Cos(phi0-phi) + Const_0*sum*TMath::Sin(phi0-phi) ) / r ;
1830 r = r - ( Const_0*sum*TMath::Cos(phi0-phi) - Const_1*(-1*sum)*TMath::Sin(phi0-phi) ) ;
1833 if (usingCartesian) Polar2Cart(r,phi,Xprime);
1834 else { Xprime[0] = r; Xprime[1] = phi; }
1888 const Int_t ROWS = 5401 ;
1889 const Int_t COLUMNS = 10001 ;
1890 const Int_t SAVEDCOLUMNS = 1001 ;
1891 const Int_t GPPMM = 20 ;
1892 const Int_t nMAPS = 4 ;
1893 const Int_t nTERMS = 800 ;
1894 const Int_t nZDRIFT = 20 ;
1895 const Int_t nYARRAY = 61 ;
1896 const Int_t ORDER = 2 ;
1897 const Float_t BOX = (COLUMNS-1)/(GPPMM*10.0) ;
1898 const Float_t PI = TMath::Pi() ;
1905 static Int_t ilow = 0, jlow = 0 ;
1906 static Float_t Bn[nTERMS+1] ;
1907 static Float_t SumArray[nMAPS][nZDRIFT][nYARRAY] ;
1908 static Float_t ZDriftArray[nZDRIFT] = {0,1,2,3,4,5,7.5,10,12.5,15,17.5,20,22.5,25,30,50,75,100,210,220} ;
1910 static Float_t YArray[nYARRAY] = { 50.0, 75.0, 100.0, 102.0,
1911 103.5, 104.0, 104.5, 105.0,
1912 105.8, 106.6, 107.4, 108.2,
1913 108.7, 109.2, 109.7, 110.15,
1914 110.95, 111.75, 112.5, 113.3,
1915 113.9, 114.4, 114.9, 115.75,
1916 116.55, 117.35, 118.15, 118.9,
1917 119.25, 119.6, 119.8, 120.0,
1918 120.25, 120.5, 120.75, 121.0,
1919 121.5, 122.1, 122.6, 124.2,
1920 125.2, 126.2, 127.195, 128.2,
1921 129.2, 130.2, 131.195, 132.2,
1922 133.2, 134.2, 135.195, 136.2,
1923 137.2, 139.2, 140.2, 142.2,
1924 144.2, 146.2, 150.0, 198.,
1927 Float_t r, phi, phi0, totalSum ;
1928 Float_t cosPhi0Phi, sinPhi0Phi ;
1929 Float_t y, z, save_sum[nMAPS][3] ;
1930 Float_t Zdrift, sum ;
1931 Float_t MapSum[nMAPS] ;
1932 Float_t DataInTheGap[SAVEDCOLUMNS] ;
1933 static Bool_t DoOnceLocal = true ;
1935 if (fTpcVolts) {DoOnceLocal = UpdateTPCHighVoltages() ;}
1938 cout <<
"StMagUtilities::PadRow40 Filling tables ..." << endl ;
1939 Int_t OFFSET = (COLUMNS-SAVEDCOLUMNS)/2 ;
1941 for ( Int_t MapID = 0 ; MapID < nMAPS ; MapID++ )
1944 GetGLWallData ( MapID, DataInTheGap ) ;
1945 for ( Int_t n = 1 ; n <= nTERMS ; n++ )
1947 Float_t COEFFICIENT = n*PI/(COLUMNS-1) ;
1948 sum = DataInTheGap[0] * TMath::Sin(COEFFICIENT*OFFSET) ;
1949 for ( Int_t i = 1 ; i < SAVEDCOLUMNS ; i+=2 ) sum += 4 * DataInTheGap[i] * TMath::Sin(COEFFICIENT*(i+OFFSET)) ;
1950 for ( Int_t i = 2 ; i < SAVEDCOLUMNS ; i+=2 ) sum += 2 * DataInTheGap[i] * TMath::Sin(COEFFICIENT*(i+OFFSET)) ;
1951 sum -= DataInTheGap[SAVEDCOLUMNS-1]*TMath::Sin(COEFFICIENT*(SAVEDCOLUMNS-1+OFFSET)) ;
1952 sum *= 2.0/(3.0*(COLUMNS-1)) ;
1956 for ( Int_t i = 0 ; i < nZDRIFT ; i++ )
1958 z = ZDriftArray[i] ;
1959 if ( z > (ROWS-1)/(10.0*GPPMM) ) z = (ROWS-1)/(10.0*GPPMM) ;
1960 for ( Int_t j = 0; j < nYARRAY ; j++ )
1964 SumArray[MapID][i][j] = 0.0 ;
1965 if ( y <= GAPRADIUS-BOX/2 || y >= GAPRADIUS+BOX/2 ) continue ;
1966 for ( Int_t n = 1 ; n <= nTERMS ; n++ )
1968 sum += ( Bn[n]/StarMagE ) * ( 1. - TMath::Exp((-1*n*PI*z)/BOX) ) * TMath::Cos((n*PI*(y-GAPRADIUS+BOX/2))/BOX) ;
1970 SumArray[MapID][i][j] = sum ;
1975 DoOnceLocal = false ;
1978 if (usingCartesian) Cart2Polar(x,r,phi) ;
1979 else { r = x[0] ; phi = x[1] ; }
1980 phi0 = ( (Int_t)((TMath::Abs(phi)+PiOver12)/PiOver6 + 6.0 ) - 6.0 ) * PiOver6 ;
1981 if ( phi < 0 ) phi0 *= -1. ;
1982 cosPhi0Phi = TMath::Cos( phi0 - phi ) ;
1983 sinPhi0Phi = TMath::Sin( phi0 - phi ) ;
1984 y = r * cosPhi0Phi ;
1985 z = LimitZ( Sector, x ) ;
1986 Zdrift = TPC_Z0 - TMath::Abs(z) ;
1988 Search ( nZDRIFT, ZDriftArray, Zdrift, ilow ) ;
1989 Search ( nYARRAY, YArray, y, jlow ) ;
1991 if ( ilow < 0 ) ilow = 0 ;
1992 if ( jlow < 0 ) jlow = 0 ;
1993 if ( ilow + ORDER >= nZDRIFT - 1 ) ilow = nZDRIFT - 1 - ORDER ;
1994 if ( jlow + ORDER >= nYARRAY - 1 ) jlow = nYARRAY - 1 - ORDER ;
1996 for ( Int_t MapID = 0 ; MapID < nMAPS ; MapID++ )
1998 for ( Int_t i = ilow ; i < ilow + ORDER + 1 ; i++ )
2000 save_sum[MapID][i-ilow] = Interpolate( &YArray[jlow], &SumArray[MapID][i][jlow], ORDER, y ) ;
2002 MapSum[MapID] = Interpolate( &ZDriftArray[ilow], save_sum[MapID], ORDER, Zdrift ) ;
2005 if ( Inner_GLW_Voltage[Sector-1] >= 100.0 ) totalSum = MapSum[0] ;
2006 else totalSum = MapSum[3] - ((Outer_GLW_Voltage[Sector-1]+113.0)/100.0)*MapSum[1] + ((Inner_GLW_Voltage[Sector-1]+113.0)/100.0)*MapSum[2] ;
2007 totalSum = -totalSum;
2011 phi = phi - ( Const_1*(-1*totalSum)*cosPhi0Phi + Const_0*totalSum*sinPhi0Phi ) / r ;
2012 r = r - ( Const_0*totalSum*cosPhi0Phi - Const_1*(-1*totalSum)*sinPhi0Phi ) ;
2015 if (usingCartesian) Polar2Cart(r,phi,Xprime) ;
2016 else { Xprime[0] = r ; Xprime[1] = phi ; }
2021 void StMagUtilities::GetGLWallData (
const Int_t select, Float_t DataInTheGap[] )
2048 Float_t output_nowall_oldTPC[1001] =
2049 { -0.0399262, -0.0363446, -0.0319311, -0.026963, -0.0217547, -0.0166374, -0.0119395, -0.00796491, -0.00497478, -0.00317029,
2050 -0.0026798, -0.00355035, -0.00574406, -0.00913993, -0.0135409, -0.0186854, -0.0242637, -0.0299365, -0.0353564, -0.0401892,
2051 -0.0441351, -0.0469478, -0.0484502, -0.0485461, -0.0472267, -0.0445716, -0.0407452, -0.0359861, -0.0305932, -0.0249076,
2052 -0.0192911, -0.0141041, -0.0096827, -0.0063177, -0.00423598, -0.00358583, -0.00442701, -0.00672641, -0.0103594, -0.0151168,
2053 -0.0207176, -0.0268257, -0.0330708, -0.0390711, -0.0444573, -0.0488955, -0.0521087, -0.0538948, -0.054139, -0.0528226,
2054 -0.0500243, -0.0459166, -0.0407551, -0.0348637, -0.0286146, -0.0224057, -0.0166353, -0.0116778, -0.00785986, -0.00543934,
2055 -0.00458888, -0.00538426, -0.00779888, -0.0117046, -0.0168789, -0.0230179, -0.029755, -0.0366829, -0.0433793, -0.0494324,
2056 -0.0544676, -0.0581708, -0.0603086, -0.0607437, -0.0594445, -0.0564874, -0.0520536, -0.0464185, -0.0399351, -0.0330126,
2057 -0.0260916, -0.0196165, -0.0140078, -0.00963553, -0.00679589, -0.00569172, -0.00641918, -0.00896067, -0.0131849, -0.0188543,
2058 -0.0256384, -0.0331341, -0.0408904, -0.0484353, -0.0553062, -0.0610783, -0.0653917, -0.0679743, -0.0686603, -0.0674007,
2059 -0.0642684, -0.0594539, -0.0532552, -0.0460599, -0.0383218, -0.0305332, -0.0231947, -0.0167834, -0.0117234, -0.00835846,
2060 -0.0069298, -0.00755992, -0.0102438, -0.0148479, -0.0211172, -0.0286902, -0.0371199, -0.0459016, -0.0545032, -0.0623987,
2061 -0.0691008, -0.0741922, -0.077352, -0.0783771, -0.077196, -0.0738746, -0.0686137, -0.0617378, -0.0536762, -0.0449369,
2062 -0.0360764, -0.0276647, -0.0202502, -0.0143252, -0.0102942, -0.008448, -0.00894392, -0.0117948, -0.0168663, -0.0238839,
2063 -0.0324485, -0.0420595, -0.0521457, -0.0620996, -0.0713153, -0.0792263, -0.0853413, -0.089276, -0.0907782, -0.0897453,
2064 -0.0862323, -0.0804503, -0.0727554, -0.0636277, -0.0536432, -0.0434389, -0.0336736, -0.0249868, -0.0179592, -0.0130748,
2065 -0.0106902, -0.0110102, -0.0140735, -0.0197479, -0.0277365, -0.0375944, -0.0487542, -0.0605603, -0.0723086, -0.0832897,
2066 -0.0928335, -0.100351, -0.105374, -0.107582, -0.10683, -0.103155, -0.0967769, -0.0880903, -0.077639, -0.0660858,
2067 -0.0541717, -0.0426718, -0.0323462, -0.0238935, -0.0179052, -0.0148282, -0.0149342, -0.0183005, -0.0248024, -0.034118,
2068 -0.0457455, -0.0590319, -0.0732118, -0.087453, -0.100908, -0.112767, -0.122308, -0.128944, -0.132262, -0.132051,
2069 -0.128318, -0.121292, -0.111413, -0.0993066, -0.0857491, -0.0716216, -0.0578565, -0.0453796, -0.035053, -0.0276199,
2070 -0.0236565, -0.023533, -0.0273873, -0.0351128, -0.0463603, -0.0605565, -0.0769351, -0.0945821, -0.11249, -0.129618,
2071 -0.14496, -0.157603, -0.166788, -0.171963, -0.172815, -0.169302, -0.161653, -0.150362, -0.136165, -0.119994,
2072 -0.102924, -0.0861121, -0.0707229, -0.0578581, -0.0484864, -0.0433804, -0.0430651, -0.0477797, -0.0574575, -0.0717239,
2073 -0.0899131, -0.111104, -0.134173, -0.157857, -0.180833, -0.201798, -0.219549, -0.233061, -0.241558, -0.244559,
2074 -0.241926, -0.23387, -0.220952, -0.204054, -0.184333, -0.163154, -0.142012, -0.122441, -0.105921, -0.0937841,
2075 -0.0871277, -0.0867421, -0.0930531, -0.106087, -0.125457, -0.150382, -0.179718, -0.212025, -0.245647, -0.278808,
2076 -0.309722, -0.3367, -0.358262, -0.373232, -0.380818, -0.380678, -0.37295, -0.35826, -0.337697, -0.312761,
2077 -0.285283, -0.257322, -0.231047, -0.208605, -0.191988, -0.182909, -0.182683, -0.192137, -0.211543, -0.240579,
2078 -0.278333, -0.323338, -0.373643, -0.426914, -0.480568, -0.531922, -0.578351, -0.617456, -0.647223, -0.666155,
2079 -0.673394, -0.668795, -0.652969, -0.627279, -0.593792, -0.555182, -0.514604, -0.475519, -0.441506, -0.416053,
2080 -0.402343, -0.40305, -0.420159, -0.454811, -0.507198, -0.5765, -0.660889, -0.757582, -0.862956, -0.972725,
2081 -1.08215, -1.18631, -1.28036, -1.35986, -1.42104, -1.4611, -1.47841, -1.47276, -1.44542, -1.39926,
2082 -1.33871, -1.26967, -1.19934, -1.13599, -1.08863, -1.06669, -1.0796, -1.13642, -1.24539, -1.41362,
2083 -1.64666, -1.94827, -2.3202, -2.76204, -3.27118, -3.84283, -4.47024, -5.14482, -5.85654, -6.59422,
2084 -7.346, -8.09977, -8.8436, -9.56622, -10.2574, -10.9083, -11.5118, -12.0628, -12.558, -12.9966,
2085 -13.3795, -13.7095, -13.9914, -14.2311, -14.4357, -14.6129, -14.7708, -14.9172, -15.0598, -15.205,
2086 -15.3587, -15.5251, -15.7073, -15.9066, -16.1231, -16.3554, -16.6011, -16.8566, -17.1176, -17.3795,
2087 -17.6373, -17.8864, -18.1224, -18.3417, -18.5414, -18.7197, -18.8758, -19.0099, -19.1233, -19.2181,
2088 -19.2974, -19.3646, -19.4237, -19.4785, -19.5331, -19.5909, -19.655, -19.7277, -19.8107, -19.9046,
2089 -20.0092, -20.1237, -20.2462, -20.3744, -20.5055, -20.6364, -20.7639, -20.8848, -20.9965, -21.0966,
2090 -21.1834, -21.2557, -21.3133, -21.3565, -21.3864, -21.4048, -21.4139, -21.4164, -21.415, -21.4125,
2091 -21.4118, -21.4153, -21.4249, -21.4422, -21.468, -21.5025, -21.5453, -21.5953, -21.651, -21.7103,
2092 -21.7709, -21.8303, -21.8861, -21.9359, -21.9777, -22.0098, -22.0312, -22.0413, -22.0401, -22.0284,
2093 -22.0072, -21.9782, -21.9434, -21.9052, -21.8658, -21.8275, -21.7926, -21.7629, -21.7397, -21.7241,
2094 -21.7165, -21.7167, -21.724, -21.7372, -21.7548, -21.7748, -21.7951, -21.8136, -21.8281, -21.8367,
2095 -21.8379, -21.8305, -21.8137, -21.7873, -21.7518, -21.708, -21.6571, -21.6009, -21.5413, -21.4804,
2096 -21.4203, -21.363, -21.3105, -21.2642, -21.2251, -21.1938, -21.1706, -21.1548, -21.1457, -21.1419,
2097 -21.1418, -21.1434, -21.1446, -21.1435, -21.1381, -21.1267, -21.1079, -21.0807, -21.0446, -20.9997,
2098 -20.9462, -20.8853, -20.8182, -20.7465, -20.6722, -20.5971, -20.5233, -20.4525, -20.3865, -20.3264,
2099 -20.2731, -20.2272, -20.1886, -20.1568, -20.1311, -20.1101, -20.0924, -20.0762, -20.0599, -20.0417,
2100 -20.0201, -19.9938, -19.9617, -19.9233, -19.8783, -19.8271, -19.7703, -19.7088, -19.6441, -19.5776,
2101 -19.5111, -19.4462, -19.3845, -19.3273, -19.2759, -19.2309, -19.1927, -19.1613, -19.136, -19.1161,
2102 -19.1003, -19.087, -19.0746, -19.0614, -19.0455, -19.0255, -19, -18.968, -18.9288, -18.8823,
2103 -18.8287, -18.7688, -18.7036, -18.6345, -18.5632, -18.4915, -18.4213, -18.3541, -18.2918, -18.2354,
2104 -18.186, -18.144, -18.1094, -18.0817, -18.0601, -18.0432, -18.0294, -18.0169, -18.0038, -17.9879,
2105 -17.9675, -17.9411, -17.9072, -17.865, -17.8143, -17.755, -17.6878, -17.6137, -17.5343, -17.4513,
2106 -17.3669, -17.2829, -17.2016, -17.1248, -17.0541, -16.9908, -16.9356, -16.8886, -16.8497, -16.8178,
2107 -16.7916, -16.7693, -16.7489, -16.7279, -16.704, -16.6748, -16.6383, -16.5926, -16.5366, -16.4694,
2108 -16.3909, -16.3015, -16.2022, -16.0948, -15.9812, -15.8638, -15.7452, -15.6281, -15.515, -15.4081,
2109 -15.3092, -15.2196, -15.1401, -15.0704, -15.0098, -14.9568, -14.9094, -14.8649, -14.8202, -14.7721,
2110 -14.7173, -14.6527, -14.5755, -14.4834, -14.3746, -14.2485, -14.1049, -13.9447, -13.7696, -13.582,
2111 -13.3851, -13.1825, -12.9781, -12.7758, -12.5793, -12.3921, -12.2168, -12.0554, -11.9086, -11.7762,
2112 -11.6568, -11.5477, -11.4452, -11.3446, -11.2402, -11.126, -10.9953, -10.8418, -10.6592, -10.442,
2113 -10.1853, -9.88567, -9.54092, -9.15041, -8.71517, -8.23789, -7.72289, -7.17603, -6.60446, -6.01647,
2114 -5.42109, -4.82787, -4.24645, -3.68621, -3.15595, -2.66358, -2.21577, -1.81777, -1.47323, -1.18411,
2115 -0.950602, -0.771239, -0.642961, -0.56131, -0.520656, -0.514466, -0.535612, -0.576679, -0.630282, -0.689371,
2116 -0.747497, -0.799048, -0.839434, -0.865214, -0.874163, -0.865283, -0.83875, -0.795814, -0.738644, -0.670148,
2117 -0.593756, -0.513198, -0.432273, -0.354635, -0.283593, -0.221943, -0.17184, -0.134707, -0.111193, -0.101173,
2118 -0.103795, -0.117558, -0.140431, -0.169992, -0.203579, -0.238455, -0.271968, -0.301693, -0.325563, -0.34197,
2119 -0.34984, -0.348664, -0.338509, -0.319988, -0.294198, -0.262639, -0.227107, -0.18958, -0.152089, -0.116601,
2120 -0.0849043, -0.0585064, -0.0385575, -0.0257941, -0.0205093, -0.0225507, -0.031344, -0.0459406, -0.0650864, -0.087305,
2121 -0.110993, -0.134518, -0.156316, -0.174983, -0.189353, -0.19856, -0.202077, -0.199742, -0.191754, -0.17865,
2122 -0.161263, -0.140664, -0.118092, -0.0948737, -0.0723401, -0.0517494, -0.0342117, -0.0206266, -0.0116346, -0.00758542,
2123 -0.00852463, -0.0141983, -0.0240761, -0.0373896, -0.053185, -0.0703843, -0.0878536, -0.104472, -0.119201, -0.13114,
2124 -0.139583, -0.144052, -0.144324, -0.140436, -0.132677, -0.121565, -0.107814, -0.0922781, -0.075906, -0.0596757,
2125 -0.0445367, -0.0313527, -0.0208515, -0.0135836, -0.00989249, -0.00989815, -0.0134945, -0.0203604, -0.0299833, -0.0416941,
2126 -0.0547112, -0.0681894, -0.0812731, -0.0931478, -0.103089, -0.110504, -0.114966, -0.116236, -0.114275, -0.109241,
2127 -0.101476, -0.0914837, -0.0798936, -0.0674217, -0.054825, -0.042855, -0.0322117, -0.0235022, -0.0172056, -0.0136455,
2128 -0.0129735, -0.0151626, -0.0200122, -0.0271629, -0.0361212, -0.0462918, -0.0570154, -0.0676102, -0.0774133, -0.0858214,
2129 -0.0923259, -0.096543, -0.0982346, -0.0973208, -0.0938823, -0.0881527, -0.0805021, -0.0714122, -0.0614456, -0.0512097,
2130 -0.0413194, -0.0323593, -0.0248484, -0.0192095, -0.0157447, -0.0146188, -0.0158514, -0.0193175, -0.0247573, -0.0317941,
2131 -0.0399591, -0.048721, -0.0575196, -0.0658003, -0.073048, -0.0788181, -0.0827625, -0.0846499, -0.0843778, -0.0819777,
2132 -0.077611, -0.0715572, -0.0641956, -0.0559808, -0.0474136, -0.0390097, -0.0312686, -0.0246423, -0.0195091, -0.0161507,
2133 -0.0147367, -0.0153153, -0.0178121, -0.0220358, -0.0276921, -0.0344022, -0.0417275, -0.049197, -0.0563367, -0.0626988,
2134 -0.0678883, -0.0715871, -0.0735725, -0.0737295, -0.0720566, -0.0686649, -0.0637696, -0.0576758, -0.0507586, -0.0434396,
2135 -0.03616, -0.0293538, -0.0234211, -0.0187039, -0.0154656, -0.0138762, -0.0140023, -0.0158046, -0.0191413, -0.0237779,
2136 -0.0294027, -0.0356467, -0.0421067, -0.0483709, -0.0540438, -0.0587705, -0.062258, -0.0642926, -0.0647525, -0.0636136,
2137 -0.0609507, -0.0569314, -0.0518052, -0.045887, -0.0395375, -0.0331401, -0.0270776, -0.0217081, -0.0173439, -0.0142322,
2138 -0.0125409, -0.0123488, -0.0136418, -0.0163146, -0.0201781, -0.0249717, -0.0303802, -0.0360537, -0.0416292, -0.0467532,
2139 -0.0511037, -0.054409, -0.0564645, -0.0571446, -0.0564096, -0.0543071, -0.0509689, -0.0466015, -0.0414731, -0.0358963,
2140 -0.0302087, -0.0247517, -0.0198494, -0.0157891, -0.0128035, -0.0110573, -0.0106373, -0.011548, -0.0137118, -0.0169747,
2141 -0.0211162, -0.0258637, -0.0309095, -0.0359305, -0.0406075, -0.0446455, -0.0477914, -0.0498492, -0.0506919, -0.0502681,
2142 -0.0486055, -0.0458076, -0.0420474, -0.0375557, -0.0326067, -0.0275007, -0.0225453, -0.0180362, -0.0142395, -0.0113752,
2143 -0.00960402, -0.00901811, -0.00963575, -0.0114009, -0.014187, -0.0178059, -0.0220192, -0.0265538, -0.0311189, -0.0354242,
2144 -0.0391975, -0.0422017, -0.044249, -0.0452121, -0.0450316, -0.0437194, -0.0413569, -0.0380902, -0.0341199, -0.0296889,
2145 -0.0250665, -0.0205323, -0.0163581, -0.0127917, -0.0100418, -0.00826492, -0.00755652, -0.00794541, -0.00939243, -0.0117934,
2146 -0.0149858, -0.0187596, -0.0228703, -0.0270542, -0.0310451, -0.0345904, -0.037467, -0.0394947, -0.0405469, -0.0405584,
2147 -0.0395287, -0.0375216, -0.0346612, -0.0311235, -0.0271251, -0.0229098, -0.0187333, -0.014847, -0.0114832, -0.0088399,
2148 -0.00706962, -0.00626981, -0.0064773, -0.00766644, -0.00975103, -0.0125901, -0.0159969, -0.0197507, -0.023611, -0.027332,
2151 Float_t output_100_diff[1001] =
2152 { 0.0194742, 0.0226635, 0.0255329, 0.0279047, 0.0296322, 0.0306094, 0.0307774, 0.0301287, 0.0287076, 0.0266079,
2153 0.0239667, 0.0209559, 0.0177718, 0.0146223, 0.0117138, 0.00923819, 0.00736053, 0.00620849, 0.00586409, 0.00635833,
2154 0.00766905, 0.00972214, 0.0123961, 0.0155295, 0.0189313, 0.0223927, 0.0257005, 0.0286504, 0.0310605, 0.0327824,
2155 0.0337111, 0.0337919, 0.0330238, 0.0314599, 0.0292042, 0.0264052, 0.0232465, 0.0199355, 0.01669, 0.0137244,
2156 0.0112363, 0.00939359, 0.00832371, 0.00810544, 0.00876369, 0.0102677, 0.0125327, 0.0154252, 0.0187715, 0.0223682,
2157 0.0259953, 0.0294301, 0.0324614, 0.0349031, 0.0366059, 0.0374674, 0.0374382, 0.0365259, 0.0347943, 0.03236,
2158 0.0293848, 0.0260659, 0.0226232, 0.0192853, 0.0162755, 0.0137967, 0.0120191, 0.0110687, 0.0110195, 0.0118882,
2159 0.0136331, 0.0161564, 0.0193104, 0.0229065, 0.0267274, 0.0305404, 0.0341126, 0.0372256, 0.0396897, 0.0413558,
2160 0.0421254, 0.0419571, 0.0408694, 0.0389397, 0.0362999, 0.0331279, 0.0296366, 0.0260602, 0.0226399, 0.0196076,
2161 0.0171719, 0.0155041, 0.0147272, 0.0149078, 0.0160518, 0.0181034, 0.0209486, 0.0244222, 0.028318, 0.0324023,
2162 0.0364281, 0.0401512, 0.0433461, 0.04582, 0.0474259, 0.0480718, 0.047727, 0.0464247, 0.0442595, 0.0413825,
2163 0.037991, 0.0343169, 0.0306111, 0.0271284, 0.0241106, 0.021771, 0.0202804, 0.0197566, 0.0202555, 0.0217682,
2164 0.0242204, 0.0274775, 0.0313525, 0.0356181, 0.0400206, 0.0442966, 0.0481899, 0.0514674, 0.0539351, 0.0554501,
2165 0.0559306, 0.0553613, 0.0537946, 0.051348, 0.0481966, 0.0445623, 0.0406998, 0.0368802, 0.0333737, 0.0304322,
2166 0.028273, 0.0270645, 0.0269149, 0.0278649, 0.0298851, 0.0328767, 0.0366782, 0.0410747, 0.0458122, 0.0506133,
2167 0.0551952, 0.0592878, 0.0626516, 0.0650928, 0.0664771, 0.0667377, 0.0658811, 0.0639866, 0.0612017, 0.0577335,
2168 0.0538351, 0.0497898, 0.0458932, 0.0424335, 0.0396734, 0.0378325, 0.0370727, 0.0374875, 0.0390953, 0.0418373,
2169 0.0455815, 0.05013, 0.0552316, 0.0605978, 0.0659214, 0.070896, 0.0752362, 0.0786965, 0.0810874, 0.0822884,
2170 0.0822565, 0.0810294, 0.0787243, 0.0755305, 0.0716978, 0.0675211, 0.063321, 0.0594237, 0.0561392, 0.0537413,
2171 0.0524494, 0.0524131, 0.0537023, 0.0563016, 0.0601105, 0.0649493, 0.0705697, 0.0766708, 0.082918, 0.0889646,
2172 0.094475, 0.099146, 0.102728, 0.10504, 0.105987, 0.105562, 0.103849, 0.10102, 0.0973253, 0.0930751,
2173 0.0886221, 0.084338, 0.0805894, 0.077714, 0.0759974, 0.0756539, 0.0768116, 0.0795024, 0.083659, 0.0891179,
2174 0.0956288, 0.10287, 0.110469, 0.118025, 0.125139, 0.131433, 0.136584, 0.14034, 0.142541, 0.143129,
2175 0.142156, 0.13978, 0.136258, 0.13193, 0.127198, 0.122499, 0.118282, 0.114969, 0.112937, 0.112485,
2176 0.113814, 0.117015, 0.122057, 0.128789, 0.136949, 0.146174, 0.15603, 0.166031, 0.175678, 0.184488,
2177 0.192026, 0.197939, 0.201976, 0.204014, 0.204063, 0.20227, 0.198916, 0.194397, 0.189201, 0.183878,
2178 0.179007, 0.175154, 0.172839, 0.172495, 0.17444, 0.178854, 0.185761, 0.195022, 0.206345, 0.219293,
2179 0.233317, 0.247785, 0.262021, 0.275354, 0.287156, 0.296894, 0.304162, 0.308715, 0.310491, 0.309623,
2180 0.306432, 0.301419, 0.295232, 0.288633, 0.282447, 0.277514, 0.274632, 0.274502, 0.277678, 0.284526,
2181 0.295193, 0.309586, 0.327372, 0.347989, 0.370676, 0.394512, 0.418475, 0.441504, 0.462566, 0.48073,
2182 0.495234, 0.505539, 0.511383, 0.512807, 0.510173, 0.504155, 0.495709, 0.486034, 0.476501, 0.468576,
2183 0.463731, 0.463352, 0.468643, 0.480538, 0.499632, 0.526115, 0.559745, 0.599836, 0.64527, 0.694553,
2184 0.745879, 0.797231, 0.846496, 0.891594, 0.930618, 0.961972, 0.984496, 0.997585, 1.00127, 0.996305,
2185 0.984152, 0.967007, 0.947739, 0.929802, 0.917114, 0.913906, 0.924547, 0.95335, 1.00438, 1.08124,
2186 1.18692, 1.32363, 1.49263, 1.69421, 1.92758, 2.19094, 2.48144, 2.79537, 3.12823, 3.47493,
2187 3.83, 4.1878, 4.54274, 4.88954, 5.22342, 5.54028, 5.83687, 6.1109, 6.36112, 6.58732,
2188 6.7903, 6.97181, 7.13439, 7.28127, 7.41612, 7.54291, 7.66564, 7.78817, 7.91402, 8.04619,
2189 8.18706, 8.33825, 8.50058, 8.67409, 8.85806, 9.05105, 9.2511, 9.45577, 9.6624, 9.86821,
2190 10.0705, 10.2668, 10.4551, 10.6338, 10.8019, 10.9591, 11.1058, 11.2427, 11.3713, 11.4936,
2191 11.6118, 11.7283, 11.8454, 11.9656, 12.0908, 12.2228, 12.3626, 12.5111, 12.6682, 12.8334,
2192 13.0059, 13.1839, 13.3659, 13.5496, 13.7329, 13.9136, 14.0898, 14.2598, 14.4222, 14.5763,
2193 14.7217, 14.8587, 14.988, 15.1109, 15.2288, 15.3437, 15.4577, 15.5729, 15.6911, 15.8143,
2194 15.9437, 16.0805, 16.2252, 16.3777, 16.5376, 16.7039, 16.8751, 17.0496, 17.2254, 17.4004,
2195 17.5726, 17.7402, 17.9016, 18.0556, 18.2015, 18.3389, 18.4683, 18.5902, 18.7061, 18.8175,
2196 18.9262, 19.0343, 19.1438, 19.2567, 19.3747, 19.4991, 19.631, 19.7707, 19.9181, 20.0728,
2197 20.2335, 20.3988, 20.5669, 20.7356, 20.9028, 21.0665, 21.2247, 21.3758, 21.5185, 21.6523,
2198 21.7768, 21.8924, 22, 22.101, 22.1972, 22.2905, 22.3832, 22.4775, 22.5755, 22.679,
2199 22.7893, 22.9076, 23.0341, 23.1686, 23.3105, 23.4584, 23.6105, 23.7647, 23.9186, 24.0698,
2200 24.2158, 24.3545, 24.484, 24.603, 24.7108, 24.8071, 24.8925, 24.968, 25.0352, 25.0963,
2201 25.1537, 25.2099, 25.2676, 25.3291, 25.3965, 25.4715, 25.5552, 25.6477, 25.7489, 25.8576,
2202 25.9721, 26.09, 26.2088, 26.3253, 26.4363, 26.5389, 26.6303, 26.7082, 26.7708, 26.8171,
2203 26.847, 26.861, 26.8606, 26.8479, 26.8257, 26.7969, 26.7651, 26.7334, 26.7051, 26.6829,
2204 26.6687, 26.6637, 26.6684, 26.6819, 26.7026, 26.7279, 26.7543, 26.7777, 26.7935, 26.797,
2205 26.7836, 26.749, 26.6896, 26.6027, 26.4865, 26.3405, 26.1654, 25.9632, 25.737, 25.4908,
2206 25.2293, 24.9577, 24.6813, 24.4048, 24.1326, 23.8679, 23.6124, 23.3666, 23.1287, 22.8955,
2207 22.6617, 22.4203, 22.163, 21.88, 21.5611, 21.1958, 20.7738, 20.2859, 19.724, 19.0824,
2208 18.3574, 17.5485, 16.658, 15.6916, 14.6582, 13.5698, 12.4411, 11.2892, 10.1332, 8.99309,
2209 7.88962, 6.84309, 5.87276, 4.99609, 4.22806, 3.58054, 3.06181, 2.67627, 2.42417, 2.30165,
2210 2.30085, 2.41025, 2.61505, 2.89783, 3.23917, 3.61844, 4.01456, 4.40684, 4.77568, 5.10333,
2211 5.3745, 5.57682, 5.70123, 5.74218, 5.69772, 5.56935, 5.36185, 5.08291, 4.74262, 4.35297,
2212 3.92722, 3.47929, 3.02313, 2.57212, 2.13855, 1.73318, 1.36484, 1.04025, 0.763831, 0.537724,
2213 0.361867, 0.23419, 0.150889, 0.106761, 0.0955836, 0.110522, 0.144525, 0.190704, 0.242677, 0.294848,
2214 0.342624, 0.382561, 0.412424, 0.431185, 0.438934, 0.436745, 0.426488, 0.4106, 0.391852, 0.373098,
2215 0.357048, 0.346059, 0.341968, 0.345968, 0.358534, 0.379407, 0.407624, 0.441596, 0.479235, 0.518097,
2216 0.555558, 0.588994, 0.615952, 0.63432, 0.642455, 0.639294, 0.624419, 0.598079, 0.561177, 0.515207,
2217 0.462166, 0.404426, 0.344597, 0.285364, 0.229338, 0.178898, 0.136065, 0.10239, 0.0788746, 0.0659282,
2218 0.0633603, 0.070408, 0.085797, 0.107832, 0.134511, 0.163648, 0.193018, 0.220482, 0.244117, 0.262323,
2219 0.273909, 0.278153, 0.274828, 0.264201, 0.247001, 0.224357, 0.197714, 0.168738, 0.1392, 0.110863,
2220 0.0853726, 0.0641537, 0.0483293, 0.038655, 0.0354816, 0.0387418, 0.0479645, 0.0623139, 0.0806502, 0.101608,
2221 0.123687, 0.145351, 0.165123, 0.181677, 0.193922, 0.20106, 0.202636, 0.198557, 0.189088, 0.174836,
2222 0.156695, 0.135794, 0.113415, 0.0909116, 0.0696226, 0.0507877, 0.0354711, 0.0244977, 0.0184055, 0.0174166,
2223 0.0214286, 0.0300267, 0.0425149, 0.0579648, 0.0752775, 0.0932556, 0.11068, 0.126388, 0.139344, 0.148702,
2224 0.15386, 0.154489, 0.150551, 0.142299, 0.130255, 0.115172, 0.0979861, 0.0797511, 0.061574, 0.0445429,
2225 0.029659, 0.0177749, 0.00954321, 0.00537846, 0.00543411, 0.00959667, 0.0174965, 0.0285347, 0.0419235, 0.0567387,
2226 0.0719795, 0.0866324, 0.0997348, 0.110435, 0.118045, 0.12208, 0.122288, 0.118661, 0.111433, 0.101062,
2227 0.0882014, 0.0736502, 0.0583066, 0.0431076, 0.0289707, 0.0167357, 0.00711347, 0.000642535, -0.00234263, -0.00173019,
2228 0.00236465, 0.00960958, 0.0194785, 0.031287, 0.044237, 0.0574682, 0.070113, 0.0813509, 0.0904597, 0.0968595,
2229 0.100148, 0.100122, 0.0967917, 0.0903737, 0.0812775, 0.0700766, 0.0574713, 0.0442432, 0.0312048, 0.0191487,
2230 0.00879776, 0.000760401, -0.0045064, -0.0067242, -0.00580874, -0.00187325, 0.00478068, 0.013684, 0.0242313, 0.0357204,
2231 0.047397, 0.058503, 0.0683241, 0.076234, 0.0817328, 0.0844778, 0.0843024, 0.081226, 0.0754499, 0.067344,
2232 0.0574208, 0.0463027, 0.0346815, 0.0232744, 0.0127786, 0.00382772, -0.00304721, -0.00745043, -0.00914488, -0.00806568,
2233 -0.00432264, 0.00180825, 0.00990365, 0.0194208, 0.0297326, 0.040168, 0.0500541, 0.0587587, 0.0657299, 0.0705293,
2234 0.0728594, 0.0725803, 0.0697176, 0.0644593, 0.0571422, 0.0482301, 0.038283, 0.0279213, 0.0177867, 0.00850097,
2235 0.000627825, -0.00536248, -0.00912172, -0.0104446, -0.0092804, -0.00573501,-6.28809e-05, 0.00735011, 0.016012, 0.0253563,
2236 0.0347787, 0.0436746, 0.0514775, 0.0576943, 0.0619356, 0.0639392, 0.0635858, 0.060905, 0.0560728, 0.0493988,
2237 0.0413064, 0.0323047, 0.0229562, 0.0138409, 0.00552003, -0.00149942, -0.00679589, -0.0100591, -0.0111089, -0.00990526,
2238 -0.00654997, -0.00127922, 0.00555238, 0.0134948, 0.0220318, 0.0306133, 0.0386909, 0.0457517, 0.0513506, 0.0551376,
2239 0.0568791, 0.0564717, 0.0539479, 0.0494731, 0.0433339, 0.0359198, 0.0276975, 0.0191815, 0.010901, 0.00336685,
2240 -0.00296033, -0.00769898, -0.0105699, -0.011413, -0.0101966, -0.00701865, -0.00209951, 0.00423263, 0.0115631, 0.0194173,
2241 0.0272912, 0.0346826, 0.0411234, 0.0462082, 0.0496194, 0.051147, 0.0507005, 0.0483141, 0.0441439, 0.0384569,
2242 0.0316137, 0.0240453, 0.0162255, 0.00864104, 0.00176056, -0.0039943, -0.00827528, -0.0108289, -0.0115116, -0.0102984,
2243 -0.00728351, -0.00267348, 0.00322611, 0.0100307, 0.0173012, 0.0245722, 0.031381, 0.037297, 0.0419483, 0.0450446,
2244 0.0463954, 0.0459208, 0.0436561, 0.0397491, 0.0344497, 0.0280939, 0.021082, 0.0138532, 0.00685788, 0.000528868,
2245 -0.0047452, -0.00864425, -0.0109365, -0.0114922, -0.0102918, -0.00742571, -0.00308868, 0.00243339, 0.00878187, 0.0155483,
2246 0.0223003, 0.028609, 0.0340759, 0.0383574, 0.0411868, 0.0423897, 0.0418952, 0.0397393, 0.0360623, 0.0310991,
2247 0.0251644, 0.0186319, 0.011911, 0.00542069,-0.000436952, -0.00530164, -0.00887746, -0.010951, -0.0114045, -0.0102224,
2248 -0.00749195, -0.00339729, 0.00179289, 0.00774252, 0.0140697, 0.0203707, 0.026246, 0.0313246, 0.0352877, 0.0378883,
2249 0.0389664, 0.0384582, 0.0364002, 0.032926, 0.0282574, 0.0226901, 0.016575, 0.0102952, 0.00424248, -0.00120778,
2250 -0.00571991, -0.00901886, -0.0109072, -0.011277, -0.0101164, -0.00750965, -0.00363127, 0.00126513, 0.00686331, 0.0128046,
2251 0.0187105, 0.0242068, 0.0289469, 0.0326332, 0.0350361, 0.0360076, 0.0354904, 0.0335211, 0.0302271, 0.0258188,
2254 Float_t output_diff_100[1001] =
2255 { -0.0366669, -0.0258292, -0.016614, -0.0096185, -0.00530268, -0.00395997, -0.00569845, -0.0104336, -0.0178929, -0.0276325,
2256 -0.0390646, -0.0514936, -0.0641597, -0.0762858, -0.0871267, -0.0960156, -0.102406, -0.105907, -0.106306, -0.103586,
2257 -0.0979244, -0.089684, -0.0793898, -0.0676974, -0.055352, -0.0431418, -0.0318492, -0.0222018, -0.0148261, -0.0102086,
2258 -0.00866466, -0.0103183, -0.0150937, -0.0227198, -0.032746, -0.0445704, -0.0574769, -0.0706795, -0.0833714, -0.0947753,
2259 -0.104192, -0.111046, -0.114919, -0.115579, -0.112993, -0.107334, -0.0989644, -0.0884205, -0.0763749, -0.0635967,
2260 -0.0509033, -0.0391096, -0.0289766, -0.0211641, -0.0161889, -0.014392, -0.0159168, -0.0206993, -0.0284718, -0.0387788,
2261 -0.0510049, -0.0644129, -0.0781895, -0.0914956, -0.103519, -0.113526, -0.120906, -0.125211, -0.126184, -0.123776,
2262 -0.118149, -0.109668, -0.0988764, -0.0864653, -0.0732289, -0.0600161, -0.0476773, -0.0370116, -0.028716, -0.0233424,
2263 -0.0212618, -0.0226412, -0.0274319, -0.0353724, -0.0460035, -0.0586965, -0.0726922, -0.0871474, -0.101188, -0.113961,
2264 -0.124694, -0.132735, -0.1376, -0.139004, -0.136872, -0.131356, -0.122816, -0.111802, -0.0990245, -0.0853026,
2265 -0.0715197, -0.058566, -0.0472832, -0.0384122, -0.0325462, -0.0300932, -0.0312502, -0.0359897, -0.0440611, -0.0550052,
2266 -0.0681823, -0.0828112, -0.0980178, -0.112889, -0.12653, -0.138118, -0.146958, -0.152523, -0.154488, -0.152751,
2267 -0.147443, -0.138917, -0.127727, -0.114598, -0.100377, -0.0859826, -0.0723494, -0.0603674, -0.0508275, -0.0443721,
2268 -0.0414542, -0.0423087, -0.0469367, -0.0551046, -0.0663581, -0.0800491, -0.0953759, -0.111432, -0.127262, -0.141922,
2269 -0.154537, -0.164356, -0.1708, -0.173498, -0.172309, -0.167337, -0.15892, -0.147612, -0.134151, -0.11941,
2270 -0.104348, -0.0899484, -0.0771561, -0.0668223, -0.0596491, -0.056146, -0.056597, -0.0610427, -0.0692763, -0.0808555,
2271 -0.0951292, -0.111277, -0.128359, -0.145374, -0.161321, -0.175262, -0.186381, -0.194031, -0.197782, -0.19744,
2272 -0.193068, -0.184979, -0.173717, -0.160028, -0.144809, -0.129061, -0.11382, -0.100097, -0.0888135, -0.0807473,
2273 -0.0764799, -0.0763617, -0.080489, -0.088697, -0.100569, -0.115459, -0.132536, -0.150826, -0.169281, -0.186836,
2274 -0.202477, -0.215307, -0.224596, -0.229833, -0.230753, -0.227357, -0.219917, -0.208951, -0.195204, -0.179594,
2275 -0.163159, -0.146997, -0.132194, -0.119762, -0.11057, -0.105296, -0.104381, -0.108003, -0.116063, -0.128189,
2276 -0.143759, -0.161938, -0.181728, -0.202028, -0.221701, -0.239642, -0.25485, -0.266486, -0.273927, -0.276805,
2277 -0.275031, -0.268801, -0.25859, -0.245116, -0.229308, -0.212241, -0.195079, -0.178998, -0.165121, -0.154444,
2278 -0.147776, -0.145691, -0.148488, -0.156176, -0.168466, -0.184793, -0.204342, -0.226101, -0.24892, -0.27158,
2279 -0.292864, -0.311634, -0.326898, -0.337874, -0.344034, -0.345143, -0.341272, -0.332796, -0.320373, -0.304912,
2280 -0.287511, -0.2694, -0.251866, -0.236173, -0.223492, -0.214827, -0.210952, -0.212366, -0.219262, -0.231509,
2281 -0.24866, -0.269977, -0.294471, -0.320959, -0.348132, -0.374635, -0.399143, -0.42044, -0.437497, -0.449528,
2282 -0.456044, -0.456881, -0.452216, -0.442559, -0.428727, -0.411799, -0.393058, -0.373914, -0.35583, -0.340234,
2283 -0.328436, -0.32155, -0.320433, -0.325625, -0.33732, -0.355348, -0.379181, -0.407957, -0.440529, -0.475521,
2284 -0.51141, -0.546603, -0.579532, -0.608742, -0.632975, -0.651245, -0.662897, -0.667653, -0.665634, -0.65736,
2285 -0.643734, -0.626, -0.605679, -0.584499, -0.564303, -0.546954, -0.534236, -0.527756, -0.528857, -0.538539,
2286 -0.557397, -0.585579, -0.622767, -0.668181, -0.720604, -0.778435, -0.839755, -0.90242, -0.964162, -1.0227,
2287 -1.07587, -1.12171, -1.15862, -1.1854, -1.20141, -1.20658, -1.20146, -1.18726, -1.16585, -1.13967,
2288 -1.11175, -1.08557, -1.06499, -1.0541, -1.05711, -1.0782, -1.12137, -1.19032, -1.28827, -1.41789,
2289 -1.58114, -1.77925, -2.01258, -2.28066, -2.58213, -2.9148, -3.2757, -3.66113, -4.06679, -4.48791,
2290 -4.91936, -5.35585, -5.79208, -6.22285, -6.64332, -7.04905, -7.43624, -7.80174, -8.14321, -8.45915,
2291 -8.74892, -9.01278, -9.25178, -9.46778, -9.66327, -9.84132, -10.0054, -10.1593, -10.3068, -10.4516,
2292 -10.5975, -10.7475, -10.9044, -11.0703, -11.2468, -11.4345, -11.6336, -11.8434, -12.0626, -12.2895,
2293 -12.5219, -12.7571, -12.9924, -13.2252, -13.4528, -13.6728, -13.8833, -14.0829, -14.2706, -14.4461,
2294 -14.6098, -14.7626, -14.906, -15.0418, -15.1725, -15.3005, -15.4284, -15.5588, -15.6941, -15.8364,
2295 -15.9872, -16.1477, -16.3185, -16.4995, -16.69, -16.8888, -17.0943, -17.3043, -17.5165, -17.7286,
2296 -17.9379, -18.1422, -18.3395, -18.5283, -18.7073, -18.8761, -19.0347, -19.1837, -19.3245, -19.4586,
2297 -19.5883, -19.7159, -19.844, -19.975, -20.1114, -20.2552, -20.4081, -20.5713, -20.7454, -20.9303,
2298 -21.1253, -21.3294, -21.5407, -21.7571, -21.9762, -22.1955, -22.4123, -22.6243, -22.8294, -23.0261,
2299 -23.2131, -23.39, -23.557, -23.7148, -23.8648, -24.009, -24.1496, -24.2892, -24.4306, -24.5766,
2300 -24.7295, -24.8917, -25.0647, -25.2498, -25.4473, -25.6572, -25.8785, -26.1098, -26.349, -26.5937,
2301 -26.8411, -27.0885, -27.333, -27.5721, -27.8037, -28.026, -28.2381, -28.4397, -28.6311, -28.8136,
2302 -28.9888, -29.1591, -29.3273, -29.4964, -29.6696, -29.8498, -30.04, -30.2425, -30.4589, -30.6904,
2303 -30.9372, -31.199, -31.4745, -31.7618, -32.0583, -32.3611, -32.667, -32.9727, -33.275, -33.5712,
2304 -33.859, -34.1365, -34.4031, -34.6586, -34.9038, -35.1404, -35.3708, -35.5981, -35.8258, -36.0576,
2305 -36.2974, -36.5488, -36.8152, -37.0993, -37.403, -37.7275, -38.073, -38.4386, -38.8227, -39.2227,
2306 -39.6353, -40.057, -40.4835, -40.911, -41.3355, -41.7537, -42.1629, -42.5612, -42.9478, -43.323,
2307 -43.6882, -44.0459, -44.3995, -44.7534, -45.1126, -45.4822, -45.8678, -46.2745, -46.707, -47.1693,
2308 -47.6643, -48.1938, -48.7584, -49.3572, -49.9881, -50.6478, -51.3319, -52.0352, -52.7521, -53.4768,
2309 -54.2036, -54.9274, -55.6442, -56.3509, -57.0462, -57.7303, -58.4053, -59.0753, -59.7461, -60.4252,
2310 -61.1216, -61.8452, -62.607, -63.4179, -64.2888, -65.2298, -66.2496, -67.3555, -68.5522, -69.8422,
2311 -71.2249, -72.6966, -74.2507, -75.8771, -77.5628, -79.2922, -81.0469, -82.8069, -84.5504, -86.255,
2312 -87.8979, -89.4568, -90.9107, -92.2403, -93.4286, -94.4617, -95.3293, -96.0246, -96.5453, -96.893,
2313 -97.074, -97.0983, -96.9802, -96.7375, -96.3909, -95.9638, -95.4813, -94.9695, -94.4549, -93.9632,
2314 -93.5189, -93.1441, -92.8581, -92.6764, -92.6103, -92.6664, -92.8465, -93.1468, -93.5588, -94.0687,
2315 -94.6583, -95.305, -95.983, -96.6637, -97.3166, -97.9107, -98.415, -98.7999, -99.0381, -99.1056,
2316 -98.9825, -98.6539, -98.11, -97.3472, -96.3676, -95.1793, -93.7961, -92.2368, -90.5247, -88.6867,
2317 -86.7522, -84.752, -82.717, -80.6773, -78.6605, -76.6913, -74.7897, -72.971, -71.2445, -69.6137,
2318 -68.0761, -66.6229, -65.2402, -63.909, -62.6065, -61.3069, -59.9831, -58.6072, -57.1528, -55.5954,
2319 -53.9144, -52.0937, -50.1227, -47.9971, -45.7193, -43.2983, -40.7499, -38.0957, -35.3632, -32.5841,
2320 -29.7935, -27.0286, -24.3272, -21.7264, -19.2611, -16.9628, -14.8584, -12.9693, -11.3105, -9.89043,
2321 -8.7105, -7.76534, -7.0432, -6.52648, -6.19271, -6.01548, -5.9657, -6.01276, -6.12587, -6.27521,
2322 -6.4331, -6.575, -6.68028, -6.73284, -6.72148, -6.64001, -6.48715, -6.26619, -5.98445, -5.6526,
2323 -5.28383, -4.89298, -4.49556, -4.1069, -3.74129, -3.41119, -3.12671, -2.89507, -2.72042, -2.60368,
2324 -2.54267, -2.53239, -2.56539, -2.63231, -2.72249, -2.82462, -2.9274, -3.02018, -3.09355, -3.13978,
2325 -3.15323, -3.13054, -3.07077, -2.9753, -2.84771, -2.69345, -2.51943, -2.33361, -2.14442, -1.96032,
2326 -1.78924, -1.63816, -1.5127, -1.41688, -1.35287, -1.32096, -1.31959, -1.34551, -1.39399, -1.45918,
2327 -1.53446, -1.61289, -1.68759, -1.75218, -1.80112, -1.83003, -1.83591, -1.81729, -1.77425, -1.7084,
2328 -1.62273, -1.52139, -1.4094, -1.29233, -1.17599, -1.06599, -0.96748, -0.884812, -0.821306, -0.779075,
2329 -0.758919, -0.760313, -0.781463, -0.819449, -0.87042, -0.929852, -0.992835, -1.05438, -1.10971, -1.15458,
2330 -1.18547, -1.19983, -1.19619, -1.17423, -1.13477, -1.07972, -1.01191, -0.934906, -0.852802, -0.769928,
2331 -0.690599, -0.618846, -0.558175, -0.511364, -0.480303, -0.465901, -0.468043, -0.485616, -0.516596, -0.558189,
2332 -0.607011, -0.659312, -0.711203, -0.758904, -0.79897, -0.828495, -0.845283, -0.847969, -0.836087, -0.810089,
2333 -0.771295, -0.721801, -0.664338, -0.602092, -0.538502, -0.477045, -0.421021, -0.37335, -0.336404, -0.311866,
2334 -0.300634, -0.302783, -0.317571, -0.343493, -0.378389, -0.419587, -0.464069, -0.508668, -0.550262, -0.585967,
2335 -0.61331, -0.630374, -0.635911, -0.62941, -0.611114, -0.581998, -0.543697, -0.498398, -0.448695, -0.397424,
2336 -0.347486, -0.301666, -0.262463, -0.231937, -0.211591, -0.202284, -0.204184, -0.216768, -0.238865, -0.268733,
2337 -0.304176, -0.342684, -0.381593, -0.418251, -0.450182, -0.475236, -0.49172, -0.498496, -0.495042, -0.481482,
2338 -0.458568, -0.427624, -0.390461, -0.349256, -0.306417, -0.264424, -0.225684, -0.192374, -0.166315, -0.148858,
2339 -0.14081, -0.142389, -0.153213, -0.172336, -0.19831, -0.229279, -0.2631, -0.297481, -0.33012, -0.35885,
2340 -0.381774, -0.397376, -0.404612, -0.402975, -0.39251, -0.373816, -0.347997, -0.316592, -0.28147, -0.244716,
2341 -0.208495, -0.174921, -0.145923, -0.123128, -0.107764, -0.100584, -0.101831, -0.11122, -0.127967, -0.150839,
2342 -0.178235, -0.20829, -0.238996, -0.268322, -0.294346, -0.315369, -0.330022, -0.337347, -0.336849, -0.32853,
2343 -0.312876, -0.290829, -0.263718, -0.233179, -0.201045, -0.169236, -0.139632, -0.113963, -0.0936982, -0.0799579,
2344 -0.0734487, -0.0744216, -0.0826608, -0.0975011, -0.117873, -0.142374, -0.169357, -0.197035, -0.223596, -0.247312,
2345 -0.266649, -0.28036, -0.287559, -0.287774, -0.280974, -0.267564, -0.24836, -0.224532, -0.197529, -0.168988,
2346 -0.140626, -0.11414, -0.0910929, -0.0728244, -0.0603648, -0.0543754, -0.0551096, -0.0624006, -0.0756752, -0.0939923,
2347 -0.116105, -0.14054, -0.165692, -0.189924, -0.21167, -0.229531, -0.242363, -0.249344, -0.250023, -0.244347,
2348 -0.23266, -0.21568, -0.194446, -0.17026, -0.144595, -0.119007, -0.0950368, -0.0741122, -0.0574613, -0.0460371,
2349 -0.0404601, -0.040982, -0.0474728, -0.0594317, -0.0760208, -0.0961202, -0.1184, -0.141403, -0.163639, -0.183677,
2350 -0.200234, -0.212253, -0.218969, -0.219952, -0.215136, -0.204816, -0.189629, -0.170512, -0.148642, -0.125358,
2351 -0.10208, -0.0802181, -0.0610843, -0.0458119, -0.0352857, -0.0300887, -0.0304677, -0.0363199, -0.0472022, -0.0623604,
2352 -0.0807791, -0.101246, -0.122431, -0.142965, -0.161533, -0.176949, -0.188232, -0.194667, -0.195847, -0.191694,
2353 -0.182467, -0.168737, -0.151355, -0.131395, -0.110084, -0.0887271, -0.0686232, -0.0509856, -0.036866, -0.0270899,
2354 -0.0222065, -0.0224563, -0.027758, -0.0377154, -0.0516442, -0.0686165, -0.0875206, -0.107131, -0.126186, -0.143466,
2357 Float_t output_113_113[1001] =
2358 { -0.0252237, -0.0236161, -0.021104, -0.0178422, -0.0140333, -0.00991613, -0.00575021, -0.00179993, 0.00168219, 0.00447149,
2359 0.00638547, 0.0072956, 0.00713573, 0.00590671, 0.00367676, 0.000577596, -0.00320339, -0.00743462, -0.0118546, -0.016188,
2360 -0.0201635, -0.0235301, -0.0260735, -0.0276301, -0.0280971, -0.0274397, -0.0256936, -0.0229637, -0.0194171, -0.015274,
2361 -0.0107934, -0.00625743, -0.00195376, 0.00184263, 0.00488678, 0.00697941, 0.00797976, 0.00781493, 0.00648491, 0.00406311,
2362 0.000692253, -0.00342408, -0.00803389, -0.0128522, -0.017579, -0.0219175, -0.0255936, -0.0283726, -0.0300749, -0.030587,
2363 -0.02987, -0.0279617, -0.0249756, -0.0210938, -0.016556, -0.0116453, -0.00667038, -0.00194609, 0.00222596, 0.00557674,
2364 0.00788691, 0.00900103, 0.00883786, 0.00739601, 0.00475468, 0.00106925, -0.00343786, -0.0084907, -0.0137768, -0.0189664,
2365 -0.0237333, -0.0277753, -0.0308334, -0.0327086, -0.0332748, -0.0324871, -0.0303859, -0.0270943, -0.0228117, -0.0178012,
2366 -0.0123742, -0.00687059, -0.00163808, 0.00298998, 0.00671557, 0.00929508, 0.0105551, 0.0104039, 0.00883804, 0.00594309,
2367 0.00188918, -0.00307905, -0.00865725, -0.0145, -0.0202421, -0.0255217, -0.0300027, -0.0333965, -0.0354803, -0.036112,
2368 -0.0352395, -0.0329051, -0.0292435, -0.0244742, -0.0188883, -0.012831, -0.00668025,-0.000823154, 0.00436811, 0.00856002,
2369 0.011479, 0.0129293, 0.0128058, 0.0111017, 0.00790989, 0.00341769, -0.00210369, -0.00831562, -0.0148327, -0.0212466,
2370 -0.0271516, -0.0321698, -0.035976, -0.0383176, -0.039032, -0.0380577, -0.0354384, -0.0313222, -0.0259527, -0.0196549,
2371 -0.0128153, -0.00585844, 0.000779598, 0.00667855, 0.0114605, 0.0148142, 0.0165152, 0.016441, 0.0145797, 0.011032,
2372 0.00600597, -0.00019502, -0.00719025, -0.0145449, -0.0217968, -0.0284853, -0.0341802, -0.0385092, -0.0411819, -0.0420092,
2373 -0.0409158, -0.0379467, -0.0332654, -0.027145, -0.0199522, -0.0121251, -0.00414674, 0.00348496, 0.0102884, 0.0158289,
2374 0.0197463, 0.0217786, 0.0217797, 0.0197297, 0.0157381, 0.0100376, 0.00297157, -0.0050264, -0.0134589, -0.0217951,
2375 -0.0295037, -0.0360866, -0.0411103, -0.0442343, -0.0452334, -0.0440133, -0.0406185, -0.0352308, -0.0281597, -0.0198246,
2376 -0.0107299, -0.00143405, 0.00748456, 0.015464, 0.0219941, 0.0266494, 0.029117, 0.0292176, 0.0269186, 0.0223381,
2377 0.0157394, 0.00751661, -0.00182869, -0.0117173, -0.021528, -0.0306357, -0.0384511, -0.0444579, -0.0482465, -0.0495413,
2378 -0.0482198, -0.0443225, -0.0380529, -0.0297668, -0.0199523, -0.00920106, 0.00182733, 0.0124459, 0.021983, 0.0298246,
2379 0.0354541, 0.0384858, 0.0386915, 0.0360172, 0.0305888, 0.0227077, 0.0128338, 0.00155976, -0.0104245, -0.0223736,
2380 -0.0335324, -0.0431836, -0.0506931, -0.0555517, -0.0574087, -0.056097, -0.0516465, -0.0442855, -0.0344301, -0.0226614,
2381 -0.00969122, 0.00367928, 0.016609, 0.0282687, 0.0378932, 0.0448306, 0.0485852, 0.0488515, 0.0455369, 0.038771,
2382 0.0289025, 0.0164806, 0.0022248, -0.013017, -0.0283195, -0.042736, -0.0553582, -0.065373, -0.0721164, -0.0751174,
2383 -0.0741324, -0.0691642, -0.0604682, -0.0485408, -0.0340943, -0.0180165, -0.00131922, 0.0149222, 0.0296351, 0.041816,
2384 0.0505956, 0.0552956, 0.0554756, 0.0509666, 0.041887, 0.0286425, 0.0119078, -0.00740993, -0.0282215, -0.0493195,
2385 -0.0694512, -0.0873951, -0.102039, -0.112451, -0.117946, -0.118129, -0.112935, -0.102637, -0.0878392, -0.0694524,
2386 -0.0486418, -0.026764, -0.0052882, 0.0142918, 0.0305497, 0.0422169, 0.0482661, 0.0479823, 0.0410166, 0.0274194,
2387 0.00764954, -0.0174407, -0.0466459, -0.0784708, -0.111213, -0.143062, -0.172207, -0.196949, -0.215814, -0.227646,
2388 -0.231699, -0.227692, -0.215851, -0.196912, -0.172101, -0.143076, -0.111849, -0.0806766, -0.0519344, -0.0279793,
2389 -0.0110062, -0.00290916, -0.00515341, -0.0186681, -0.0437655, -0.0800935, -0.126625, -0.181686, -0.243025, -0.307916,
2390 -0.373302, -0.435956, -0.492669, -0.540441, -0.576682, -0.599388, -0.607307, -0.600073, -0.578292, -0.543596,
2391 -0.498636, -0.44703, -0.393258, -0.342507, -0.30048, -0.273167, -0.266594, -0.286556, -0.338361, -0.426573,
2392 -0.55479, -0.725455, -0.939712, -1.19731, -1.4966, -1.83451, -2.20668, -2.60761, -3.03082, -3.46912,
2393 -3.91487, -4.36025, -4.7976, -5.21966, -5.61984, -5.99248, -6.33304, -6.6382, -6.90597, -7.13574,
2394 -7.32817, -7.48519, -7.60977, -7.70581, -7.77785, -7.83088, -7.87008, -7.90054, -7.92705, -7.95392,
2395 -7.98474, -8.02231, -8.06851, -8.12429, -8.18968, -8.26385, -8.34521, -8.43156, -8.52022, -8.60825,
2396 -8.69263, -8.77041, -8.83893, -8.89591, -8.93963, -8.96897, -8.98345, -8.98327, -8.96926, -8.9428,
2397 -8.90574, -8.86029, -8.80884, -8.75387, -8.69778, -8.64274, -8.59062, -8.54285, -8.50038, -8.46363,
2398 -8.43249, -8.40636, -8.38417, -8.3645, -8.34566, -8.32577, -8.30295, -8.27538, -8.24143, -8.19973,
2399 -8.1493, -8.08954, -8.02032, -7.94193, -7.85506, -7.76079, -7.66048, -7.5557, -7.44815, -7.3395,
2400 -7.23137, -7.12519, -7.02213, -6.92304, -6.8284, -6.73832, -6.65253, -6.57041, -6.49103, -6.41323,
2401 -6.33568, -6.25698, -6.17573, -6.09064, -6.00059, -5.90472, -5.80244, -5.6935, -5.57798, -5.45631,
2402 -5.3292, -5.19764, -5.06279, -4.92597, -4.78851, -4.65173, -4.51684, -4.38487, -4.25663, -4.13264,
2403 -4.01311, -3.89796, -3.78682, -3.67903, -3.57372, -3.46986, -3.36631, -3.26192, -3.15557, -3.04627,
2404 -2.9332, -2.81575, -2.6936, -2.56669, -2.43525, -2.29978, -2.16102, -2.01992, -1.87754, -1.73506,
2405 -1.59363, -1.45438, -1.31831, -1.18624, -1.05878, -0.936292, -0.818871, -0.706345, -0.59829, -0.494057,
2406 -0.392817, -0.293606, -0.195392, -0.0971273, 0.00218312, 0.103428, 0.207336, 0.314439, 0.425038, 0.539191,
2407 0.656711, 0.777178, 0.899964, 1.02427, 1.14917, 1.27366, 1.39673, 1.51738, 1.63471, 1.74796,
2408 1.85651, 1.95994, 2.05803, 2.15079, 2.23839, 2.32119, 2.3997, 2.47451, 2.54629, 2.61572,
2409 2.68344, 2.75004, 2.81596, 2.88155, 2.94696, 3.0122, 3.07709, 3.14132, 3.20443, 3.26586,
2410 3.32495, 3.38105, 3.43347, 3.4816, 3.52487, 3.56284, 3.5952, 3.62176, 3.64251, 3.65758,
2411 3.6672, 3.67176, 3.67171, 3.66757, 3.65988, 3.64915, 3.63587, 3.62043, 3.60314, 3.58416,
2412 3.56352, 3.54112, 3.51671, 3.4899, 3.4602, 3.42702, 3.38971, 3.34757, 3.29992, 3.24609,
2413 3.18547, 3.11753, 3.04186, 2.95816, 2.86628, 2.76622, 2.65813, 2.54231, 2.41921, 2.28942,
2414 2.15367, 2.01278, 1.86772, 1.7195, 1.56924, 1.41813, 1.26739, 1.11829, 0.972119, 0.830169,
2415 0.693732, 0.564068, 0.442392, 0.32985, 0.227495, 0.136263, 0.0569418, -0.00985115, -0.0636982, -0.104405,
2416 -0.13202, -0.146851, -0.149477, -0.140743, -0.12176, -0.0938823, -0.058685, -0.0179251, 0.0265032, 0.072622,
2417 0.118429, 0.161963, 0.201365, 0.23495, 0.261261, 0.279126, 0.287705, 0.286519, 0.275471, 0.254854,
2418 0.225333, 0.187926, 0.143955, 0.094994, 0.0428023, -0.0107534, -0.0637787, -0.114435, -0.161022, -0.202049,
2419 -0.236308, -0.262921, -0.281383, -0.291582, -0.293804, -0.288716, -0.277333, -0.260968, -0.241169, -0.219639,
2420 -0.198158, -0.178487, -0.162283, -0.151018, -0.145895, -0.147792, -0.157206, -0.174225, -0.198513, -0.229317,
2421 -0.265494, -0.305552, -0.347714, -0.389991, -0.430267, -0.466387, -0.496251, -0.517908, -0.529637, -0.530024,
2422 -0.518021, -0.492996, -0.45476, -0.403576, -0.340148, -0.265599, -0.181418, -0.0894111, 0.00837995, 0.10975,
2423 0.212417, 0.31411, 0.412648, 0.50602, 0.592452, 0.670462, 0.738906, 0.797, 0.844334, 0.880869,
2424 0.906916, 0.923104, 0.930337, 0.92974, 0.922601, 0.91031, 0.894292, 0.875954, 0.85662, 0.837492,
2425 0.819605, 0.803797, 0.790696, 0.780706, 0.774014, 0.770599, 0.770259, 0.772634, 0.777243, 0.783518,
2426 0.790843, 0.798592, 0.806164, 0.813011, 0.818666, 0.822764, 0.82505, 0.825387, 0.823759, 0.820256,
2427 0.815069, 0.808471, 0.800794, 0.792412, 0.783713, 0.775083, 0.76688, 0.75942, 0.752959, 0.747687,
2428 0.743717, 0.741088, 0.739763, 0.739637, 0.740546, 0.74228, 0.744596, 0.747231, 0.749921, 0.752411,
2429 0.75447, 0.755903, 0.756559, 0.756334, 0.755179, 0.753097, 0.75014, 0.746406, 0.742031, 0.737176,
2430 0.732026, 0.726771, 0.721598, 0.716683, 0.712181, 0.708218, 0.704886, 0.702237, 0.700288, 0.699014,
2431 0.698354, 0.698215, 0.69848, 0.699011, 0.699659, 0.700273, 0.700706, 0.700823, 0.700509, 0.699673,
2432 0.698255, 0.696222, 0.693577, 0.690353, 0.68661, 0.682435, 0.677935, 0.673227, 0.66844, 0.663699,
2433 0.659126, 0.654826, 0.650889, 0.647379, 0.644337, 0.641773, 0.639671, 0.637987, 0.636652, 0.635578,
2434 0.634661, 0.633787, 0.63284, 0.631705, 0.630279, 0.62847, 0.62621, 0.623452, 0.620174, 0.616383,
2435 0.612112, 0.607417, 0.602376, 0.597083, 0.591645, 0.586172, 0.580776, 0.57556, 0.570615, 0.566012,
2436 0.561804, 0.558015, 0.554647, 0.551672, 0.54904, 0.546678, 0.544497, 0.542394, 0.540262, 0.537992,
2437 0.535481, 0.53264, 0.529396, 0.525698, 0.521519, 0.516859, 0.511744, 0.506225, 0.500374, 0.494283,
2438 0.488054, 0.481798, 0.475623, 0.469634, 0.463922, 0.458561, 0.453603, 0.449077, 0.444982, 0.441293,
2439 0.437959, 0.434909, 0.432051, 0.429285, 0.426501, 0.423592, 0.420458, 0.417009, 0.413178, 0.408916,
2440 0.404203, 0.399043, 0.393469, 0.387537, 0.381326, 0.374932, 0.36846, 0.362022, 0.355729, 0.349683,
2441 0.34397, 0.338659, 0.333795, 0.329396, 0.325453, 0.321932, 0.318772, 0.315892, 0.313195, 0.310576,
2442 0.307923, 0.305129, 0.302098, 0.298748, 0.295017, 0.29087, 0.286297, 0.281317, 0.275972, 0.270333,
2443 0.264485, 0.258533, 0.252586, 0.246756, 0.241152, 0.235867, 0.230979, 0.226545, 0.222594, 0.219129,
2444 0.216125, 0.213531, 0.211273, 0.209259, 0.207385, 0.205538, 0.20361, 0.201495, 0.199106, 0.196372,
2445 0.193245, 0.189706, 0.185762, 0.181447, 0.176822, 0.171967, 0.166981, 0.161971, 0.15705, 0.152327,
2446 0.147902, 0.143858, 0.140259, 0.137141, 0.134516, 0.132366, 0.130647, 0.129292, 0.128212, 0.127306,
2447 0.126461, 0.125567, 0.124517, 0.123216, 0.121588, 0.119579, 0.117161, 0.114334, 0.111123, 0.107584,
2448 0.10379, 0.0998372, 0.0958296, 0.0918792, 0.0880958, 0.0845814, 0.0814234, 0.0786893, 0.0764226, 0.0746399,
2449 0.0733301, 0.072455, 0.0719515, 0.0717357, 0.071708, 0.0717592, 0.0717772, 0.0716541, 0.0712927, 0.0706126,
2450 0.0695551, 0.0680867, 0.0662014, 0.0639205, 0.0612921, 0.0583873, 0.0552961, 0.0521216, 0.0489739, 0.0459631,
2451 0.043192, 0.0407506, 0.0387097, 0.0371167, 0.0359931, 0.0353325, 0.0351014, 0.0352412, 0.0356718, 0.0362965,
2452 0.0370081, 0.0376956, 0.0382507, 0.0385752, 0.0385863, 0.0382221, 0.0374456, 0.0362465, 0.034642, 0.032676,
2453 0.0304156, 0.0279478, 0.0253731, 0.0228, 0.020338, 0.0180903, 0.0161482, 0.0145846, 0.0134502, 0.0127697,
2454 0.0125406, 0.0127334, 0.0132933, 0.0141437, 0.0151909, 0.01633, 0.0174514, 0.0184475, 0.0192196, 0.0196839,
2455 0.0197766, 0.0194581, 0.0187149, 0.017561, 0.0160368, 0.0142063, 0.0121533, 0.009976, 0.00778115, 0.00567701,
2456 0.00376676, 0.00214209, 0.000877411, 2.52803e-05, -0.0003869,-0.000358562, 8.40117e-05, 0.000889335, 0.00198395, 0.00327702,
2457 0.00466598, 0.00604298, 0.00730153, 0.00834322, 0.00908375, 0.0094582, 0.00942507, 0.00896876, 0.00810059, 0.00685796,
2459 const size_t thousandFloats = 1001*
sizeof(Float_t);
2461 case 3 : memcpy(DataInTheGap,output_113_113 ,thousandFloats);
break;
2462 case 2 : memcpy(DataInTheGap,output_diff_100,thousandFloats);
break;
2463 case 1 : memcpy(DataInTheGap,output_100_diff,thousandFloats);
break;
2465 default : memcpy(DataInTheGap,output_nowall_oldTPC,thousandFloats);
2487 Double_t r, phi, z ;
2489 if (usingCartesian) Cart2Polar(x,r,phi);
2490 else { r = x[0]; phi = x[1]; }
2492 z = LimitZ( Sector, x ) ;
2494 if ( z < 0 ) phi += EASTCLOCKERROR/1000. ;
2495 if ( z > 0 ) phi += WESTCLOCKERROR/1000. ;
2498 if (usingCartesian) Polar2Cart(r,phi,Xprime);
2499 else { Xprime[0] = r; Xprime[1] = phi; }
2515 cout <<
"StMagUtilities::UndoMembrane This routine was made obosolete on 10/1/2009. Do not use it." << endl ;
2559 cout <<
"StMagUtilities::UndoEndcap This routine was made obosolete on 10/1/2009. Do not use it." << endl ;
2608 Float_t Er_integral, Ephi_integral ;
2609 Double_t r, phi, z ;
2611 const Int_t ORDER = 1 ;
2615 cout <<
"StMagUtilities::IFCShift Please wait for the tables to fill ... ~5 seconds" << endl ;
2616 Int_t Nterms = 100 ;
2617 Double_t Denominator[100];
2618 memset(Denominator,0,100*
sizeof(Double_t));
2619 for ( Int_t i = 0 ; i < EMap_nZ ; ++i )
2621 z = TMath::Abs( eZList[i] ) ;
2622 for ( Int_t j = 0 ; j < EMap_nR ; ++j )
2625 shiftEr[i][j] = 0.0 ;
2626 if (r < IFCRadius)
continue;
2627 if (r > OFCRadius)
continue;
2628 if (z > TPC_Z0)
continue;
2629 Double_t IntegralOverZ = 0.0 ;
2630 for ( Int_t n = 1 ; n < Nterms ; ++n )
2632 Double_t k = (2*n-1) * TMath::Pi() / TPC_Z0 ;
2633 Double_t Cn = -4.0 * IFCShift / ( k * TPC_Z0 ) ;
2634 Double_t Numerator =
2635 TMath::BesselK0( k*OFCRadius ) * TMath::BesselI1( k*r ) +
2636 TMath::BesselK1( k*r ) * TMath::BesselI0( k*OFCRadius ) ;
2637 if (Denominator[n] == 0) Denominator[n] =
2638 TMath::BesselK0( k*OFCRadius ) * TMath::BesselI0( k*IFCRadius ) -
2639 TMath::BesselK0( k*IFCRadius ) * TMath::BesselI0( k*OFCRadius ) ;
2640 Double_t zterm = 1 + TMath::Cos( k*z ) ;
2641 Double_t qwe = Numerator / Denominator[n] ;
2642 IntegralOverZ += Cn * zterm * qwe ;
2643 if ( n>10 && fabs(IntegralOverZ)*1.e-10 > fabs(qwe) )
break;
2645 if ( eZList[i] < 0 ) IntegralOverZ = -1 * IntegralOverZ ;
2646 shiftEr[i][j] = IntegralOverZ ; }
2650 if (usingCartesian) Cart2Polar(x,r,phi);
2651 else { r = x[0]; phi = x[1]; }
2652 if ( phi < 0 ) phi += TMath::TwoPi() ;
2653 z = LimitZ( Sector, x ) ;
2655 Interpolate2DEdistortion( ORDER, r, z, shiftEr, Er_integral ) ;
2656 Ephi_integral = 0.0 ;
2661 phi = phi - ( Const_0*Ephi_integral - Const_1*Er_integral ) / r ;
2662 r = r - ( Const_0*Er_integral + Const_1*Ephi_integral ) ;
2665 if (usingCartesian) Polar2Cart(r,phi,Xprime);
2666 else { Xprime[0] = r; Xprime[1] = phi; }
2682 int active_options = mDistortionMode & (kSpaceCharge | kSpaceChargeR2 | kSpaceChargeFXT);
2683 if (active_options & (active_options-1)) {
2684 cout <<
"StMagUtilities ERROR **** Do not use multiple SpaceCharge modes at the same time" << endl ;
2685 cout <<
"StMagUtilities ERROR **** These routines have overlapping functionality." << endl ;
2689 if (mDistortionMode & kSpaceCharge) {
2691 }
else if (mDistortionMode & kSpaceChargeR2) {
2693 }
else if (mDistortionMode & kSpaceChargeFXT) {
2713 Float_t Er_integral, Ephi_integral ;
2714 Double_t r, phi, z ;
2716 const Int_t ORDER = 1 ;
2720 Int_t Nterms = 100 ;
2721 Double_t Denominator[100];
2722 memset(Denominator,0,100*
sizeof(Double_t));
2723 for ( Int_t i = 0 ; i < EMap_nZ ; ++i )
2725 z = TMath::Abs( eZList[i] ) ;
2726 for ( Int_t j = 0 ; j < EMap_nR ; ++j )
2729 spaceEr[i][j] = 0.0 ;
2730 if (r < IFCRadius)
continue;
2731 if (r > OFCRadius)
continue;
2732 if (z > TPC_Z0)
continue;
2733 Double_t IntegralOverZ = 0.0 ;
2734 for ( Int_t n = 1 ; n < Nterms ; ++n )
2736 Double_t k = n * TMath::Pi() / TPC_Z0 ;
2737 Double_t zterm = TMath::Power(-1,(n+1)) * ( 1.0 - TMath::Cos( k * ( TPC_Z0 - z ) ) ) ;
2740 Double_t Cn = -4.0 / ( k*k*k * TPC_Z0 * StarMagE ) ;
2741 Double_t Numerator =
2742 TMath::BesselI1( k*r ) * TMath::BesselK0( k*OFCRadius ) -
2743 TMath::BesselI1( k*r ) * TMath::BesselK0( k*IFCRadius ) +
2744 TMath::BesselK1( k*r ) * TMath::BesselI0( k*OFCRadius ) -
2745 TMath::BesselK1( k*r ) * TMath::BesselI0( k*IFCRadius ) ;
2746 if (Denominator[n] == 0) Denominator[n] =
2747 TMath::BesselK0( k*OFCRadius ) * TMath::BesselI0( k*IFCRadius ) -
2748 TMath::BesselK0( k*IFCRadius ) * TMath::BesselI0( k*OFCRadius ) ;
2749 Double_t qwe = Numerator / Denominator[n] ;
2750 IntegralOverZ += Cn * zterm * qwe ;
2751 if ( n>10 && fabs(IntegralOverZ)*1.e-10 > fabs(qwe) )
break;
2753 spaceEr[i][j] = IntegralOverZ ;
2758 if (usingCartesian) Cart2Polar(x,r,phi);
2759 else { r = x[0]; phi = x[1]; }
2760 if ( phi < 0 ) phi += TMath::TwoPi() ;
2761 z = LimitZ( Sector, x ) ;
2763 Interpolate2DEdistortion( ORDER, r, z, spaceEr, Er_integral ) ;
2764 Ephi_integral = 0.0 ;
2768 if (fSpaceCharge) GetSpaceCharge();
2773 Float_t Weight = SpaceCharge * (doingDistortion ? SmearCoefSC : 1.0);
2774 phi = phi - Weight * ( Const_0*Ephi_integral - Const_1*Er_integral ) / r ;
2775 r = r - Weight * ( Const_0*Er_integral + Const_1*Ephi_integral ) ;
2778 if (usingCartesian) Polar2Cart(r,phi,Xprime);
2779 else { Xprime[0] = r; Xprime[1] = phi; }
2805 const Int_t ORDER = 1 ;
2807 Float_t Er_integral, Ephi_integral ;
2808 Double_t r, phi, z ;
2810 if (fSpaceChargeR2) { GetSpaceChargeR2();}
2814 cout <<
"StMagUtilities::UndoSpace Please wait for the tables to fill ... ~5 seconds" << endl ;
2815 const Int_t ROWS = 257 ;
2816 const Int_t COLUMNS = 129 ;
2817 const Int_t ITERATIONS = 100 ;
2818 const Double_t GRIDSIZER = (OFCRadius-IFCRadius) / (ROWS-1) ;
2819 const Double_t GRIDSIZEZ = TPC_Z0 / (COLUMNS-1) ;
2820 TMatrix ArrayV(ROWS,COLUMNS), Charge(ROWS,COLUMNS) ;
2821 TMatrix EroverEz(ROWS,COLUMNS) ;
2822 Float_t Rlist[ROWS], Zedlist[COLUMNS] ;
2825 for ( Int_t j = 0 ; j < COLUMNS ; j++ )
2827 Double_t zed = j*GRIDSIZEZ ;
2829 for ( Int_t i = 0 ; i < ROWS ; i++ )
2831 Double_t Radius = IFCRadius + i*GRIDSIZER ;
2838 for ( Int_t j = 1 ; j < COLUMNS-1 ; j++ )
2840 Double_t zed = j*GRIDSIZEZ ;
2841 for ( Int_t i = 1 ; i < ROWS-1 ; i++ )
2843 Double_t Radius = IFCRadius + i*GRIDSIZER ;
2844 Double_t zterm = (TPC_Z0-zed) * (OFCRadius*OFCRadius - IFCRadius*IFCRadius) / TPC_Z0 ;
2857 Charge(i,j) = zterm * SpaceChargeRadialDependence(Radius) ;
2873 PoissonRelaxation( ArrayV, Charge, EroverEz, ITERATIONS ) ;
2876 Int_t ilow=0, jlow=0 ;
2877 Float_t save_Er[2] ;
2878 for ( Int_t i = 0 ; i < EMap_nZ ; ++i )
2880 z = TMath::Abs( eZList[i] ) ;
2881 for ( Int_t j = 0 ; j < EMap_nR ; ++j )
2884 Search( ROWS, Rlist, r, ilow ) ;
2885 Search( COLUMNS, Zedlist, z, jlow ) ;
2886 if ( ilow < 0 ) ilow = 0 ;
2887 if ( jlow < 0 ) jlow = 0 ;
2888 if ( ilow + 1 >= ROWS - 1 ) ilow = ROWS - 2 ;
2889 if ( jlow + 1 >= COLUMNS - 1 ) jlow = COLUMNS - 2 ;
2890 save_Er[0] = EroverEz(ilow,jlow) + (EroverEz(ilow,jlow+1)-EroverEz(ilow,jlow))*(z-Zedlist[jlow])/GRIDSIZEZ ;
2891 save_Er[1] = EroverEz(ilow+1,jlow) + (EroverEz(ilow+1,jlow+1)-EroverEz(ilow+1,jlow))*(z-Zedlist[jlow])/GRIDSIZEZ ;
2892 spaceR2Er[i][j] = save_Er[0] + (save_Er[1]-save_Er[0])*(r-Rlist[ilow])/GRIDSIZER ;
2898 if (usingCartesian) Cart2Polar(x,r,phi);
2899 else { r = x[0]; phi = x[1]; }
2900 if ( phi < 0 ) phi += TMath::TwoPi() ;
2901 z = LimitZ( Sector, x ) ;
2903 Interpolate2DEdistortion( ORDER, r, z, spaceR2Er, Er_integral ) ;
2904 Ephi_integral = 0.0 ;
2913 double Weight = SpaceChargeR2 * (doingDistortion ? SmearCoefSC : 1.0);
2914 if ( z < 0) Weight *= SpaceChargeEWRatio ;
2915 phi = phi - Weight * ( Const_0*Ephi_integral - Const_1*Er_integral ) / r ;
2916 r = r - Weight * ( Const_0*Er_integral + Const_1*Ephi_integral ) ;
2919 if (usingCartesian) Polar2Cart(r,phi,Xprime);
2920 else { Xprime[0] = r; Xprime[1] = phi; }
2938 const Int_t ORDER = 1 ;
2940 Float_t Er_integral, Ephi_integral ;
2941 Double_t r, phi, z ;
2943 if (fSpaceChargeR2) { GetSpaceChargeR2();}
2947 cout <<
"StMagUtilities::UndoSpace Please wait for the tables to fill ... ~10 seconds" << endl ;
2948 const Int_t ROWS = 257 ;
2949 const Int_t COLUMNS = 129 ;
2950 const Int_t ITERATIONS = 100 ;
2951 const Double_t GRIDSIZER = (OFCRadius-IFCRadius) / (ROWS-1) ;
2952 const Double_t GRIDSIZEZ = TPC_Z0 / (COLUMNS-1) ;
2953 TMatrix ArrayVE(ROWS,COLUMNS), ChargeE(ROWS,COLUMNS) ;
2954 TMatrix EroverEzE(ROWS,COLUMNS) ;
2955 TMatrix ArrayVW(ROWS,COLUMNS), ChargeW(ROWS,COLUMNS) ;
2956 TMatrix EroverEzW(ROWS,COLUMNS) ;
2957 Float_t Rlist[ROWS], Zedlist[COLUMNS] ;
2960 for ( Int_t j = 0 ; j < COLUMNS ; j++ )
2962 Double_t zed = j*GRIDSIZEZ ;
2964 for ( Int_t i = 0 ; i < ROWS ; i++ )
2966 Double_t Radius = IFCRadius + i*GRIDSIZER ;
2975 for ( Int_t j = 1 ; j < COLUMNS-1 ; j++ )
2977 Double_t zed = j*GRIDSIZEZ ;
2978 for ( Int_t i = 1 ; i < ROWS-1 ; i++ )
2980 Double_t Radius = IFCRadius + i*GRIDSIZER ;
2981 Double_t zterm = (TPC_Z0-zed) * (OFCRadius*OFCRadius - IFCRadius*IFCRadius) / TPC_Z0 ;
2982 ChargeE(i,j) = zterm * SpaceChargeFXTRadialDependenceEast(Radius) ;
2983 ChargeW(i,j) = zterm * SpaceChargeFXTRadialDependenceWest(Radius) ;
2987 PoissonRelaxation( ArrayVE, ChargeE, EroverEzE, ITERATIONS ) ;
2988 PoissonRelaxation( ArrayVW, ChargeW, EroverEzW, ITERATIONS ) ;
2991 Int_t ilow=0, jlow=0 ;
2992 Float_t save_Er[2] ;
2993 for ( Int_t i = 0 ; i < EMap_nZ ; ++i )
2995 z = TMath::Abs( eZList[i] ) ;
2996 TMatrix& EroverEz = ( eZList[i] < 0 ? EroverEzE : EroverEzW );
2997 for ( Int_t j = 0 ; j < EMap_nR ; ++j )
3000 Search( ROWS, Rlist, r, ilow ) ;
3001 Search( COLUMNS, Zedlist, z, jlow ) ;
3002 if ( ilow < 0 ) ilow = 0 ;
3003 if ( jlow < 0 ) jlow = 0 ;
3004 if ( ilow + 1 >= ROWS - 1 ) ilow = ROWS - 2 ;
3005 if ( jlow + 1 >= COLUMNS - 1 ) jlow = COLUMNS - 2 ;
3006 save_Er[0] = EroverEz(ilow,jlow) + (EroverEz(ilow,jlow+1)-EroverEz(ilow,jlow))*(z-Zedlist[jlow])/GRIDSIZEZ ;
3007 save_Er[1] = EroverEz(ilow+1,jlow) + (EroverEz(ilow+1,jlow+1)-EroverEz(ilow+1,jlow))*(z-Zedlist[jlow])/GRIDSIZEZ ;
3008 spaceR2Er[i][j] = save_Er[0] + (save_Er[1]-save_Er[0])*(r-Rlist[ilow])/GRIDSIZER ;
3014 if (usingCartesian) Cart2Polar(x,r,phi);
3015 else { r = x[0]; phi = x[1]; }
3016 if ( phi < 0 ) phi += TMath::TwoPi() ;
3017 z = LimitZ( Sector, x ) ;
3019 Interpolate2DEdistortion( ORDER, r, z, spaceR2Er, Er_integral ) ;
3020 Ephi_integral = 0.0 ;
3025 double Weight = SpaceChargeR2 * (doingDistortion ? SmearCoefSC : 1.0);
3026 if ( z < 0) Weight *= SpaceChargeEWRatio ;
3027 phi = phi - Weight * ( Const_0*Ephi_integral - Const_1*Er_integral ) / r ;
3028 r = r - Weight * ( Const_0*Er_integral + Const_1*Ephi_integral ) ;
3031 if (usingCartesian) Polar2Cart(r,phi,Xprime);
3032 else { Xprime[0] = r; Xprime[1] = phi; }
3080 const Int_t ORDER = 1 ;
3082 static Bool_t DoOnceLocal = true ;
3083 Float_t Er_integral, Ephi_integral ;
3084 Double_t r, phi, z ;
3085 const Int_t TIMEBINS = 20 ;
3087 static Float_t abortR2Er[TIMEBINS][EMap_nZ][EMap_nR] ;
3090 const Int_t ROWS = 257 ;
3091 const Int_t COLUMNS = 129 ;
3095 Float_t AbortGapCharge = 0.0;
3096 Int_t AbortGapChargeSize = 1;
3098 Bool_t testCase = (TimeSinceDeposition >= 0);
3100 if (fAbortGapCharge) GetAbortGapCharge();
3101 AbortGapChargeSize = AbortGapCharges->GetSize();
3102 if (AbortGapChargeSize == 0) { memcpy(Xprime,x,threeFloats); return ; }
3106 Double_t fullDriftTime = TPC_Z0 / IonDriftVel;
3111 cout <<
"StMagUtilities::UndoAbortGap Please wait for the tables to fill ... ~5 seconds" << endl;
3112 const Int_t ITERATIONS = 100 ;
3113 const Double_t GRIDSIZER = (OFCRadius-IFCRadius) / (ROWS-1) ;
3114 const Double_t GRIDSIZEZ = TPC_Z0 / (COLUMNS-1) ;
3115 Float_t Rlist[ROWS], Zedlist[COLUMNS] ;
3117 Double_t timeBinLength = TPC_Z0 / TIMEBINS;
3118 Double_t zterm = OFCRadius*OFCRadius - IFCRadius*IFCRadius ;
3120 TMatrix ArrayV(ROWS,COLUMNS), Charge(ROWS,COLUMNS), EroverEz(ROWS,COLUMNS) ;
3122 for ( Int_t k = 0 ; k < TIMEBINS ; k++ )
3124 Double_t driftZ = k * timeBinLength ;
3126 for ( Int_t j = 0 ; j < COLUMNS ; j++ )
3128 Double_t zed = j*GRIDSIZEZ ;
3130 for ( Int_t i = 0 ; i < ROWS ; i++ )
3132 Double_t Radius = IFCRadius + i*GRIDSIZER ;
3140 for ( Int_t j = 1 ; j < COLUMNS-1 ; j++ )
3142 Double_t zed = j*GRIDSIZEZ ;
3143 if ( zed <= TPC_Z0 - driftZ )
3145 for ( Int_t i = 1 ; i < ROWS-1 ; i++ )
3147 Double_t Radius = IFCRadius + i*GRIDSIZER ;
3148 Charge(i,j) = zterm * SpaceChargeRadialDependence(Radius) ;
3151 for ( Int_t i = 1 ; i < ROWS-1 ; i++ ) Charge(i,j) = 0.0 ;
3156 PoissonRelaxation( ArrayV, Charge, EroverEz, ITERATIONS ) ;
3158 Int_t ilow=0, jlow=0 ;
3159 Float_t save_Er[2] ;
3160 for ( Int_t i = 0 ; i < EMap_nZ ; ++i )
3162 z = TMath::Abs( eZList[i] ) ;
3163 for ( Int_t j = 0 ; j <EMap_nR ; ++j )
3166 Search( ROWS, Rlist, r, ilow ) ;
3167 Search( COLUMNS, Zedlist, z, jlow ) ;
3168 if ( ilow < 0 ) ilow = 0 ;
3169 if ( jlow < 0 ) jlow = 0 ;
3170 if ( ilow + 1 >= ROWS - 1 ) ilow = ROWS - 2 ;
3171 if ( jlow + 1 >= COLUMNS - 1 ) jlow = COLUMNS - 2 ;
3172 save_Er[0] = EroverEz(ilow,jlow) + (EroverEz(ilow,jlow+1)-EroverEz(ilow,jlow))*(z-Zedlist[jlow])/GRIDSIZEZ ;
3173 save_Er[1] = EroverEz(ilow+1,jlow) + (EroverEz(ilow+1,jlow+1)-EroverEz(ilow+1,jlow))*(z-Zedlist[jlow])/GRIDSIZEZ ;
3174 abortR2Er[k][i][j] = save_Er[0] + (save_Er[1]-save_Er[0])*(r-Rlist[ilow])/GRIDSIZER ;
3179 DoOnceLocal = false ;
3182 if (usingCartesian) Cart2Polar(x,r,phi);
3183 else { r = x[0]; phi = x[1]; }
3184 if ( phi < 0 ) phi += TMath::TwoPi() ;
3185 z = LimitZ( Sector, x ) ;
3187 Double_t timeBinDuration = fullDriftTime / TIMEBINS;
3189 for ( Int_t i=0 ; i<AbortGapChargeSize ; i++)
3193 if (fSpaceChargeR2) { GetSpaceChargeR2();}
3194 AbortGapCharge = SpaceChargeR2;
3196 TimeSinceDeposition = (Float_t) (*AbortGapTimes)[i];
3197 if ( TimeSinceDeposition >= fullDriftTime )
continue;
3198 AbortGapCharge = AbortGapChargeCoef * (*AbortGapCharges)[i];
3202 Int_t timeSlice = TIMEBINS;
3203 Double_t closestTime = 25.0 ;
3205 for ( Int_t k = 0 ; k < TIMEBINS ; k++ )
3207 Double_t timeSliceTime = k * timeBinDuration;
3208 Double_t timeDiff = abs(TimeSinceDeposition - timeSliceTime ) ;
3209 if ( timeDiff < abs(TimeSinceDeposition - closestTime) )
3211 closestTime = timeSliceTime;
3216 if ( timeSlice >=0 && timeSlice < TIMEBINS)
3219 Interpolate2DEdistortion( ORDER, r, z, abortR2Er[timeSlice], Er_integral ) ;
3220 Ephi_integral = 0.0 ;
3224 double Weight = AbortGapCharge * (doingDistortion ? SmearCoefSC : 1.0);
3225 phi = phi - Weight * ( Const_0*Ephi_integral - Const_1*Er_integral ) / r ;
3226 r = r - Weight * ( Const_0*Er_integral + Const_1*Ephi_integral ) ;
3231 if (usingCartesian) Polar2Cart(r,phi,Xprime);
3232 else { Xprime[0] = r; Xprime[1] = phi; }
3267 const Int_t ORDER = 1 ;
3269 static Bool_t DoOnceLocal = true ;
3270 static Int_t NumberOfEastInnerShorts = 0, NumberOfEastOuterShorts = 0 , NumberOfWestInnerShorts = 0, NumberOfWestOuterShorts = 0 ;
3272 Float_t Er_integral, Ephi_integral ;
3273 Double_t r, phi, z ;
3275 if (fTpcVolts) {DoOnceLocal = UpdateShortedRing() ;}
3280 cout <<
"StMagUtilities::UndoShort Please wait for the tables to fill ... ~5 seconds" << endl ;
3284 const Float_t Z01 = 1.225 ;
3286 Float_t EastInnerMissingSum = 0, EastOuterMissingSum = 0, WestInnerMissingSum = 0, WestOuterMissingSum = 0 ;
3287 Float_t EastInnerExtraSum = 0, EastOuterExtraSum = 0, WestInnerExtraSum = 0, WestOuterExtraSum = 0 ;
3288 Float_t EastInnerMissingOhms[10], EastOuterMissingOhms[10], WestInnerMissingOhms[10], WestOuterMissingOhms[10] ;
3289 Float_t EastInnerShortZ[10], EastOuterShortZ[10], WestInnerShortZ[10], WestOuterShortZ[10] ;
3291 NumberOfEastInnerShorts = 0; NumberOfEastOuterShorts = 0 ; NumberOfWestInnerShorts = 0; NumberOfWestOuterShorts = 0 ;
3293 for ( Int_t i = 0 ; i < ShortTableRows ; i++ )
3295 if ( ( Side[i] + Cage[i] + Ring[i] + TMath::Abs(MissingResistance[i]) + TMath::Abs(Resistor[i]) ) == 0.0 ) continue ;
3296 if ( Side[i] == 0 && Cage[i] == 0 )
3297 { NumberOfEastInnerShorts++ ; EastInnerMissingSum += MissingResistance[i] ; EastInnerExtraSum += Resistor[i] ;
3298 EastInnerMissingOhms[NumberOfEastInnerShorts-1] = MissingResistance[i] ;
3299 EastInnerShortZ[NumberOfEastInnerShorts-1] = TPC_Z0 - ( Z01 + (Ring[i]-1)*RPitch) ; }
3300 if ( Side[i] == 0 && Cage[i] == 1 )
3301 { NumberOfEastOuterShorts++ ; EastOuterMissingSum += MissingResistance[i] ; EastOuterExtraSum += Resistor[i] ;
3302 EastOuterMissingOhms[NumberOfEastOuterShorts-1] = MissingResistance[i] ;
3303 EastOuterShortZ[NumberOfEastOuterShorts-1] = TPC_Z0 - ( Z01 + (Ring[i]-1)*RPitch) ; }
3304 if ( Side[i] == 1 && Cage[i] == 0 )
3305 { NumberOfWestInnerShorts++ ; WestInnerMissingSum += MissingResistance[i] ; WestInnerExtraSum += Resistor[i] ;
3306 WestInnerMissingOhms[NumberOfWestInnerShorts-1] = MissingResistance[i] ;
3307 WestInnerShortZ[NumberOfWestInnerShorts-1] = TPC_Z0 - ( Z01 + (Ring[i]-1)*RPitch) ; }
3308 if ( Side[i] == 1 && Cage[i] == 1 )
3309 { NumberOfWestOuterShorts++ ; WestOuterMissingSum += MissingResistance[i] ; WestOuterExtraSum += Resistor[i] ;
3310 WestOuterMissingOhms[NumberOfWestOuterShorts-1] = MissingResistance[i] ;
3311 WestOuterShortZ[NumberOfWestOuterShorts-1] = TPC_Z0 - ( Z01 + (Ring[i]-1)*RPitch) ; }
3315 if ( (NumberOfEastInnerShorts + NumberOfEastOuterShorts + NumberOfWestInnerShorts + NumberOfWestOuterShorts) == 0 )
3316 { memcpy(Xprime,x,threeFloats); return ; }
3318 Float_t EastInnerRtot = Rtot + EastInnerExtraSum - EastInnerMissingSum ;
3319 Float_t EastOuterRtot = Rtot + EastOuterExtraSum - EastOuterMissingSum ;
3320 Float_t WestInnerRtot = Rtot + WestInnerExtraSum - WestInnerMissingSum ;
3321 Float_t WestOuterRtot = Rtot + WestOuterExtraSum - WestOuterMissingSum ;
3326 Int_t Nterms = 100 ;
3327 Double_t Denominator[100];
3328 memset(Denominator,0,100*
sizeof(Double_t));
3329 for ( Int_t i = 0 ; i < EMap_nZ ; ++i )
3332 for ( Int_t j = 0 ; j < EMap_nR ; ++j )
3335 shortEr[i][j] = 0.0 ;
3336 if (r < IFCRadius)
continue;
3337 if (r > OFCRadius)
continue;
3338 if (TMath::Abs(z) > TPC_Z0)
continue;
3339 Double_t IntegralOverZ = 0.0 ;
3340 for ( Int_t n = 1 ; n < Nterms ; ++n )
3342 Double_t k = n * TMath::Pi() / TPC_Z0 ;
3349 for ( Int_t m = 0 ; m < NumberOfEastInnerShorts ; m++ )
3350 sum += ( 1 - Rfrac - TMath::Cos(k*EastInnerShortZ[m]) ) * EastInnerMissingOhms[m] ;
3351 Ein = 2 * ( Rfrac*EastInnerExtraSum + sum ) / (k*EastInnerRtot) ;
3353 for ( Int_t m = 0 ; m < NumberOfEastOuterShorts ; m++ )
3354 sum += ( 1 - Rfrac - TMath::Cos(k*EastOuterShortZ[m]) ) * EastOuterMissingOhms[m] ;
3355 Eout = 2 * ( Rfrac*EastOuterExtraSum + sum ) / (k*EastOuterRtot) ;
3357 if ( z == 0 ) continue ;
3361 for ( Int_t m = 0 ; m < NumberOfWestInnerShorts ; m++ )
3362 sum += ( 1 - Rfrac - TMath::Cos(k*WestInnerShortZ[m]) ) * WestInnerMissingOhms[m] ;
3363 Ein = 2 * ( Rfrac*WestInnerExtraSum + sum ) / (k*WestInnerRtot) ;
3365 for ( Int_t m = 0 ; m < NumberOfWestOuterShorts ; m++ )
3366 sum += ( 1 - Rfrac - TMath::Cos(k*WestOuterShortZ[m]) ) * WestOuterMissingOhms[m] ;
3367 Eout = 2 * ( Rfrac*WestOuterExtraSum + sum ) / (k*WestOuterRtot) ;
3377 Double_t An = Ein * TMath::BesselK0( k*OFCRadius ) - Eout * TMath::BesselK0( k*IFCRadius ) ;
3378 Double_t Bn = Eout * TMath::BesselI0( k*IFCRadius ) - Ein * TMath::BesselI0( k*OFCRadius ) ;
3379 Double_t Numerator =
3380 An * TMath::BesselI1( k*r ) - Bn * TMath::BesselK1( k*r ) ;
3381 if (Denominator[n] == 0) Denominator[n] =
3382 TMath::BesselK0( k*OFCRadius ) * TMath::BesselI0( k*IFCRadius ) -
3383 TMath::BesselK0( k*IFCRadius ) * TMath::BesselI0( k*OFCRadius ) ;
3384 Double_t zterm = TMath::Cos( k*(TPC_Z0-TMath::Abs(z)) ) - 1 ;
3385 Double_t qwe = Numerator / Denominator[n] ;
3386 IntegralOverZ += zterm * qwe ;
3387 if ( n>10 && fabs(IntegralOverZ)*1.e-10 > fabs(qwe) )
break;
3389 shortEr[i][j] = IntegralOverZ ;
3392 DoOnceLocal = false ;
3395 if ( SectorSide( Sector, x ) > 0 ) {
3396 if ( (NumberOfWestInnerShorts + NumberOfWestOuterShorts) == 0 )
3397 { memcpy(Xprime,x,threeFloats); return ; }
3398 }
else if ( (NumberOfEastInnerShorts + NumberOfEastOuterShorts) == 0 )
3399 { memcpy(Xprime,x,threeFloats); return ; }
3401 if (usingCartesian) Cart2Polar(x,r,phi);
3402 else { r = x[0]; phi = x[1]; }
3403 if ( phi < 0 ) phi += TMath::TwoPi() ;
3404 z = LimitZ( Sector, x ) ;
3406 Interpolate2DEdistortion( ORDER, r, z, shortEr, Er_integral ) ;
3407 Ephi_integral = 0.0 ;
3412 phi = phi - ( Const_0*Ephi_integral - Const_1*Er_integral ) / r ;
3413 r = r - ( Const_0*Er_integral + Const_1*Ephi_integral ) ;
3416 if (usingCartesian) Polar2Cart(r,phi,Xprime);
3417 else { Xprime[0] = r; Xprime[1] = phi; }
3439 const Int_t ORDER = 1 ;
3441 Float_t Er_integral, Ephi_integral ;
3442 Double_t r, phi, z ;
3449 cout <<
"StMagUtilities::UndoGG VE Please wait for the tables to fill ... ~5 seconds" << endl ;
3451 Int_t Nterms = 100 ;
3452 Double_t Denominator[100];
3453 memset(Denominator,0,100*
sizeof(Double_t));
3454 for ( Int_t i = 0 ; i < EMap_nZ ; ++i )
3457 for ( Int_t j = 0 ; j < EMap_nR ; ++j )
3460 GGVoltErrorEr[i][j] = 0.0 ;
3461 if (r < IFCRadius)
continue;
3462 if (r > OFCRadius)
continue;
3463 if (TMath::Abs(z) > TPC_Z0)
continue;
3464 Double_t IntegralOverZ = 0.0 ;
3465 for ( Int_t n = 1 ; n < Nterms ; ++n )
3467 if ( z == 0 ) continue ;
3468 Double_t k = n * TMath::Pi() / TPC_Z0 ;
3469 Double_t Ein = -2.0 * (z < 0 ? deltaVGGEast : deltaVGGWest) * deltaGGeffectiveness /
3470 (k * (CathodeV - GGideal));
3471 Double_t Eout = Ein ;
3473 Double_t An = Ein * TMath::BesselK0( k*OFCRadius ) - Eout * TMath::BesselK0( k*IFCRadius ) ;
3474 Double_t Bn = Eout * TMath::BesselI0( k*IFCRadius ) - Ein * TMath::BesselI0( k*OFCRadius ) ;
3475 Double_t Numerator =
3476 An * TMath::BesselI1( k*r ) - Bn * TMath::BesselK1( k*r ) ;
3477 if (Denominator[n] == 0) Denominator[n] =
3478 TMath::BesselK0( k*OFCRadius ) * TMath::BesselI0( k*IFCRadius ) -
3479 TMath::BesselK0( k*IFCRadius ) * TMath::BesselI0( k*OFCRadius ) ;
3480 Double_t zterm = TMath::Cos( k*(TPC_Z0-TMath::Abs(z)) ) - 1 ;
3481 Double_t qwe = Numerator / Denominator[n] ;
3482 IntegralOverZ += zterm * qwe ;
3483 if ( n>10 && fabs(IntegralOverZ)*1.e-10 > fabs(qwe) )
break;
3485 GGVoltErrorEr[i][j] = IntegralOverZ ;
3490 if ( deltaVGGEast == 0.0 && deltaVGGWest == 0 )
3491 { memcpy(Xprime,x,threeFloats) ; return ; }
3493 if (usingCartesian) Cart2Polar(x,r,phi);
3494 else { r = x[0]; phi = x[1]; }
3495 if ( phi < 0 ) phi += TMath::TwoPi() ;
3496 z = LimitZ( Sector, x ) ;
3498 Interpolate2DEdistortion( ORDER, r, z, GGVoltErrorEr, Er_integral ) ;
3499 Ephi_integral = 0.0 ;
3504 phi = phi - ( Const_0*Ephi_integral - Const_1*Er_integral ) / r ;
3505 r = r - ( Const_0*Er_integral + Const_1*Ephi_integral ) ;
3508 if (usingCartesian) Polar2Cart(r,phi,Xprime);
3509 else { Xprime[0] = r; Xprime[1] = phi; }
3517 Float_t StMagUtilities::Interpolate2DTable(
const Int_t ORDER,
const Float_t x,
const Float_t y,
const Int_t nx,
const Int_t ny,
3518 const Float_t XV[],
const Float_t YV[],
const TMatrix &Array )
3521 static Int_t jlow = 0, klow = 0 ;
3522 Float_t save_Array[ORDER+1] ;
3524 Search( nx, XV, x, jlow ) ;
3525 Search( ny, YV, y, klow ) ;
3526 if ( jlow < 0 ) jlow = 0 ;
3527 if ( klow < 0 ) klow = 0 ;
3528 if ( jlow + ORDER >= nx - 1 ) jlow = nx - 1 - ORDER ;
3529 if ( klow + ORDER >= ny - 1 ) klow = ny - 1 - ORDER ;
3531 for ( Int_t j = jlow ; j < jlow + ORDER + 1 ; j++ )
3533 Float_t *ajkl = &((TMatrix&)Array)(j,klow);
3534 save_Array[j-jlow] = Interpolate( &YV[klow], ajkl , ORDER, y ) ;
3537 return( Interpolate( &XV[jlow], save_Array, ORDER, x ) ) ;
3542 void StMagUtilities::Interpolate2DEdistortion(
const Int_t ORDER,
const Float_t r,
const Float_t z,
3543 const Float_t Er[EMap_nZ][EMap_nR], Float_t &Er_value )
3546 static Int_t jlow = 0, klow = 0 ;
3547 Float_t save_Er[ORDER+1] ;
3549 Search( EMap_nZ, eZList, z, jlow ) ;
3550 Search( EMap_nR, eRList, r, klow ) ;
3551 if ( jlow < 0 ) jlow = 0 ;
3552 if ( klow < 0 ) klow = 0 ;
3553 if ( jlow + ORDER >= EMap_nZ - 1 ) jlow = EMap_nZ - 1 - ORDER ;
3554 if ( klow + ORDER >= EMap_nR - 1 ) klow = EMap_nR - 1 - ORDER ;
3556 for ( Int_t j = jlow ; j < jlow + ORDER + 1 ; j++ )
3558 save_Er[j-jlow] = Interpolate( &eRList[klow], &Er[j][klow], ORDER, r ) ;
3560 Er_value = Interpolate( &eZList[jlow], save_Er, ORDER, z ) ;
3565 void StMagUtilities::Interpolate3DEdistortion(
const Int_t ORDER,
const Float_t r,
const Float_t phi,
const Float_t z,
3566 const Float_t Er[EMap_nZ][EMap_nPhi][EMap_nR],
const Float_t Ephi[EMap_nZ][EMap_nPhi][EMap_nR],
3567 Float_t &Er_value, Float_t &Ephi_value )
3570 static Int_t ilow = 0, jlow = 0, klow = 0 ;
3571 Float_t save_Er[ORDER+1], saved_Er[ORDER+1] ;
3572 Float_t save_Ephi[ORDER+1], saved_Ephi[ORDER+1] ;
3574 Search( EMap_nZ, eZList, z, ilow ) ;
3575 Search( EMap_nPhi, ePhiList, phi, jlow ) ;
3576 Search( EMap_nR, eRList, r, klow ) ;
3577 if ( ilow < 0 ) ilow = 0 ;
3578 if ( jlow < 0 ) jlow = 0 ;
3579 if ( klow < 0 ) klow = 0 ;
3581 if ( ilow + ORDER >= EMap_nZ - 1 ) ilow = EMap_nZ - 1 - ORDER ;
3582 if ( jlow + ORDER >= EMap_nPhi - 1 ) jlow = EMap_nPhi - 1 - ORDER ;
3583 if ( klow + ORDER >= EMap_nR - 1 ) klow = EMap_nR - 1 - ORDER ;
3585 for ( Int_t i = ilow ; i < ilow + ORDER + 1 ; i++ )
3587 for ( Int_t j = jlow ; j < jlow + ORDER + 1 ; j++ )
3589 save_Er[j-jlow] = Interpolate( &eRList[klow], &Er[i][j][klow], ORDER, r ) ;
3590 save_Ephi[j-jlow] = Interpolate( &eRList[klow], &Ephi[i][j][klow], ORDER, r ) ;
3592 saved_Er[i-ilow] = Interpolate( &ePhiList[jlow], save_Er, ORDER, phi ) ;
3593 saved_Ephi[i-ilow] = Interpolate( &ePhiList[jlow], save_Ephi, ORDER, phi ) ;
3595 Er_value = Interpolate( &eZList[ilow], saved_Er, ORDER, z ) ;
3596 Ephi_value = Interpolate( &eZList[ilow], saved_Ephi, ORDER, z ) ;
3602 Float_t StMagUtilities::Interpolate3DTable (
const Int_t ORDER,
const Float_t x,
const Float_t y,
const Float_t z,
3603 const Int_t nx,
const Int_t ny,
const Int_t nz,
3604 const Float_t XV[],
const Float_t YV[],
const Float_t ZV[],
3605 TMatrix **ArrayofArrays )
3608 static Int_t ilow = 0, jlow = 0, klow = 0 ;
3609 Float_t save_Array[ORDER+1], saved_Array[ORDER+1] ;
3611 Search( nx, XV, x, ilow ) ;
3612 Search( ny, YV, y, jlow ) ;
3613 Search( nz, ZV, z, klow ) ;
3615 if ( ilow < 0 ) ilow = 0 ;
3616 if ( jlow < 0 ) jlow = 0 ;
3617 if ( klow < 0 ) klow = 0 ;
3619 if ( ilow + ORDER >= nx - 1 ) ilow = nx - 1 - ORDER ;
3620 if ( jlow + ORDER >= ny - 1 ) jlow = ny - 1 - ORDER ;
3621 if ( klow + ORDER >= nz - 1 ) klow = nz - 1 - ORDER ;
3623 for ( Int_t k = klow ; k < klow + ORDER + 1 ; k++ )
3625 TMatrix &Table = *ArrayofArrays[k] ;
3626 for ( Int_t i = ilow ; i < ilow + ORDER + 1 ; i++ )
3628 save_Array[i-ilow] = Interpolate( &YV[jlow], &Table(i,jlow), ORDER, y ) ;
3630 saved_Array[k-klow] = Interpolate( &XV[ilow], save_Array, ORDER, x ) ;
3632 return( Interpolate( &ZV[klow], saved_Array, ORDER, z ) ) ;
3659 void StMagUtilities::PoissonRelaxation( TMatrix &ArrayVM, TMatrix &ChargeM, TMatrix &EroverEzM,
3660 const Int_t ITERATIONS )
3663 const Int_t ROWS = ArrayVM.GetNrows() ;
3664 const Int_t COLUMNS = ArrayVM.GetNcols() ;
3665 const Float_t GRIDSIZER = (OFCRadius-IFCRadius) / (ROWS-1) ;
3666 const Float_t GRIDSIZEZ = TPC_Z0 / (COLUMNS-1) ;
3667 const Float_t Ratio = GRIDSIZER*GRIDSIZER / (GRIDSIZEZ*GRIDSIZEZ) ;
3672 if ( !IsPowerOfTwo(ROWS-1) )
3673 { cout <<
"StMagUtilities::PoissonRelaxation - Error in the number of ROWS. Must be 2**M - 1" << endl ; exit(1) ; }
3674 if ( !IsPowerOfTwo(COLUMNS-1) )
3675 { cout <<
"StMagUtilities::PoissonRelaxation - Error in the number of COLUMNS. Must be 2**N - 1" << endl ; exit(1) ; }
3678 Float_t *ArrayE,*ArrayV,*
Charge,*SumCharge,*EroverEz ;
3680 TMatrix ArrayEM(ROWS,COLUMNS) ;
3681 TMatrix SumChargeM(ROWS,COLUMNS) ;
3682 ArrayE = ArrayEM.GetMatrixArray() ;
3683 SumCharge = SumChargeM.GetMatrixArray() ;
3684 ArrayV = ArrayVM.GetMatrixArray() ;
3685 Charge = ChargeM.GetMatrixArray() ;
3686 EroverEz = EroverEzM.GetMatrixArray() ;
3692 Int_t i_one = (ROWS-1)/4 ;
3693 Int_t j_one = (COLUMNS-1)/4 ;
3694 Int_t loops = 1 + (int) ( 0.5 + TMath::Log2( (
double) TMath::Max(i_one,j_one) ) ) ;
3696 Float_t coef1[ROWS],coef2[ROWS];
3697 memset(coef1,0,ROWS*
sizeof(Float_t));
3698 memset(coef2,0,ROWS*
sizeof(Float_t));
3700 for ( Int_t count = 0 ; count < loops ; count++ ) {
3706 Int_t one__ = i_one*COLUMNS;
3707 Int_t half__ = one__ / 2;
3708 Int_t half = j_one / 2;
3710 Float_t tempGRIDSIZER = GRIDSIZER * i_one ;
3711 Float_t tempRatio = Ratio * i_one * i_one / ( j_one * j_one ) ;
3712 Float_t tempFourth = 1.0 / (2.0 + 2.0*tempRatio) ;
3714 for ( Int_t i = i_one ; i < ROWS-1 ; i+=i_one ) {
3715 Float_t Radius = IFCRadius + i*GRIDSIZER ;
3716 coef1[i] = 1.0 + tempGRIDSIZER/(2*Radius);
3717 coef2[i] = 1.0 - tempGRIDSIZER/(2*Radius);
3720 for ( Int_t i = i_one ; i < ROWS-1 ; i += i_one ) {
3721 Int_t i__ = i*COLUMNS;
3722 Float_t Radius = IFCRadius + i*GRIDSIZER ;
3723 for ( Int_t j = j_one ; j < COLUMNS-1 ; j += j_one ) {
3724 Int_t i__j = i__ + j;
3725 if ( i_one == 1 && j_one == 1 ) SumCharge[i__j] = Charge[i__j] ;
3727 Float_t weight = 0.0 ;
3729 SumCharge[i__j]= 0.0 ;
3730 for ( Int_t ii = i-i_one/2 ; ii <= i+i_one/2 ; ii++ ) {
3731 for ( Int_t jj = j-j_one/2 ; jj <= j+j_one/2 ; jj++ ) {
3732 if ( ii == i-i_one/2 || ii == i+i_one/2 || jj == j-j_one/2 || jj == j+j_one/2 ) weight = 0.5*Radius ;
3733 else weight = Radius ;
3734 SumCharge[i__j] += Charge[ii*COLUMNS+jj]*weight ;
3738 SumCharge[i__j] /= sum ;
3740 SumCharge[i__j] *= tempGRIDSIZER*tempGRIDSIZER;
3744 for ( Int_t k = 1 ; k <= ITERATIONS; k++ ) {
3747 Float_t OverRelaxM1 = TMath::Sqrt( TMath::Cos( (k*TMath::PiOver2())/ITERATIONS ) ) ;
3748 Float_t OverRelaxtempFourth = (1.0 + OverRelaxM1) * tempFourth ;
3750 for ( Int_t i = i_one ; i < ROWS-1 ; i += i_one ) {
3751 Int_t i__ = i*COLUMNS;
3752 for ( Int_t j = j_one ; j < COLUMNS-1 ; j += j_one ) {
3753 Int_t i__j = i__ + j;
3755 ArrayV[i__j] = ( coef2[i] * ArrayV[i__j - one__]
3756 + tempRatio * ( ArrayV[i__j - j_one] + ArrayV[i__j + j_one] )
3757 + coef1[i] * ArrayV[i__j + one__]
3759 ) * OverRelaxtempFourth
3760 - OverRelaxM1 * ArrayV[i__j] ;
3767 if ( i_one > 1 || j_one > 1 ) {
3768 for ( Int_t i = i_one ; i < ROWS-1 ; i += i_one ) {
3769 Int_t i__ = i*COLUMNS;
3770 for ( Int_t j = j_one ; j < COLUMNS-1 ; j += j_one ) {
3771 Int_t i__j = i__ + j;
3774 ArrayV[i__j + half__] = ( ArrayV[i__j + one__] + ArrayV[i__j] ) / 2 ;
3776 ArrayV[i__j - half__] = ( ArrayV[j] + ArrayV[one__ + j] ) / 2 ;
3779 ArrayV[i__j + half] = ( ArrayV[i__j + j_one] + ArrayV[i__j] ) / 2 ;
3781 ArrayV[i__j - half] = ( ArrayV[i__] + ArrayV[i__ + j_one] ) / 2 ;
3784 ArrayV[i__j + half__ + half] = ( ArrayV[i__j + one__ + j_one] + ArrayV[i__j] ) / 2 ;
3786 ArrayV[i__j - half__ - half] = ( ArrayV[j - j_one] + ArrayV[one__ + j] ) / 2 ;
3788 ArrayV[i__j - half__ - half] = ( ArrayV[i__ - one__] + ArrayV[i__ + j_one] ) / 2 ;
3811 if (i_one > j_one / 2) { i_one = i_one / 2 ;
if ( i_one < 1 ) i_one = 1 ; }
3812 if (j_one > i_one ) { j_one = j_one / 2 ;
if ( j_one < 1 ) j_one = 1 ; }
3816 Float_t coef10,coef12 ;
3817 coef10 = -0.5 / GRIDSIZER ;
3818 coef12 = (GRIDSIZEZ/3.0) / (-1*StarMagE) ;
3820 Int_t one__ = COLUMNS;
3825 for ( Int_t j = COLUMNS-1 ; j >= 0 ; j-- )
3828 for ( Int_t i = 1 ; i < ROWS-1 ; i++ ) {
3829 Int_t i__j = i*COLUMNS + j;
3830 ArrayE[i__j] = coef10 * ( ArrayV[i__j + one__] - ArrayV[i__j - one__] ) ;
3832 Int_t r__j = (ROWS-1)*one__ + j;
3833 ArrayE[j] = coef10 * ( - ArrayV[2*one__ + j] + 4*ArrayV[one__ + j] - 3*ArrayV[j] ) ;
3834 ArrayE[r__j] = coef10 * ( 3*ArrayV[r__j] - 4*ArrayV[r__j - one__] + ArrayV[r__j - 2*one__] ) ;
3836 for ( Int_t i = 0 ; i < ROWS ; i++ )
3839 Int_t i__ = i*COLUMNS;
3840 Int_t i__j = i__ + j;
3841 Int_t i__c = i__ + COLUMNS-1;
3842 EroverEz[i__j] = 0.0 ;
3843 if ( j == COLUMNS-2 ) EroverEz[i__j] = coef12 * 1.5 * (ArrayE[i__c - 1] + ArrayE[i__c]) ;
3844 else if ( j != COLUMNS-1 )
3846 for ( Int_t k = i__j ; k <= i__c ; k++ )
3848 EroverEz[i__j] += Index*ArrayE[k] ;
3849 if ( Index != 4 ) Index = 4;
else Index = 2 ;
3851 if ( Index == 4 ) EroverEz[i__j] -= ArrayE[i__c] ;
3852 else if ( Index == 2 ) EroverEz[i__j] += (0.5*ArrayE[i__c - 1] - 2.5*ArrayE[i__c]) ;
3853 EroverEz[i__j] *= coef12 ;
3879 void StMagUtilities::Poisson3DRelaxation( TMatrix **ArrayofArrayV, TMatrix **ArrayofCharge, TMatrix **ArrayofEroverEz,
3880 TMatrix **ArrayofEPhioverEz,
3881 const Int_t PHISLICES,
const Float_t DELTAPHI,
3882 const Int_t ITERATIONS,
const Int_t SYMMETRY )
3885 const Int_t ROWS = ArrayofArrayV[0]->GetNrows();
3886 const Int_t COLUMNS = ArrayofArrayV[0]->GetNcols();
3887 const Float_t GRIDSIZEPHI = DELTAPHI ;
3888 const Float_t GRIDSIZER = (OFCRadius-IFCRadius) / (ROWS-1) ;
3889 const Float_t GRIDSIZEZ = TPC_Z0 / (COLUMNS-1) ;
3890 const Float_t RatioPhi = GRIDSIZER*GRIDSIZER / (GRIDSIZEPHI*GRIDSIZEPHI) ;
3891 const Float_t RatioZ = GRIDSIZER*GRIDSIZER / (GRIDSIZEZ*GRIDSIZEZ) ;
3894 if ( !IsPowerOfTwo((ROWS-1)) )
3895 { cout <<
"StMagUtilities::Poisson3DRelaxation - Error in the number of ROWS. Must be 2**M - 1" << endl ; exit(1) ; }
3896 if ( !IsPowerOfTwo((COLUMNS-1)) )
3897 { cout <<
"StMagUtilities::Poisson3DRelaxation - Error in the number of COLUMNS. Must be 2**N - 1" << endl ; exit(1) ; }
3898 if ( PHISLICES <= 3 )
3899 { cout <<
"StMagUtilities::Poisson3DRelaxation - Error in the number of PHISLICES. Must be larger than 3" << endl ; exit(1) ; }
3902 Float_t *ArrayE,*ArrayV,*ArrayVM,*ArrayVP,*Charge,*SumCharge,*EroverEz,*EPhioverEz ;
3904 TMatrix ArrayEM(ROWS,COLUMNS) ;
3905 ArrayE = ArrayEM.GetMatrixArray() ;
3911 Int_t loops, m_plus, m_minus ;
3912 Int_t i_one = (ROWS-1)/4 ;
3913 Int_t j_one = (COLUMNS-1)/4 ;
3914 loops = TMath::Max(i_one, j_one) ;
3915 loops = 1 + (int) ( 0.5 + TMath::Log2((
double)loops) ) ;
3917 TMatrix *ArrayofSumCharge[1000] ;
3918 if ( PHISLICES > 1000 ) { cout <<
"StMagUtilities::Poisson3D PHISLICES > 1000 is not allowed (nor wise) " << endl ; exit(1) ; }
3919 for ( Int_t i = 0 ; i < PHISLICES ; i++ ) { ArrayofSumCharge[i] =
new TMatrix(ROWS,COLUMNS) ; }
3920 Float_t OverRelaxers[ITERATIONS] ;
3921 for ( Int_t k = 1 ; k <= ITERATIONS; k++ ) {
3922 OverRelaxers[k-1] = 1.0 + TMath::Sqrt( TMath::Cos( (k*TMath::PiOver2())/ITERATIONS ) ) ;
3925 Float_t coef1[ROWS],coef2[ROWS],coef3[ROWS],coef4[ROWS],OverRelaxcoef4[ROWS];
3926 memset(coef1,0,ROWS*
sizeof(Float_t));
3927 memset(coef2,0,ROWS*
sizeof(Float_t));
3928 memset(coef3,0,ROWS*
sizeof(Float_t));
3929 memset(coef4,0,ROWS*
sizeof(Float_t));
3930 memset(OverRelaxcoef4,0,ROWS*
sizeof(Float_t));
3932 for ( Int_t count = 0 ; count < loops ; count++ ) {
3938 Int_t one__ = i_one*COLUMNS;
3939 Int_t half__ = one__ / 2;
3940 Int_t half = j_one / 2;
3941 Bool_t iWillDivide = (i_one > 1 && i_one >= j_one / 2) ;
3942 Bool_t jWillDivide = (j_one > 1 && j_one >= i_one / 2) ;
3944 Float_t tempGRIDSIZER = GRIDSIZER * i_one ;
3945 Float_t tempRatioPhi = RatioPhi * i_one * i_one ;
3946 Float_t tempRatioZ = RatioZ * i_one * i_one / ( j_one * j_one ) ;
3948 for ( Int_t i = i_one ; i < ROWS-1 ; i+=i_one ) {
3949 Float_t Radius = IFCRadius + i*GRIDSIZER ;
3950 coef1[i] = 1.0 + tempGRIDSIZER/(2*Radius);
3951 coef2[i] = 1.0 - tempGRIDSIZER/(2*Radius);
3952 coef3[i] = tempRatioPhi/(Radius*Radius);
3953 coef4[i] = 0.5 / (1.0 + tempRatioZ + coef3[i]);
3956 for ( Int_t m = 0 ; m < PHISLICES ; m++ ) {
3957 Charge = ArrayofCharge[m]->GetMatrixArray() ;
3958 SumCharge = ArrayofSumCharge[m]->GetMatrixArray() ;
3959 for ( Int_t i = i_one ; i < ROWS-1 ; i += i_one ) {
3960 Int_t i__ = i*COLUMNS;
3961 Float_t Radius = IFCRadius + i*GRIDSIZER ;
3962 for ( Int_t j = j_one ; j < COLUMNS-1 ; j += j_one ) {
3963 Int_t i__j = i__ + j;
3964 if ( i_one == 1 && j_one == 1 ) SumCharge[i__j] = Charge[i__j] ;
3966 Float_t weight = 0.0 ;
3968 SumCharge[i__j]= 0.0 ;
3969 for ( Int_t ii = i-i_one/2 ; ii <= i+i_one/2 ; ii++ ) {
3970 for ( Int_t jj = j-j_one/2 ; jj <= j+j_one/2 ; jj++ ) {
3971 if ( ii == i-i_one/2 || ii == i+i_one/2 || jj == j-j_one/2 || jj == j+j_one/2 ) weight = 0.5*Radius ;
3972 else weight = Radius ;
3973 SumCharge[i__j] += Charge[ii*COLUMNS+jj]*weight ;
3977 SumCharge[i__j] /= sum ;
3979 SumCharge[i__j] *= tempGRIDSIZER*tempGRIDSIZER;
3984 for ( Int_t k = 1 ; k <= ITERATIONS; k++ ) {
3985 Float_t OverRelax = OverRelaxers[k-1];
3986 Float_t OverRelaxM1 = OverRelax - 1.0 ;
3987 for ( Int_t i = i_one ; i < ROWS-1 ; i+=i_one ) {
3988 OverRelaxcoef4[i] = OverRelax * coef4[i] ;
3992 for ( Int_t m = 0 ; m < PHISLICES ; m++ ) {
3997 if ( m_plus > PHISLICES-1 ) m_plus = PHISLICES - 2 ;
3998 if ( m_minus < 0 ) m_minus = 1 ;
4001 if ( m_plus > PHISLICES-1 ) m_plus = 0 ;
4002 if ( m_minus < 0 ) m_minus = PHISLICES - 1 ;
4004 ArrayV = ArrayofArrayV[m]->GetMatrixArray();
4005 ArrayVP = ArrayofArrayV[m_plus]->GetMatrixArray();
4006 ArrayVM = ArrayofArrayV[m_minus]->GetMatrixArray();
4007 SumCharge = ArrayofSumCharge[m]->GetMatrixArray() ;
4008 for ( Int_t i = i_one ; i < ROWS-1 ; i+=i_one ) {
4009 Int_t i__ = i*COLUMNS;
4010 for ( Int_t j = j_one ; j < COLUMNS-1 ; j+=j_one ) {
4011 Int_t i__j = i__ + j;
4013 ArrayV[i__j] = ( coef2[i] * ArrayV[i__j - one__]
4014 + tempRatioZ * ( ArrayV[i__j - j_one] + ArrayV[i__j + j_one] )
4015 + coef1[i] * ArrayV[i__j + one__]
4016 + coef3[i] * ( ArrayVP[i__j] + ArrayVM[i__j] )
4018 ) * OverRelaxcoef4[i]
4019 - OverRelaxM1*ArrayV[i__j] ;
4024 if ( k == ITERATIONS && ( iWillDivide || jWillDivide ) ) {
4025 for ( Int_t i = i_one ; i < ROWS-1 ; i+=i_one ) {
4026 Int_t i__ = i*COLUMNS;
4027 for ( Int_t j = j_one ; j < COLUMNS-1 ; j+=j_one ) {
4028 Int_t i__j = i__ + j;
4030 if ( iWillDivide ) {
4031 ArrayV[i__j + half__] = ( ArrayV[i__j + one__] + ArrayV[i__j] ) / 2 ;
4033 ArrayV[i__j - half__] = ( ArrayV[j] + ArrayV[one__ + j] ) / 2 ;
4036 if ( jWillDivide ) {
4037 ArrayV[i__j + half] = ( ArrayV[i__j + j_one] + ArrayV[i__j] ) / 2 ;
4039 ArrayV[i__j - half] = ( ArrayV[i__] + ArrayV[i__ + j_one] ) / 2 ;
4041 if ( iWillDivide ) {
4042 ArrayV[i__j + half__ + half] = ( ArrayV[i__j + one__ + j_one] + ArrayV[i__j]
4043 + ArrayV[i__j + one__ ] + ArrayV[i__j + j_one] ) / 4 ;
4045 ArrayV[i__j - half__ - half] = ( ArrayV[j - j_one] + ArrayV[i__j]
4046 + ArrayV[j] + ArrayV[i__j - j_one] ) / 4 ;
4047 else if ( j == j_one )
4048 ArrayV[i__j - half__ - half] = ( ArrayV[i__ - one__] + ArrayV[i__j]
4049 + ArrayV[i__] + ArrayV[i__j - one__] ) / 4 ;
4062 if (iWillDivide) { i_one = i_one / 2 ;
if ( i_one < 1 ) i_one = 1 ; }
4063 if (jWillDivide) { j_one = j_one / 2 ;
if ( j_one < 1 ) j_one = 1 ; }
4068 Float_t coef10,coef11,coef12 ;
4069 coef10 = -0.5 / GRIDSIZER ;
4070 coef11 = -0.5 / GRIDSIZEPHI ;
4071 coef12 = (GRIDSIZEZ/3.0) / (-1*StarMagE) ;
4073 Int_t one__ = COLUMNS;
4078 for ( Int_t m = 0 ; m < PHISLICES ; m++ )
4080 ArrayV = ArrayofArrayV[m]->GetMatrixArray();
4081 EroverEz = ArrayofEroverEz[m]->GetMatrixArray();
4083 for ( Int_t j = COLUMNS-1 ; j >= 0 ; j-- )
4086 for ( Int_t i = 1 ; i < ROWS-1 ; i++ ) {
4087 Int_t i__j = i*COLUMNS + j;
4088 ArrayE[i__j] = coef10 * ( ArrayV[i__j + one__] - ArrayV[i__j - one__] ) ;
4090 Int_t r__j = (ROWS-1)*one__ + j;
4091 ArrayE[j] = coef10 * ( - ArrayV[2*one__ + j] + 4*ArrayV[one__ + j] - 3*ArrayV[j] ) ;
4092 ArrayE[r__j] = coef10 * ( 3*ArrayV[r__j] - 4*ArrayV[r__j - one__] + ArrayV[r__j - 2*one__] ) ;
4094 for ( Int_t i = 0 ; i < ROWS ; i++ )
4097 Int_t i__ = i*COLUMNS;
4098 Int_t i__j = i__ + j;
4099 Int_t i__c = i__ + COLUMNS-1;
4100 EroverEz[i__j] = 0.0;
4101 if ( j == COLUMNS-2 ) EroverEz[i__j] = coef12 * 1.5 * (ArrayE[i__c - 1] + ArrayE[i__c]) ;
4102 else if ( j != COLUMNS-1 )
4104 for ( Int_t k = i__j ; k <= i__c ; k++ )
4106 EroverEz[i__j] += Index*ArrayE[k] ;
4107 if ( Index != 4 ) Index = 4;
else Index = 2 ;
4109 if ( Index == 4 ) EroverEz[i__j] -= ArrayE[i__c] ;
4110 else if ( Index == 2 ) EroverEz[i__j] += (0.5*ArrayE[i__c - 1] - 2.5*ArrayE[i__c]) ;
4111 EroverEz[i__j] *= coef12 ;
4122 for ( Int_t m = 0 ; m < PHISLICES ; m++ )
4127 if ( m_plus > PHISLICES-1 ) m_plus = PHISLICES - 2 ;
4128 if ( m_minus < 0 ) m_minus = 1 ;
4131 if ( m_plus > PHISLICES-1 ) m_plus = 0 ;
4132 if ( m_minus < 0 ) m_minus = PHISLICES - 1 ;
4134 ArrayVP = ArrayofArrayV[m_plus]->GetMatrixArray() ;
4135 ArrayVM = ArrayofArrayV[m_minus]->GetMatrixArray() ;
4136 EPhioverEz = ArrayofEPhioverEz[m]->GetMatrixArray() ;
4137 for ( Int_t j = COLUMNS-1 ; j >= 0 ; j-- )
4140 for ( Int_t i = 0 ; i < ROWS ; i++ )
4142 Int_t i__j = i*COLUMNS + j;
4143 Float_t Radius = IFCRadius + i*GRIDSIZER ;
4144 ArrayE[i__j] = coef11 * ( ArrayVP[i__j] - ArrayVM[i__j] ) / Radius ;
4147 for ( Int_t i = 0 ; i < ROWS ; i++ )
4149 Int_t i__ = i*COLUMNS;
4150 Int_t i__j = i__ + j;
4151 Int_t i__c = i__ + COLUMNS-1;
4153 EPhioverEz[i__j] = 0.0 ;
4154 if ( j == COLUMNS-2 ) EPhioverEz[i__j] = coef12 * 1.5 * (ArrayE[i__c - 1] + ArrayE[i__c]) ;
4155 else if ( j != COLUMNS-1 )
4157 for ( Int_t k = i__j ; k <= i__c ; k++ )
4159 EPhioverEz[i__j] += Index*ArrayE[k] ;
4160 if ( Index != 4 ) Index = 4;
else Index = 2 ;
4162 if ( Index == 4 ) EPhioverEz[i__j] -= ArrayE[i__c] ;
4163 else if ( Index == 2 ) EPhioverEz[i__j] += (0.5*ArrayE[i__c - 1] - 2.5*ArrayE[i__c]) ;
4164 EPhioverEz[i__j] *= coef12 ;
4173 for ( Int_t m = 0 ; m < PHISLICES ; m++ )
4175 ArrayofSumCharge[m] -> Delete() ;
4215 const Prime PrimaryOrGlobal, Float_t x_new[3], Float_t p_new[3],
4216 const unsigned int RowMask1 ,
const unsigned int RowMask2 ,
const Float_t VertexError)
4219 x_new[0] = x[0] ; x_new[1] = x[1] ; x_new[2] = x[2] ;
4220 p_new[0] = p[0] ; p_new[1] = p[1] ; p_new[2] = p[2] ;
4222 SectorNumber( sec, x ) ;
4225 if ( finite((
double)Charge)*finite(x[0])*finite(x[1])*finite(x[2])*finite(p[0])*finite(p[1])*finite(p[2]) == 0 ) return ;
4227 const Int_t ROWS = TPCROWS[sec-1] ;
4228 const Float_t TestRadius = 77.00 ;
4231 Float_t B[3], Rotation, Direction, xx[3], xxprime[3] ;
4232 Double_t Xtrack[ROWS], Ytrack[ROWS], Ztrack[ROWS] ;
4233 Double_t Xtrack1[ROWS], Ytrack1[ROWS], Ztrack1[ROWS] ;
4234 Double_t R[ROWS], dX[ROWS], dY[ROWS], C0, X0, Y0, R0, Pt, R2, theta, theta0, DeltaTheta ;
4235 Double_t Xprime[ROWS+1], Yprime[ROWS+1], eX[ROWS+1], eY[ROWS+1] ;
4238 ChargeB = Charge * TMath::Sign((
int)1,(
int)(B[2]*1000)) ;
4239 Pt = TMath::Sqrt( p[0]*p[0] + p[1]*p[1] ) ;
4240 R0 = TMath::Abs( 1000.0 * Pt / ( 0.299792 * B[2] ) ) ;
4241 X0 = x[0] + ChargeB * p[1] * R0 / Pt ;
4242 Y0 = x[1] - ChargeB * p[0] * R0 / Pt ;
4243 Rotation = TMath::Sign( (
double)1.0, (x[0]-X0)*p[1] - (x[1]-Y0)*p[0] ) ;
4245 memcpy(R,TPCROWR[sec-1],ROWS*
sizeof(Double_t));
4248 if (Y0 == 0.0) Direction = TMath::Sign((
float)1.0,p[1]) ;
4252 R2 = TestRadius * TestRadius ;
4253 C0 = ( R2 - R0*R0 + X0*X0 + Y0*Y0 ) ;
4254 Double_t X1 = 0.5 * ( C0*X0 - TMath::Sqrt( TMath::Abs( C0*C0*X0*X0 - (C0*C0 - 4*Y0*Y0*R2) * (X0*X0+Y0*Y0) ) ) ) / (X0*X0 + Y0*Y0) ;
4255 Double_t Y1 = ( R2 - R0*R0 - 2*X0*X1 + X0*X0 + Y0*Y0 ) / ( 2 * Y0 ) ;
4256 Double_t X2 = 0.5 * ( C0*X0 + TMath::Sqrt( TMath::Abs( C0*C0*X0*X0 - (C0*C0 - 4*Y0*Y0*R2) * (X0*X0+Y0*Y0) ) ) ) / (X0*X0 + Y0*Y0) ;
4257 Double_t Y2 = ( R2 - R0*R0 - 2*X0*X2 + X0*X0 + Y0*Y0 ) / ( 2 * Y0 ) ;
4258 if ( X2*p[0] + Y2*p[1] < X1*p[0] + Y1*p[1] ) Direction = -1.0 ;
4261 theta0 = TMath::ATan2( (x[1]-Y0) , (x[0]-X0) ) ;
4262 Xprime[0] = theta0 ;
4265 eY[0] = VertexError ;
4268 unsigned int OneBit = 1 ;
4269 for ( Int_t i = 0 ; i < ROWS ; i++ )
4271 if ( ( i < 24 ) && ( ( RowMask1 & OneBit<<(i+8) ) == 0 ) ) continue ;
4272 if ( ( i >= 24 ) && ( ( RowMask2 & OneBit<<(i-24) ) == 0 ) ) continue ;
4274 C0 = ( R[i]*R[i] - R0*R0 + X0*X0 + Y0*Y0 ) ;
4275 if (Y0 == 0.0) Xtrack[index] = 0.5 * C0 / X0 ;
4276 else Xtrack[index] = 0.5*( C0*X0 + Direction*TMath::Sqrt( TMath::Abs( C0*C0*X0*X0 - (C0*C0-4*Y0*Y0*R[i]*R[i])*(X0*X0+Y0*Y0) )) )
4278 if (Y0 == 0.0) Ytrack[index] = Direction * TMath::Sqrt( TMath::Abs( R[i]*R[i] - Xtrack[index]*Xtrack[index] ) ) ;
4279 else Ytrack[index] = ( R[i]*R[i] - R0*R0 - 2*X0*Xtrack[index] + X0*X0 + Y0*Y0 ) / ( 2 * Y0 ) ;
4280 DeltaTheta = TMath::ATan2(x[1]-Y0,x[0]-X0) - TMath::ATan2(Ytrack[index]-Y0,Xtrack[index]-X0) ;
4281 while ( DeltaTheta < -1*TMath::Pi() ) DeltaTheta += TMath::TwoPi() ;
4282 while ( DeltaTheta >= TMath::Pi() ) DeltaTheta -= TMath::TwoPi() ;
4283 Ztrack[index] = x[2] - Rotation*DeltaTheta*R0*p[2] / Pt ;
4284 xx[0] = Xtrack[index] ; xx[1] = Ytrack[index] ; xx[2] = Ztrack[index] ;
4286 xx[0] = Xtrack[index] - (xxprime[0]-xx[0]) ; xx[1] = Ytrack[index] - (xxprime[1]-xx[1]) ; xx[2] = Ztrack[index] - (xxprime[2]-xx[2]) ;
4288 Xtrack1[index] = xxprime[0] ; Ytrack1[index] = xxprime[1] ; Ztrack1[index] = xxprime[2] ;
4289 theta = TMath::ATan2( (Ytrack[index]-Y0) , (Xtrack[index]-X0) ) ;
4290 while ( (theta - theta0) < -1*TMath::Pi() ) theta = theta + 2*TMath::Pi() ;
4291 while ( (theta - theta0) >= TMath::Pi() ) theta = theta - 2*TMath::Pi() ;
4292 dX[index] = Xtrack1[index] - Xtrack[index] ;
4293 dY[index] = Ytrack1[index] - Ytrack[index] ;
4294 Xprime[index+1] = theta ;
4295 Yprime[index+1] = dY[index]*TMath::Sin(theta) + dX[index]*TMath::Cos(theta) ;
4296 eX[index+1] = 0.5 / R0 ;
4297 eY[index+1] = 0.0100 ;
4299 if ( index == -1 ) return ;
4301 TGraphErrors gre(index-PrimaryOrGlobal+2,&Xprime[PrimaryOrGlobal],&Yprime[PrimaryOrGlobal],&eX[PrimaryOrGlobal],&eY[PrimaryOrGlobal]) ;
4302 TF1 FIT(
"myFIT",
"[0] + [1]*sin(x) + [2]*cos(x)" );
4303 FIT.SetParameter( 0, 0. );
4304 FIT.SetParameter( 1, 0. );
4305 FIT.SetParameter( 2, 0. );
4306 gre.Fit(
"myFIT",
"NQ") ;
4326 Double_t X0_new = X0 + FIT.GetParameter( 2 ) ;
4327 Double_t Y0_new = Y0 + FIT.GetParameter( 1 ) ;
4328 Double_t R0_new = R0 + FIT.GetParameter( 0 ) ;
4329 Double_t Pt_new = TMath::Abs( R0_new * 0.299792 * B[2] / 1000. ) ;
4331 if ( TMath::Sqrt( x[0]*x[0]+x[1]*x[1] ) <= IFCRadius )
4332 { x_new[0] = x[0] ; x_new[1] = x[1] ; x_new[2] = x[2] ; }
4336 xx[0] = x[0] - (xxprime[0]-x[0]) ; xx[1] = x[1] - (xxprime[1]-x[1]) ; xx[2] = x[2] - (xxprime[2]-x[2]) ;
4340 Int_t count = 0 ; p_new[2] = 0.0 ;
4341 for ( Int_t i = 0 ; i < index+1 ; i++ )
4343 DeltaTheta = (TMath::ATan2(Ytrack1[i]-Y0_new,Xtrack1[i]-X0_new)-TMath::ATan2(x_new[1]-Y0_new,x_new[0]-X0_new)) ;
4344 while ( DeltaTheta < -1*TMath::Pi() ) DeltaTheta += TMath::TwoPi() ;
4345 while ( DeltaTheta >= TMath::Pi() ) DeltaTheta -= TMath::TwoPi() ;
4346 if ( DeltaTheta != 0 ) { p_new[2] += (Ztrack1[i]-x_new[2]) / DeltaTheta ; count += 1 ; }
4349 p_new[0] = Pt_new * ( x_new[1] - Y0_new ) / ( ChargeB * R0_new ) ;
4350 p_new[1] = Pt_new * ( X0_new - x_new[0] ) / ( ChargeB * R0_new ) ;
4351 p_new[2] *= Pt_new / ( Rotation * R0_new * count ) ;
4391 const Prime PrimaryOrGlobal, Int_t &new_Charge, Float_t x_new[3], Float_t p_new[3],
4392 const unsigned int RowMask1,
const unsigned int RowMask2,
const Float_t VertexError )
4395 x_new[0] = x[0] ; x_new[1] = x[1] ; x_new[2] = x[2] ;
4396 p_new[0] = p[0] ; p_new[1] = p[1] ; p_new[2] = p[2] ;
4398 SectorNumber(sec, x);
4400 if ( finite((
double)Charge)*finite(x[0])*finite(x[1])*finite(x[2])*finite(p[0])*finite(p[1])*finite(p[2]) == 0 ) return ;
4402 const Float_t InnerOuterRatio = 0.6 ;
4403 const Int_t ROWS = TPCROWS[sec-1] ;
4404 const Int_t RefIndex = 7 ;
4405 const Int_t MinHits = 15 ;
4406 const Int_t DEBUG = 0 ;
4409 Float_t B[3], Direction, xx[3], xxprime[3] ;
4410 Double_t Xreference, Yreference ;
4411 Double_t Xtrack[ROWS], Ytrack[ROWS], Ztrack[ROWS] ;
4412 Double_t R[ROWS], C0, X0, Y0, R0, Pt, R2, DeltaTheta ;
4413 Double_t Xprime[ROWS+1], Yprime[ROWS+1], Zprime[ROWS+1], dX[ROWS+1], dY[ROWS+1] ;
4414 Double_t U[ROWS+1], V[ROWS+1], eU[ROWS+1], eV[ROWS+1] ;
4419 Double_t tempSpaceChargeR2 = SpaceChargeR2 ;
4420 ManualSpaceChargeR2(sc,SpaceChargeEWRatio);
4423 ChargeB = Charge * TMath::Sign((
int)1,(
int)(B[2]*1000)) ;
4424 Pt = TMath::Sqrt( p[0]*p[0] + p[1]*p[1] ) ;
4425 R0 = TMath::Abs( 1000.0 * Pt / ( 0.299792 * B[2] ) ) ;
4426 X0 = x[0] + ChargeB * p[1] * R0 / Pt ;
4427 Y0 = x[1] - ChargeB * p[0] * R0 / Pt ;
4429 memcpy(R,TPCROWR[sec-1],ROWS*
sizeof(Double_t));
4433 if (TMath::Abs(Y0) < 0.001 ) Direction = TMath::Sign( (
float)1.0, p[1] ) ;
4437 R2 = R[RefIndex]*R[RefIndex] ;
4438 C0 = ( R2 - R0*R0 + X0*X0 + Y0*Y0 ) ;
4439 Double_t X1 = 0.5 * ( C0*X0 - TMath::Sqrt( TMath::Abs( C0*C0*X0*X0 - (C0*C0 - 4*Y0*Y0*R2) * (X0*X0+Y0*Y0) ) ) )
4441 Double_t Y1 = ( R2 - R0*R0 - 2*X0*X1 + X0*X0 + Y0*Y0 ) / ( 2 * Y0 ) ;
4442 Double_t X2 = 0.5 * ( C0*X0 + TMath::Sqrt( TMath::Abs( C0*C0*X0*X0 - (C0*C0 - 4*Y0*Y0*R2) * (X0*X0+Y0*Y0) ) ) )
4444 Double_t Y2 = ( R2 - R0*R0 - 2*X0*X2 + X0*X0 + Y0*Y0 ) / ( 2 * Y0 ) ;
4445 if ( X2*p[0] + Y2*p[1] < X1*p[0] + Y1*p[1] ) Direction = -1.0 ;
4448 Xreference = Yreference = 0.0 ;
4449 for ( Int_t i = 0 ; i < ROWS ; i++ )
4451 C0 = ( R[i]*R[i] - R0*R0 + X0*X0 + Y0*Y0 ) ;
4452 if ( TMath::Abs(Y0) < 0.001 ) Xtrack[i] = 0.5 * C0 / X0 ;
4453 else Xtrack[i] = 0.5*( C0*X0 + Direction*TMath::Sqrt( TMath::Abs( C0*C0*X0*X0 -
4454 (C0*C0-4*Y0*Y0*R[i]*R[i])*(X0*X0+Y0*Y0) )) ) / (X0*X0+Y0*Y0) ;
4455 if ( TMath::Abs(Y0) < 0.001 ) Ytrack[i] = Direction * TMath::Sqrt( TMath::Abs( R[i]*R[i] - Xtrack[i]*Xtrack[i] ) ) ;
4456 else Ytrack[i] = ( R[i]*R[i] - R0*R0 - 2*X0*Xtrack[i] + X0*X0 + Y0*Y0 ) / ( 2 * Y0 ) ;
4457 DeltaTheta = TMath::ATan2(x[1]-Y0,x[0]-X0) - TMath::ATan2(Ytrack[i]-Y0,Xtrack[i]-X0) ;
4458 while ( DeltaTheta < -1*TMath::Pi() ) DeltaTheta += TMath::TwoPi() ;
4459 while ( DeltaTheta >= TMath::Pi() ) DeltaTheta -= TMath::TwoPi() ;
4460 Ztrack[i] = x[2] + ChargeB*DeltaTheta*R0*p[2] / Pt ;
4461 xx[0] = Xtrack[i] ; xx[1] = Ytrack[i] ; xx[2] = Ztrack[i] ;
4465 Xtrack[i] = xxprime[0] ; Ytrack[i] = xxprime[1] ; Ztrack[i] = xxprime[2] ;
4467 if ( i == RefIndex )
4469 Xreference = Xtrack[i] ;
4470 Yreference = Ytrack[i] ;
4475 unsigned int OneBit = 1 ;
4476 for ( Int_t i = 0 ; i < ROWS ; i++ )
4478 if ( ( i < 24 ) && (( RowMask1 & OneBit<<(i+8) ) == 0 )) continue ;
4479 if ( ( i >= 24 ) && (( RowMask2 & OneBit<<(i-24) ) == 0 )) continue ;
4481 Xprime[Index] = Xtrack[i] ; Yprime[Index] = Ytrack[i] ; Zprime[Index] = Ztrack[i] ;
4482 dX[Index] = 0.2 ; dY[Index] = 1.0 ;
4484 if ( i < INNER[sec-1] ) { dX[Index] *= InnerOuterRatio ; dY[Index] *= InnerOuterRatio ; } ;
4488 Xprime[0] = x[0] ; Yprime[0] = x[1] ; Zprime[0] = x[2] ; dX[0] = VertexError ; dY[0] = VertexError ;
4492 for ( Int_t i = PrimaryOrGlobal ; i < Index+1 ; i++ )
4494 Double_t zero = 0.001 ;
4495 Double_t displacement2 ;
4498 (Xprime[i]-Xreference)*(Xprime[i]-Xreference) + (Yprime[i]-Yreference)*(Yprime[i]-Yreference) ;
4500 if ( displacement2 > zero )
4503 U[count] = ( Xprime[i] - Xreference ) / displacement2 ;
4504 V[count] = ( Yprime[i] - Yreference ) / displacement2 ;
4505 eU[count] = dX[i] / displacement2 ;
4506 eV[count] = dY[i] / displacement2 ;
4510 if ( count < MinHits ) return ;
4512 TGraphErrors gre( count+1, U, V, eU, eV ) ;
4513 gre.Fit(
"pol1",
"Q") ;
4514 TF1 *fit = gre.GetFunction(
"pol1" ) ;
4518 TCanvas* c1 =
new TCanvas(
"A Simple Fit",
"The Fit", 250, 10, 700, 500) ;
4522 TCanvas* c2 =
new TCanvas(
"The circles are OK",
"The circles ", 250, 800, 700, 500) ;
4523 TGraph* gra =
new TGraph( Index-PrimaryOrGlobal+1, &Xprime[PrimaryOrGlobal], &Yprime[PrimaryOrGlobal] ) ;
4525 gra -> SetMaximum(200) ;
4526 gra -> SetMinimum(-200) ;
4528 gra -> GetXaxis() -> SetLimits(-200.,200.) ;
4533 Double_t X0_new = Xreference - ( fit->GetParameter(1) / ( 2.0 * fit->GetParameter(0) ) ) ;
4534 Double_t Y0_new = Yreference + ( 1.0 / ( 2.0 * fit->GetParameter(0) ) ) ;
4535 Double_t R0_new = TMath::Sqrt( (Xreference-X0_new)*(Xreference-X0_new) + (Yreference-Y0_new)*(Yreference-Y0_new) ) ;
4536 Double_t Pt_new = TMath::Abs ( R0_new * 0.299792 * B[2] / 1000. ) ;
4542 if ( TMath::Sqrt( x[0]*x[0]+x[1]*x[1] ) <= IFCRadius )
4543 { x_new[0] = x[0] ; x_new[1] = x[1] ; x_new[2] = x[2] ; }
4549 Int_t icount = 0 ; p_new[2] = 0.0 ;
4550 for ( Int_t i = 0 ; i < Index+1 ; i++ )
4552 DeltaTheta = (TMath::ATan2(Yprime[i]-Y0_new,Xprime[i]-X0_new)-TMath::ATan2(x_new[1]-Y0_new,x_new[0]-X0_new)) ;
4553 while ( DeltaTheta < -1*TMath::Pi() ) DeltaTheta += TMath::TwoPi() ;
4554 while ( DeltaTheta >= TMath::Pi() ) DeltaTheta -= TMath::TwoPi() ;
4555 if ( DeltaTheta != 0 ) { p_new[2] += (Zprime[i]-x_new[2]) / DeltaTheta ; icount += 1 ; }
4558 p_new[0] = Pt_new * ( x_new[1] - Y0_new ) / ( ChargeB * R0_new ) ;
4559 p_new[1] = Pt_new * ( X0_new - x_new[0] ) / ( ChargeB * R0_new ) ;
4560 p_new[2] *= Pt_new / ( -1 * ChargeB * R0_new * icount ) ;
4563 Float_t change = TMath::Abs( TMath::ATan2(Y0,X0) - TMath::ATan2(Y0_new,X0_new) ) ;
4564 if ( change > 0.9*TMath::Pi() && change < 1.1*TMath::Pi() ) new_Charge = -1 * Charge ;
4565 else new_Charge = Charge ;
4568 fSpaceChargeR2 = tempfSpaceChargeR2 ;
4569 SpaceChargeR2 = tempSpaceChargeR2 ;
4580 Float_t DCA,
const unsigned int RowMask1,
const unsigned int RowMask2, Float_t &pSpace )
4582 const Int_t ROWS = TPCROWS[sec-1] ;
4583 const Int_t RefIndex = 7 ;
4584 const Int_t MinInnerTPCHits = 5 ;
4585 const Int_t MinOuterTPCHits = 10 ;
4586 const Int_t DEBUG = 0 ;
4588 unsigned int OneBit = 1 ;
4589 Int_t InnerTPCHits = 0, OuterTPCHits = 0 ;
4590 for ( Int_t i = 0 ; i < ROWS ; i++ )
4592 if ( i < INNER[sec-1] )
4594 if ( ( i < 24 ) && (( RowMask1 & OneBit<<(i+8) ) != 0 )) InnerTPCHits++ ;
4595 if ( ( i >= 24 ) && (( RowMask2 & OneBit<<(i-24) ) != 0 )) InnerTPCHits++ ;
4599 if ( ( i < 24 ) && (( RowMask1 & OneBit<<(i+8) ) != 0 )) OuterTPCHits++ ;
4600 if ( ( i >= 24 ) && (( RowMask2 & OneBit<<(i-24) ) != 0 )) OuterTPCHits++ ;
4606 if ( (Pt < 0.3) || (Pt > 2.0) )
return(1) ;
4607 if ( (VertexZ < -50) || (VertexZ > 50) )
return(2) ;
4608 if ( (PseudoRapidity < -1.0) || (PseudoRapidity > 1.0) )
return(3) ;
4609 if ( (Charge != 1) && (Charge != -1) )
return(4) ;
4610 if ( (DCA < -4.0) || (DCA > 4.0) )
return(5) ;
4611 if ( InnerTPCHits < MinInnerTPCHits )
return(6) ;
4612 if ( OuterTPCHits < MinOuterTPCHits )
return(7) ;
4615 Float_t B[3], Direction, xx[3], xxprime[3] ;
4616 Double_t Xreference, Yreference ;
4617 Double_t Xtrack[ROWS], Ytrack[ROWS], Ztrack[ROWS] ;
4618 Double_t R[ROWS], C0, X0, Y0, R0, Pz_over_Pt, Z_coef, DeltaTheta ;
4619 Double_t Xprime[ROWS+1], Yprime[ROWS+1], dX[ROWS+1], dY[ROWS+1] ;
4620 Double_t U[ROWS+1], V[ROWS+1], eU[ROWS+1], eV[ROWS+1] ;
4624 Double_t tempSpaceChargeR2 = SpaceChargeR2 ;
4625 if (!useManualSCForPredict) ManualSpaceChargeR2(0.01,SpaceChargeEWRatio);
4627 Float_t x[3] = { 0, 0, 0 } ;
4629 ChargeB = Charge * TMath::Sign((
int)1,(
int)(B[2]*1000)) ;
4630 R0 = TMath::Abs( 1000.0 * Pt / ( 0.299792 * B[2] ) ) ;
4631 X0 = ChargeB * 0.707107 * R0 ;
4632 Y0 = ChargeB * -0.707107 * R0 ;
4633 Pz_over_Pt = TMath::SinH(PseudoRapidity) ;
4634 Z_coef = ChargeB*R0*Pz_over_Pt ;
4636 memcpy(R,TPCROWR[sec-1],ROWS*
sizeof(Double_t));
4639 Float_t InnerOuterRatio = 0.0 ;
4640 Xreference = Yreference = 0.0 ;
4642 for ( Int_t i = 0 ; i < ROWS ; i++ )
4644 C0 = ( R[i]*R[i] - R0*R0 + X0*X0 + Y0*Y0 ) ;
4645 if ( TMath::Abs(Y0) < 0.001 ) Xtrack[i] = 0.5 * C0 / X0 ;
4646 else Xtrack[i] = 0.5*( C0*X0 + Direction*TMath::Sqrt( TMath::Abs( C0*C0*X0*X0 -
4647 (C0*C0-4*Y0*Y0*R[i]*R[i])*(X0*X0+Y0*Y0) )) ) / (X0*X0+Y0*Y0) ;
4648 if ( TMath::Abs(Y0) < 0.001 ) Ytrack[i] = Direction * TMath::Sqrt( TMath::Abs( R[i]*R[i] - Xtrack[i]*Xtrack[i] ) ) ;
4649 else Ytrack[i] = ( R[i]*R[i] - R0*R0 - 2*X0*Xtrack[i] + X0*X0 + Y0*Y0 ) / ( 2 * Y0 ) ;
4650 DeltaTheta = TMath::ATan2(-1*Y0,-1*X0) - TMath::ATan2(Ytrack[i]-Y0,Xtrack[i]-X0) ;
4651 while ( DeltaTheta < -1*TMath::Pi() ) DeltaTheta += TMath::TwoPi() ;
4652 while ( DeltaTheta >= TMath::Pi() ) DeltaTheta -= TMath::TwoPi() ;
4653 Ztrack[i] = VertexZ + DeltaTheta*Z_coef ;
4654 xx[0] = Xtrack[i] ; xx[1] = Ytrack[i] ; xx[2] = Ztrack[i] ;
4656 if (mDistortionMode & (kSpaceCharge | kSpaceChargeR2 | kSpaceChargeFXT)) {
4658 InnerOuterRatio = 1.3 ;
4659 for (
unsigned int j = 0 ; j < 3; ++j )
4664 if (mDistortionMode & (kGridLeak | k3DGridLeak | kFullGridLeak)) {
4666 InnerOuterRatio = 0.6 ;
4667 for (
unsigned int j = 0 ; j < 3 ; ++j )
4673 Xtrack[i] = 2*Xtrack[i] - xx[0] ;
4674 Ytrack[i] = 2*Ytrack[i] - xx[1] ;
4675 Ztrack[i] = 2*Ztrack[i] - xx[2] ;
4677 if ( i == RefIndex )
4679 Xreference = Xtrack[i] ;
4680 Yreference = Ytrack[i] ;
4685 for ( Int_t i = 0 ; i < ROWS ; i++ )
4687 if ( ( i < 24 ) && (( RowMask1 & OneBit<<(i+8) ) == 0 )) continue ;
4688 if ( ( i >= 24 ) && (( RowMask2 & OneBit<<(i-24) ) == 0 )) continue ;
4690 Xprime[Index] = Xtrack[i] ; Yprime[Index] = Ytrack[i] ;
4691 dX[Index] = 0.2 ; dY[Index] = 1.0 ;
4693 if ( i < INNER[sec-1] ) { dX[Index] *= InnerOuterRatio ; dY[Index] *= InnerOuterRatio ; } ;
4698 for ( Int_t i = 0 ; i < Index+1 ; i++ )
4700 Double_t zero = 0.001 ;
4701 Double_t displacement2 ;
4704 (Xprime[i]-Xreference)*(Xprime[i]-Xreference) + (Yprime[i]-Yreference)*(Yprime[i]-Yreference) ;
4706 if ( displacement2 > zero )
4709 U[count] = ( Xprime[i] - Xreference ) / displacement2 ;
4710 V[count] = ( Yprime[i] - Yreference ) / displacement2 ;
4711 eU[count] = dX[i] / displacement2 ;
4712 eV[count] = dY[i] / displacement2 ;
4716 TGraphErrors gre( count+1, U, V, eU, eV ) ;
4717 gre.Fit(
"pol1",
"Q") ;
4718 TF1 *fit = gre.GetFunction(
"pol1" ) ;
4722 TCanvas* c1 =
new TCanvas(
"A Simple Fit",
"The Fit", 250, 10, 700, 500) ;
4726 TCanvas* c2 =
new TCanvas(
"The circles are OK",
"The circles ", 250, 800, 700, 500) ;
4727 TGraph* gra =
new TGraph( Index+1, Xprime, Yprime ) ;
4729 gra -> SetMaximum(200) ;
4730 gra -> SetMinimum(-200) ;
4732 gra -> GetXaxis() -> SetLimits(-200.,200.) ;
4737 Double_t X0_new = Xreference - ( fit->GetParameter(1) / ( 2.0 * fit->GetParameter(0) ) ) ;
4738 Double_t Y0_new = Yreference + ( 1.0 / ( 2.0 * fit->GetParameter(0) ) ) ;
4739 Double_t R0_new = TMath::Sqrt( (Xreference-X0_new)*(Xreference-X0_new) + (Yreference-Y0_new)*(Yreference-Y0_new) ) ;
4740 Double_t DCA_new = TMath::Sqrt( X0_new*X0_new + Y0_new*Y0_new ) - R0_new ;
4742 pSpace = (DCA * ChargeB) * SpaceChargeR2 / DCA_new ;
4745 fSpaceChargeR2 = tempfSpaceChargeR2 ;
4746 SpaceChargeR2 = tempSpaceChargeR2 ;
4782 Float_t PseudoRapidity, Float_t Phi, Float_t DCA,
4783 const unsigned long long RowMask1,
4784 const unsigned long long RowMask2,
4785 Float_t RowMaskErrorR[64], Float_t RowMaskErrorRPhi[64], Float_t &pSpace )
4788 const Int_t INNERDETECTORS = 6 ;
4789 const Int_t SSDLAYERS = 1 ;
4790 const Int_t MinInnerTPCHits = 5 ;
4791 const Int_t MinOuterTPCHits = 10 ;
4792 const Int_t DEBUG = 0 ;
4794 const Int_t TPCOFFSET = INNERDETECTORS + SSDLAYERS + 1 ;
4795 const Int_t BITS = INNERDETECTORS + TPCROWS[sec-1] + SSDLAYERS + 1 ;
4797 unsigned int OneBit = 1U ;
4798 Int_t InnerTPCHits = 0, OuterTPCHits = 0 ;
4799 for ( Int_t i = 0 ; i < TPCROWS[sec-1] ; i++ )
4801 if ( i < INNER[sec-1] )
4803 if ( ( i < 24 ) && (( RowMask1 & OneBit<<(i+8) ) != 0 )) InnerTPCHits++ ;
4804 if ( ( i >= 24 ) && (( RowMask2 & OneBit<<(i-24) ) != 0 )) InnerTPCHits++ ;
4808 if ( ( i < 24 ) && (( RowMask1 & OneBit<<(i+8) ) != 0 )) OuterTPCHits++ ;
4809 if ( ( i >= 24 ) && (( RowMask2 & OneBit<<(i-24) ) != 0 )) OuterTPCHits++ ;
4815 if ( (Pt < 0.3) || (Pt > 2.0) )
return(1) ;
4816 if ( (VertexZ < -50) || (VertexZ > 50) )
return(2) ;
4817 if ( (PseudoRapidity < -1.0) || (PseudoRapidity > 1.0) )
return(3) ;
4818 if ( (Charge != 1) && (Charge != -1) )
return(4) ;
4819 if ( (DCA < -4.0) || (DCA > 4.0) )
return(5) ;
4820 if ( InnerTPCHits < MinInnerTPCHits )
return(6) ;
4821 if ( OuterTPCHits < MinOuterTPCHits )
return(7) ;
4823 Int_t ChargeB, HitSector ;
4824 Float_t B[3], xx[3], xxprime[3] ;
4825 Double_t Xtrack[BITS], Ytrack[BITS], Ztrack[BITS] ;
4826 Double_t R[BITS], X0, Y0, X0Prime, Y0Prime, R0, Pz_over_Pt, Z_coef, DeltaTheta ;
4827 Double_t Xprime[BITS+1], Yprime[BITS+1], dX[BITS+1], dY[BITS+1] ;
4828 Float_t PhiPrime, HitPhi, HitLocalPhi ;
4829 Double_t cosPhi, sinPhi, cosPhiPrime, sinPhiPrime, cosPhiMPrime, sinPhiMPrime ;
4833 Double_t tempSpaceChargeR2 = SpaceChargeR2 ;
4834 if (!useManualSCForPredict) ManualSpaceChargeR2(0.01,SpaceChargeEWRatio);
4837 xx[0] = TPCROWR[2][0]; xx[1] = 0; xx[2] = 50;
4840 Int_t tempDistortionMode = mDistortionMode;
4841 mDistortionMode = (tempDistortionMode & (kSpaceChargeR2 | kSpaceChargeFXT));
4842 if (tempDistortionMode & kFullGridLeak) mDistortionMode |= kFullGridLeak ;
4843 else if (tempDistortionMode & k3DGridLeak ) mDistortionMode |= k3DGridLeak ;
4844 else if (tempDistortionMode & kGridLeak ) mDistortionMode |= kGridLeak ;
4846 Float_t x[3] = { 0, 0, 0 } ;
4848 ChargeB = Charge * TMath::Sign((
int)1,(
int)(B[2]*1000)) ;
4849 R0 = TMath::Abs( 1000.0 * Pt / ( 0.299792 * B[2] ) ) ;
4850 X0 = ChargeB * 0.0 * R0 ;
4851 Y0 = ChargeB * -1.0 * R0 ;
4852 Pz_over_Pt = TMath::SinH(PseudoRapidity) ;
4853 Z_coef = ChargeB*R0*Pz_over_Pt ;
4856 while ( PhiPrime < 0 ) PhiPrime += PiOver6 ;
4857 PhiPrime = fmod(PhiPrime + PiOver12,PiOver6) - PiOver12;
4858 cosPhi = TMath::Cos(Phi);
4859 sinPhi = TMath::Sin(Phi);
4860 cosPhiPrime = TMath::Cos(PhiPrime);
4861 sinPhiPrime = TMath::Sin(PhiPrime);
4862 cosPhiMPrime = cosPhi*cosPhiPrime+sinPhi*sinPhiPrime;
4863 sinPhiMPrime = sinPhi*cosPhiPrime-cosPhi*sinPhiPrime;
4865 X0Prime = cosPhiPrime*X0 - sinPhiPrime*Y0;
4866 Y0Prime = sinPhiPrime*X0 + cosPhiPrime*Y0;
4883 memcpy(&(R[TPCOFFSET]),TPCROWR[sec-1],TPCROWS[sec-1]*
sizeof(Double_t));
4886 for ( Int_t i = 0 ; i < BITS ; i++ )
4889 if ( ( i > 0 && i < 32 ) && (( RowMask1 & OneBit<<(i) ) == 0 )) continue ;
4890 if ( ( i > 0 && i >= 32 ) && (( RowMask2 & OneBit<<(i-32) ) == 0 )) continue ;
4892 if ( R[i] < IFCRadius || R[i] > OFCRadius ) {
4893 Ytrack[i] = -1 * ChargeB * ( R[i]*R[i]/(2*R0) ) ;
4894 Xtrack[i] = TMath::Sqrt( R[i]*R[i] - Ytrack[i]*Ytrack[i] ) ;
4895 DeltaTheta = TMath::ATan2(-1*Y0,-1*X0) - TMath::ATan2(Ytrack[i]-Y0,Xtrack[i]-X0) ;
4896 while ( DeltaTheta < -1*TMath::Pi() ) DeltaTheta += TMath::TwoPi() ;
4897 while ( DeltaTheta >= TMath::Pi() ) DeltaTheta -= TMath::TwoPi() ;
4898 Ztrack[i] = VertexZ + DeltaTheta*Z_coef;
4905 Ytrack[i] = ChargeB * TMath::Sqrt( R0*R0 - (Xtrack[i]-X0Prime)*(Xtrack[i]-X0Prime) ) + Y0Prime ;
4906 HitLocalPhi = TMath::ATan2(Ytrack[i],Xtrack[i]);
4907 while (TMath::Abs(HitLocalPhi) > PiOver12) {
4909 PhiPrime -= TMath::Sign(PiOver6,HitLocalPhi);
4910 cosPhiPrime = TMath::Cos(PhiPrime);
4911 sinPhiPrime = TMath::Sin(PhiPrime);
4912 cosPhiMPrime = cosPhi*cosPhiPrime+sinPhi*sinPhiPrime;
4913 sinPhiMPrime = sinPhi*cosPhiPrime-cosPhi*sinPhiPrime;
4915 X0Prime = cosPhiPrime*X0 - sinPhiPrime*Y0;
4916 Y0Prime = sinPhiPrime*X0 + cosPhiPrime*Y0;
4918 Ytrack[i] = ChargeB * TMath::Sqrt( R0*R0 - (Xtrack[i]-X0Prime)*(Xtrack[i]-X0Prime) ) + Y0Prime ;
4919 HitLocalPhi = TMath::ATan2(Ytrack[i],Xtrack[i]);
4922 DeltaTheta = TMath::ATan2(-1*Y0Prime,-1*X0Prime) - TMath::ATan2(Ytrack[i]-Y0Prime,Xtrack[i]-X0Prime) ;
4923 while ( DeltaTheta < -1*TMath::Pi() ) DeltaTheta += TMath::TwoPi() ;
4924 while ( DeltaTheta >= TMath::Pi() ) DeltaTheta -= TMath::TwoPi() ;
4925 Ztrack[i] = VertexZ + DeltaTheta*Z_coef;
4928 xx[0] = cosPhiMPrime*Xtrack[i] - sinPhiMPrime*Ytrack[i];
4929 xx[1] = sinPhiMPrime*Xtrack[i] + cosPhiMPrime*Ytrack[i];
4932 if (mDistortionMode & (kGridLeak | k3DGridLeak | kFullGridLeak)) {
4933 HitPhi = TMath::ATan2(xx[1],xx[0]) ;
4934 while ( HitPhi < 0 ) HitPhi += TMath::TwoPi() ;
4935 while ( HitPhi >= TMath::TwoPi() ) HitPhi -= TMath::TwoPi() ;
4937 SectorNumber ( HitSector, HitPhi, xx[2] );
4938 if ( GLWeights[HitSector] < 0) {
4940 fSpaceChargeR2 = tempfSpaceChargeR2 ;
4941 SpaceChargeR2 = tempSpaceChargeR2 ;
4942 mDistortionMode = tempDistortionMode ;
4948 Xtrack[i] = cosPhi*xxprime[0] + sinPhi*xxprime[1] ;
4949 Ytrack[i] = -sinPhi*xxprime[0] + cosPhi*xxprime[1] ;
4950 Ztrack[i] = xxprime[2] ;
4955 Int_t TPCIndex = -1 ;
4956 for ( Int_t i = 1 ; i < BITS ; i++ )
4958 if ( ( i < 32 ) && (( RowMask1 & OneBit<<(i) ) == 0 )) continue ;
4959 if ( ( i >= 32 ) && (( RowMask2 & OneBit<<(i-32) ) == 0 )) continue ;
4960 Index++ ;
if ( i >= TPCOFFSET ) TPCIndex++ ;
4961 Xprime[Index] = Xtrack[i] ; Yprime[Index] = Ytrack[i] ;
4962 dX[Index] = RowMaskErrorR[i] ; dY[Index] = RowMaskErrorRPhi[i] ;
4966 TGraphErrors gre(Index+1,Xprime,Yprime,dX,dY) ;
4970 TF1 newDCA (
"newDCA" ,
"( [1] + [3] * sqrt( [1]**2 - x**2 + 2*x*[0] + [2]*[2] - [2]*(2*sqrt([0]**2+[1]**2))) )" );
4971 newDCA.SetParameters( X0, Y0, 0.0, ChargeB );
4972 newDCA.FixParameter(3, ChargeB);
4973 gre.Fit(
"newDCA",
"Q") ;
4974 Double_t DCA_new = newDCA.GetParameter( 2 ) ;
4979 TCanvas* c1 =
new TCanvas(
"A Simple Fit",
"The Fit", 250, 10, 700, 500) ;
4980 TCanvas* c2 =
new TCanvas(
"The circles are OK",
"The circles ", 250, 800, 700, 500) ;
4984 TGraph* gra =
new TGraph( Index+1, Xprime, Yprime ) ;
4986 gra -> SetMaximum(200) ;
4987 gra -> SetMinimum(-200) ;
4989 gra -> GetXaxis() -> SetLimits(-200.,200.) ;
4994 pSpace = (DCA * ChargeB) * SpaceChargeR2 / DCA_new ;
4997 fSpaceChargeR2 = tempfSpaceChargeR2 ;
4998 SpaceChargeR2 = tempSpaceChargeR2 ;
4999 mDistortionMode = tempDistortionMode ;
5028 Float_t PseudoRapidity, Float_t Phi, Float_t DCA,
5030 Double_t ErrorR[128],
5031 Double_t ErrorRPhi[128], Float_t &pSpace )
5034 const Int_t MinInnerTPCHits = 5 ;
5035 const Int_t MinOuterTPCHits = 10 ;
5036 const Int_t DEBUG = 0 ;
5040 if ( (Pt < 0.3) || (Pt > 2.0) )
return(1) ;
5041 if ( (VertexZ < -50) || (VertexZ > 50) )
return(2) ;
5042 if ( (PseudoRapidity < -1.0) || (PseudoRapidity > 1.0) )
return(3) ;
5043 if ( (Charge != 1) && (Charge != -1) )
return(4) ;
5044 if ( (DCA < -4.0) || (DCA > 4.0) )
return(5) ;
5046 Int_t InnerTPCHits = 0, OuterTPCHits = 0 ;
5047 for ( Int_t i = 0 ; i < NHits ; i++ )
5050 if ( R[i] > IFCRadius && R[i] < 1.035*GAPRADIUS ) InnerTPCHits++ ;
5051 if ( R[i] < OFCRadius && R[i] > 0.990*GAPRADIUS ) OuterTPCHits++ ;
5054 if ( InnerTPCHits < MinInnerTPCHits )
return(6) ;
5055 if ( OuterTPCHits < MinOuterTPCHits )
return(7) ;
5056 InnerTPCHits = 0 ; OuterTPCHits = 0 ;
5058 Int_t ChargeB, HitSector ;
5059 Float_t B[3], xx[3], xxprime[3] ;
5060 Double_t Xtrack[NHits], Ytrack[NHits], Ztrack[NHits] ;
5061 Double_t X0, Y0, R0, Pz_over_Pt, Z_coef, DeltaTheta ;
5063 Double_t cosPhi, sinPhi ;
5067 Double_t tempSpaceChargeR2 = SpaceChargeR2 ;
5068 if (!useManualSCForPredict) ManualSpaceChargeR2(0.01,SpaceChargeEWRatio);
5070 if (DoOnce) UndoDistortion(0,0,0);
5071 Int_t tempDistortionMode = mDistortionMode;
5072 mDistortionMode = (tempDistortionMode & (kSpaceCharge | kSpaceChargeR2 | kSpaceChargeFXT | kGridLeak | k3DGridLeak | kFullGridLeak));
5074 memset(xx,0,3*
sizeof(Float_t));
5076 ChargeB = Charge * TMath::Sign((
int)1,(
int)(B[2]*1000)) ;
5077 R0 = TMath::Abs( 1000.0 * Pt / ( 0.299792 * B[2] ) ) ;
5078 X0 = ChargeB * 0.0 * R0 ;
5079 Y0 = ChargeB * -1.0 * R0 ;
5080 Pz_over_Pt = TMath::SinH(PseudoRapidity) ;
5081 Z_coef = ChargeB*R0*Pz_over_Pt ;
5083 cosPhi = TMath::Cos(Phi);
5084 sinPhi = TMath::Sin(Phi);
5086 for ( Int_t i = 0 ; i < NHits ; i++ )
5089 Ytrack[i] = -1 * ChargeB * ( R[i]*R[i]/(2*R0) ) ;
5090 Xtrack[i] = TMath::Sqrt( R[i]*R[i] - Ytrack[i]*Ytrack[i] ) ;
5091 DeltaTheta = TMath::ATan2(-1*Y0,-1*X0) - TMath::ATan2(Ytrack[i]-Y0,Xtrack[i]-X0) ;
5092 while ( DeltaTheta < -1*TMath::Pi() ) DeltaTheta += TMath::TwoPi() ;
5093 while ( DeltaTheta >= TMath::Pi() ) DeltaTheta -= TMath::TwoPi() ;
5094 Ztrack[i] = VertexZ + DeltaTheta*Z_coef;
5096 if ( R[i] < IFCRadius || R[i] > OFCRadius ||
5097 TMath::Abs(Ztrack[i]) > TPC_Z0-0.5) continue ;
5107 xx[0] = cosPhi*Xtrack[i] - sinPhi*Ytrack[i];
5108 xx[1] = sinPhi*Xtrack[i] + cosPhi*Ytrack[i];
5111 HitPhi = TMath::ATan2(xx[1],xx[0]) ;
5112 while ( HitPhi < 0 ) HitPhi += TMath::TwoPi() ;
5113 while ( HitPhi >= TMath::TwoPi() ) HitPhi -= TMath::TwoPi() ;
5115 SectorNumber ( HitSector, HitPhi, xx[2] );
5117 if (mDistortionMode & (kGridLeak | k3DGridLeak | kFullGridLeak)) {
5118 if ( GLWeights[HitSector] < 0) {
5120 fSpaceChargeR2 = tempfSpaceChargeR2 ;
5121 SpaceChargeR2 = tempSpaceChargeR2 ;
5122 mDistortionMode = tempDistortionMode ;
5126 HitPhi = fmod(HitPhi + PiOver12, PiOver6) - PiOver12;
5127 if ( R[i]*TMath::Cos(HitPhi) < GAPRADIUS ) InnerTPCHits++ ;
5128 else OuterTPCHits++ ;
5131 Xtrack[i] = cosPhi*xxprime[0] + sinPhi*xxprime[1] ;
5132 Ytrack[i] = -sinPhi*xxprime[0] + cosPhi*xxprime[1] ;
5133 Ztrack[i] = xxprime[2] ;
5137 if ( InnerTPCHits < MinInnerTPCHits ) {
5139 fSpaceChargeR2 = tempfSpaceChargeR2 ;
5140 SpaceChargeR2 = tempSpaceChargeR2 ;
5141 mDistortionMode = tempDistortionMode ;
5144 if ( OuterTPCHits < MinOuterTPCHits ) {
5146 fSpaceChargeR2 = tempfSpaceChargeR2 ;
5147 SpaceChargeR2 = tempSpaceChargeR2 ;
5148 mDistortionMode = tempDistortionMode ;
5152 TGraphErrors gre(NHits,Xtrack,Ytrack,ErrorR,ErrorRPhi) ;
5156 static TF1* newDCA = 0;
5157 if (!newDCA) newDCA =
new TF1(
"newDCA" ,
"( [1] + [3] * sqrt( [1]**2 - x**2 + 2*x*[0] + [2]*[2] - [2]*(2*sqrt([0]**2+[1]**2))) )" );
5158 newDCA->SetParameters( X0, Y0, 0.0, ChargeB );
5159 newDCA->FixParameter(3, ChargeB);
5160 gre.Fit(newDCA,
"Q") ;
5161 Double_t DCA_new = newDCA->GetParameter( 2 ) ;
5166 TCanvas* c1 =
new TCanvas(
"A Simple Fit",
"The Fit", 250, 10, 700, 500) ;
5167 TCanvas* c2 =
new TCanvas(
"The circles are OK",
"The circles ", 250, 800, 700, 500) ;
5171 TGraph* gra =
new TGraph( NHits, Xtrack, Ytrack ) ;
5173 gra -> SetMaximum(200) ;
5174 gra -> SetMinimum(-200) ;
5176 gra -> GetXaxis() -> SetLimits(-200.,200.) ;
5181 pSpace = (DCA * ChargeB) * SpaceChargeR2 / DCA_new ;
5184 fSpaceChargeR2 = tempfSpaceChargeR2 ;
5185 SpaceChargeR2 = tempSpaceChargeR2 ;
5186 mDistortionMode = tempDistortionMode ;
5197 Int_t StMagUtilities::IsPowerOfTwo(Int_t i)
5200 while( i > 0 ) { j += (i&1) ; i = (i>>1) ; }
5201 if ( j == 1 )
return(1) ;
5215 void StMagUtilities::SectorNumber( Int_t& Sector ,
const Float_t x[] )
5217 if ( Sector > 0 ) return ;
5218 Float_t phi = (usingCartesian ? TMath::ATan2(x[1],x[0]) : x[1]) ;
5219 SectorNumber( Sector, phi, x[2] ) ;
5221 void StMagUtilities::SectorNumber( Int_t& Sector , Float_t phi,
const Float_t z )
5223 if ( Sector > 0 ) return ;
5224 if ( phi < 0 ) phi += TMath::TwoPi() ;
5225 Sector = ( ( 30 - (int)(phi/PiOver12) )%24 ) / 2 ;
5226 if ( z < 0 ) Sector = 24 - Sector ;
5227 else if ( Sector == 0 ) Sector = 12 ;
5231 Int_t StMagUtilities::SectorSide( Int_t& Sector ,
const Float_t x[] )
5233 return SectorSide(Sector, x[2]) ;
5235 Int_t StMagUtilities::SectorSide( Int_t& Sector ,
const Float_t z )
5237 if ( Sector <= 0 ) SectorNumber(Sector,0,z);
5238 return (Sector <= 12 ? 1 : -1) ;
5250 Float_t StMagUtilities::LimitZ( Int_t& Sector ,
const Float_t x[] )
5253 const Float_t zlimit = 0.2;
5254 Int_t sign = SectorSide(Sector, z);
5255 if (sign * z < zlimit) z = sign * zlimit;
5269 int active_options = mDistortionMode & (kGridLeak | k3DGridLeak | kFullGridLeak);
5270 if (active_options & (active_options-1)) {
5271 cout <<
"StMagUtilities ERROR **** Do not use multiple GridLeak modes at the same time" << endl ;
5272 cout <<
"StMagUtilities ERROR **** These routines have overlapping functionality." << endl ;
5276 if (mDistortionMode & kGridLeak) {
5278 }
else if (mDistortionMode & k3DGridLeak) {
5280 }
else if (mDistortionMode & kFullGridLeak) {
5299 const Int_t ORDER = 1 ;
5300 const Int_t ROWS = 513 ;
5301 const Int_t COLUMNS = 129 ;
5302 const Int_t ITERATIONS = 750 ;
5303 const Double_t GRIDSIZER = (OFCRadius-IFCRadius) / (ROWS-1) ;
5304 const Double_t GRIDSIZEZ = TPC_Z0 / (COLUMNS-1) ;
5307 static TMatrix EroverEz(ROWS,COLUMNS) ;
5308 static Float_t Rlist[ROWS], Zedlist[COLUMNS] ;
5310 Float_t Er_integral, Ephi_integral ;
5311 Double_t r, phi, z ;
5313 if ( InnerGridLeakStrength == 0 && MiddlGridLeakStrength == 0 && OuterGridLeakStrength == 0 )
5314 { memcpy(Xprime,x,threeFloats) ; return ; }
5318 cout <<
"StMagUtilities::UndoGridL Please wait for the tables to fill ... ~30 seconds" << endl ;
5319 TMatrix ArrayV(ROWS,COLUMNS), Charge(ROWS,COLUMNS) ;
5322 for ( Int_t j = 0 ; j < COLUMNS ; j++ )
5324 Double_t zed = j*GRIDSIZEZ ;
5326 for ( Int_t i = 0 ; i < ROWS ; i++ )
5328 Double_t Radius = IFCRadius + i*GRIDSIZER ;
5343 Float_t GainRatio = 3.0 ;
5344 for ( Int_t j = 1 ; j < COLUMNS-1 ; j++ )
5346 for ( Int_t i = 1 ; i < ROWS-1 ; i++ )
5348 Double_t Radius = IFCRadius + i*GRIDSIZER ;
5350 if ( (Radius >= InnerGridLeakRadius) && (Radius < MiddlGridLeakRadius) )
5351 Charge(i,j) += InnerGridLeakStrength / (1.0e-3*Radius*Radius) ;
5353 if ( (Radius >= MiddlGridLeakRadius-MiddlGridLeakWidth) && (Radius <= MiddlGridLeakRadius+MiddlGridLeakWidth) )
5354 Charge(i,j) += MiddlGridLeakStrength ;
5356 if ( (Radius >= MiddlGridLeakRadius) && (Radius < OuterGridLeakRadius) )
5357 Charge(i,j) += OuterGridLeakStrength / (GainRatio*1.0e-3*Radius*Radius) ;
5361 PoissonRelaxation( ArrayV, Charge, EroverEz, ITERATIONS ) ;
5365 if (usingCartesian) Cart2Polar(x,r,phi);
5366 else { r = x[0]; phi = x[1]; }
5367 if ( phi < 0 ) phi += TMath::TwoPi() ;
5368 z = LimitZ( Sector, x ) ;
5370 Float_t phi0 = PiOver6 * ((Int_t)(0.499 + phi/PiOver6)) ;
5371 Float_t local_y = r * TMath::Cos( phi - phi0 ) ;
5374 Er_integral = Interpolate2DTable( ORDER, local_y, TMath::Abs(z), ROWS, COLUMNS, Rlist, Zedlist, EroverEz ) ;
5375 Ephi_integral = 0.0 ;
5378 if (fSpaceChargeR2) GetSpaceChargeR2();
5383 Float_t Weight = SpaceChargeR2 * (doingDistortion ? SmearCoefSC*SmearCoefGL : 1.0);
5384 if ( z < 0) Weight *= SpaceChargeEWRatio ;
5385 phi = phi - Weight * ( Const_0*Ephi_integral - Const_1*Er_integral ) / r ;
5386 r = r - Weight * ( Const_0*Er_integral + Const_1*Ephi_integral ) ;
5389 if (usingCartesian) Polar2Cart(r,phi,Xprime);
5390 else { Xprime[0] = r; Xprime[1] = phi; }
5407 const Int_t ORDER = 1 ;
5408 const Int_t neR3D = 73 ;
5409 const Int_t ITERATIONS = 100 ;
5410 const Int_t SYMMETRY = 1 ;
5411 const Int_t ROWS = 513 ;
5412 const Int_t COLUMNS = 129 ;
5413 const Int_t PHISLICES = 8 ;
5414 const Float_t GRIDSIZEPHI = (2 - SYMMETRY) * TMath::Pi() / (12.0*(PHISLICES-1)) ;
5415 const Float_t GRIDSIZER = (OFCRadius-IFCRadius) / (ROWS-1) ;
5416 const Float_t GRIDSIZEZ = TPC_Z0 / (COLUMNS-1) ;
5418 static TMatrix *ArrayofArrayV[PHISLICES], *ArrayofCharge[PHISLICES] ;
5419 static TMatrix *ArrayofEroverEz[PHISLICES], *ArrayofEPhioverEz[PHISLICES] ;
5421 static Float_t Rlist[ROWS], Zedlist[COLUMNS] ;
5423 static TMatrix *ArrayoftiltEr[PHISLICES] ;
5424 static TMatrix *ArrayoftiltEphi[PHISLICES] ;
5427 static Float_t eRadius[neR3D] = { 50.0, 60.0, 70.0, 80.0, 90.0,
5428 100.0, 104.0, 106.5, 109.0, 111.5,
5429 114.0, 115.0, 116.0, 117.0, 118.0,
5430 118.5, 118.75, 119.0, 119.25, 119.5,
5431 119.75, 120.0, 120.25, 120.5, 120.75,
5432 121.0, 121.25, 121.5, 121.75, 122.0,
5433 122.25, 122.5, 122.75, 123.0, 123.25,
5434 123.5, 123.75, 124.0, 124.25, 124.5,
5435 124.75, 125.0, 125.25, 125.5, 126.0,
5436 126.25, 126.5, 126.75, 127.0, 127.25,
5437 127.5, 128.0, 128.5, 129.0, 129.5,
5438 130.0, 130.5, 131.0, 131.5, 132.0,
5439 133.0, 135.0, 137.5, 140.0, 145.0,
5440 150.0, 160.0, 170.0, 180.0, 190.0,
5441 195.0, 198.0, 200.0 } ;
5443 static Float_t Philist[PHISLICES] ;
5445 for ( Int_t k = 0 ; k < PHISLICES ; k++ ) Philist[k] = GRIDSIZEPHI * k ;
5447 Float_t Er_integral, Ephi_integral ;
5450 memcpy(Xprime,x,threeFloats) ;
5452 if ( MiddlGridLeakStrength == 0 ) return ;
5456 cout <<
"StMagUtilities::Undo3DGrid Please wait for the tables to fill ... ~5 seconds * PHISLICES" << endl ;
5458 for ( Int_t k = 0 ; k < PHISLICES ; k++ )
5460 ArrayoftiltEr[k] =
new TMatrix(neR3D,EMap_nZ) ;
5461 ArrayoftiltEphi[k] =
new TMatrix(neR3D,EMap_nZ) ;
5464 for ( Int_t k = 0 ; k < PHISLICES ; k++ )
5466 ArrayofArrayV[k] =
new TMatrix(ROWS,COLUMNS) ;
5467 ArrayofCharge[k] =
new TMatrix(ROWS,COLUMNS) ;
5468 ArrayofEroverEz[k] =
new TMatrix(ROWS,COLUMNS) ;
5469 ArrayofEPhioverEz[k] =
new TMatrix(ROWS,COLUMNS) ;
5472 for ( Int_t k = 0 ; k < PHISLICES ; k++ )
5474 TMatrix &ArrayV = *ArrayofArrayV[k] ;
5475 TMatrix &Charge = *ArrayofCharge[k] ;
5476 Float_t cosPhiK = TMath::Cos(Philist[k]) ;
5479 for ( Int_t i = 0 ; i < ROWS ; i++ )
5481 Rlist[i] = IFCRadius + i*GRIDSIZER ;
5482 for ( Int_t j = 0 ; j < COLUMNS ; j++ )
5484 Zedlist[j] = j * GRIDSIZEZ ;
5490 for ( Int_t i = 1 ; i < ROWS-1 ; i++ )
5492 Float_t Radius = IFCRadius + i*GRIDSIZER ;
5493 Float_t local_y_hi = (Radius+GRIDSIZER/2.0) * cosPhiK ;
5494 Float_t local_y_lo = (Radius-GRIDSIZER/2.0) * cosPhiK ;
5495 Float_t charge_y_hi = OUTERGGFirst ;
5496 Float_t charge_y_lo = INNERGGLast ;
5497 for ( Int_t j = 1 ; j < COLUMNS-1 ; j++ )
5500 Float_t top = 0, bottom = 0 ;
5502 if (local_y_hi > charge_y_lo && local_y_hi < charge_y_hi)
5505 if (local_y_lo > charge_y_lo && local_y_lo < charge_y_hi)
5506 bottom = local_y_lo ;
5508 bottom = charge_y_lo ;
5511 if (local_y_lo > charge_y_lo && local_y_lo < charge_y_hi)
5513 bottom = local_y_lo ;
5514 if (local_y_hi > charge_y_lo && local_y_hi < charge_y_hi)
5520 Float_t Weight = 1. / (local_y_hi*local_y_hi - local_y_lo*local_y_lo) ;
5523 const Float_t BackwardsCompatibilityRatio = 3.75 ;
5524 Charge(i,j) = (top*top-bottom*bottom) * Weight * MiddlGridLeakStrength*BackwardsCompatibilityRatio ;
5533 Poisson3DRelaxation( ArrayofArrayV, ArrayofCharge, ArrayofEroverEz, ArrayofEPhioverEz, PHISLICES,
5534 GRIDSIZEPHI, ITERATIONS, SYMMETRY) ;
5538 for ( Int_t k = 0 ; k < PHISLICES ; k++ )
5540 TMatrix &tiltEr = *ArrayoftiltEr[k] ;
5541 TMatrix &tiltEphi = *ArrayoftiltEphi[k] ;
5543 for ( Int_t i = 0 ; i < EMap_nZ ; i++ )
5545 z = TMath::Abs(eZList[i]) ;
5546 for ( Int_t j = 0 ; j < neR3D ; j++ )
5549 tiltEr(j,i) = Interpolate3DTable( ORDER,r,z,phi,ROWS,COLUMNS,PHISLICES,Rlist,Zedlist,Philist,ArrayofEroverEz ) ;
5550 tiltEphi(j,i) = Interpolate3DTable( ORDER,r,z,phi,ROWS,COLUMNS,PHISLICES,Rlist,Zedlist,Philist,ArrayofEPhioverEz ) ;
5555 for ( Int_t k = 0 ; k < PHISLICES ; k++ )
5557 ArrayofArrayV[k] -> Delete() ;
5558 ArrayofCharge[k] -> Delete() ;
5559 ArrayofEroverEz[k] -> Delete() ;
5560 ArrayofEPhioverEz[k] -> Delete() ;
5565 if (usingCartesian) Cart2Polar(x,r,phi);
5566 else { r = x[0]; phi = x[1]; }
5567 if ( phi < 0 ) phi += TMath::TwoPi() ;
5568 z = LimitZ( Sector, x ) ;
5570 Float_t phi_prime, local_y, r_eff, FLIP = 1.0 ;
5572 if ( SYMMETRY == 1 )
5574 Int_t N = (int)(phi/PiOver12) ;
5575 phi_prime = phi - N * PiOver12 ;
5576 if ( TMath::Power(-1,N) < 0 ) phi_prime = PiOver12 - phi_prime ;
5577 if ( TMath::Power(-1,N) < 0 ) FLIP = -1 ;
5581 local_y = r * TMath::Cos(phi_prime) ;
5582 if ( local_y > GAPRADIUS - WIREGAP && local_y < GAPRADIUS ) r_eff = (GAPRADIUS - WIREGAP) / TMath::Cos(phi_prime) ;
5583 if ( local_y < GAPRADIUS + WIREGAP && local_y > GAPRADIUS ) r_eff = (GAPRADIUS + WIREGAP) / TMath::Cos(phi_prime) ;
5586 Er_integral = Interpolate3DTable( ORDER, r_eff, TMath::Abs(z), phi_prime, neR3D, EMap_nZ, PHISLICES, eRadius, eZList, Philist, ArrayoftiltEr ) ;
5587 Ephi_integral = Interpolate3DTable( ORDER, r_eff, TMath::Abs(z), phi_prime, neR3D, EMap_nZ, PHISLICES, eRadius, eZList, Philist, ArrayoftiltEphi ) ;
5588 Ephi_integral *= FLIP ;
5590 if (fSpaceChargeR2) GetSpaceChargeR2();
5596 Float_t Weight = SpaceChargeR2 * (doingDistortion ? SmearCoefSC*SmearCoefGL : 1.0);
5597 if (GLWeights[Sector] >= 0) Weight *= GLWeights[Sector] ;
5598 if ( z < 0) Weight *= SpaceChargeEWRatio ;
5599 phi = phi - Weight * ( Const_0*Ephi_integral - Const_1*Er_integral ) / r ;
5600 r = r - Weight * ( Const_0*Er_integral + Const_1*Ephi_integral ) ;
5603 if (usingCartesian) Polar2Cart(r,phi,Xprime);
5604 else { Xprime[0] = r; Xprime[1] = phi; }
5621 const Int_t ORDER = 1 ;
5622 const Int_t neR3D = 132 ;
5623 const Int_t ITERATIONS = 80 ;
5624 const Int_t ROWS = 513 ;
5625 const Int_t COLUMNS = 65 ;
5626 const Int_t PHISLICES = 180 ;
5628 const Int_t PHISLICES1 = PHISLICES+1 ;
5629 const Float_t GRIDSIZEPHI = TMath::TwoPi() / PHISLICES ;
5630 const Float_t GRIDSIZER = (OFCRadius-IFCRadius) / (ROWS-1) ;
5631 const Float_t GRIDSIZEZ = TPC_Z0 / (COLUMNS-1) ;
5633 static TMatrix *ArrayofArrayV[PHISLICES] , *ArrayofCharge[PHISLICES] ;
5634 static TMatrix *ArrayofEroverEzW[PHISLICES], *ArrayofEPhioverEzW[PHISLICES] ;
5635 static TMatrix *ArrayofEroverEzE[PHISLICES], *ArrayofEPhioverEzE[PHISLICES] ;
5637 static Float_t Rlist[ROWS], Zedlist[COLUMNS] ;
5639 static TMatrix *ArrayoftiltEr[PHISLICES1] ;
5640 static TMatrix *ArrayoftiltEphi[PHISLICES1] ;
5643 static Float_t eRadius[neR3D] = { 50.0, 50.5, 51.0, 51.25, 51.5,
5644 51.75, 52.0, 52.25, 52.5, 52.75,
5645 53.0, 53.25, 53.5, 53.75, 54.0,
5646 54.25, 54.75, 55.0, 55.25, 55.5,
5647 56.0, 56.5, 57.0, 58.0, 60.0,
5648 62.0, 66.0, 70.0, 80.0, 90.0,
5649 100.0, 104.0, 106.5, 109.0, 111.5,
5650 114.0, 115.0, 116.0, 117.0, 118.0,
5651 118.5, 118.75, 119.0, 119.25, 119.5,
5652 119.75, 120.0, 120.25, 120.5, 120.75,
5653 121.0, 121.25, 121.5, 121.75, 122.0,
5654 122.25, 122.5, 122.75, 123.0, 123.25,
5655 123.5, 123.75, 124.0, 124.25, 124.5,
5656 124.75, 125.0, 125.25, 125.5, 125.75,
5657 126.0, 126.25, 126.5, 126.75, 127.0,
5658 127.25, 127.5, 128.0, 128.5, 129.0,
5659 129.5, 130.0, 130.5, 131.0, 131.5,
5660 132.0, 133.0, 135.0, 137.5, 140.0,
5661 145.0, 150.0, 160.0, 170.0, 180.0,
5662 184.0, 186.5, 189.0, 190.0, 190.5,
5663 190.75, 191.0, 191.25, 191.5, 191.75,
5664 192.0, 192.25, 192.5, 192.75, 193.0,
5665 193.25, 193.5, 193.75, 194.0, 194.25,
5666 194.5, 194.75, 195.0, 195.25, 195.5,
5667 196.0, 196.25, 196.5, 196.75, 197.0,
5668 197.25, 197.5, 198.0, 198.5, 199.0,
5671 static Float_t Philist [PHISLICES1] ;
5673 Float_t Er_integral, Ephi_integral ;
5676 memcpy(Xprime,x,threeFloats) ;
5678 if ( MiddlGridLeakStrength == 0 ) return ;
5682 cout <<
"StMagUtilities::UndoFullGrid Please wait for the tables to fill ... ~5 seconds * PHISLICES" << endl ;
5683 Int_t SeclistW[PHISLICES1] ;
5684 Int_t SeclistE[PHISLICES1] ;
5685 Float_t SecPhis [PHISLICES1] ;
5687 for ( Int_t i = 0 ; i < ROWS ; i++ ) Rlist[i] = IFCRadius + i*GRIDSIZER ;
5688 for ( Int_t j = 0 ; j < COLUMNS ; j++ ) Zedlist[j] = j * GRIDSIZEZ ;
5690 for ( Int_t k = 0 ; k < PHISLICES1 ; k++ )
5692 Philist[k] = GRIDSIZEPHI * k ;
5693 Int_t k2 = (12*k+PHISLICES/2)/PHISLICES;
5694 SeclistW[k] = (14-k2)%12+1;
5695 SeclistE[k] = (k2+8)%12+13;
5696 SecPhis[k] = GRIDSIZEPHI * (k2*PHISLICES/12-k) ;
5699 for ( Int_t k = 0 ; k < PHISLICES1 ; k++ )
5701 ArrayoftiltEr[k] =
new TMatrix(neR3D,EMap_nZ) ;
5702 ArrayoftiltEphi[k] =
new TMatrix(neR3D,EMap_nZ) ;
5705 for ( Int_t k = 0 ; k < PHISLICES ; k++ )
5707 ArrayofArrayV[k] =
new TMatrix(ROWS,COLUMNS) ;
5708 ArrayofCharge[k] =
new TMatrix(ROWS,COLUMNS) ;
5709 ArrayofEroverEzW[k] =
new TMatrix(ROWS,COLUMNS) ;
5710 ArrayofEPhioverEzW[k] =
new TMatrix(ROWS,COLUMNS) ;
5711 ArrayofEroverEzE[k] =
new TMatrix(ROWS,COLUMNS) ;
5712 ArrayofEPhioverEzE[k] =
new TMatrix(ROWS,COLUMNS) ;
5715 for ( Int_t m = -1; m < 2 ; m+=2 )
5717 TMatrix** ArrayofEroverEz = ( m < 0 ? ArrayofEroverEzW : ArrayofEroverEzE ) ;
5718 TMatrix** ArrayofEPhioverEz = ( m < 0 ? ArrayofEPhioverEzW : ArrayofEPhioverEzE ) ;
5719 Int_t* Seclist = ( m < 0 ? SeclistW : SeclistE ) ;
5721 for ( Int_t k = 0 ; k < PHISLICES ; k++ )
5723 TMatrix &ArrayV = *ArrayofArrayV[k] ;
5724 TMatrix &Charge = *ArrayofCharge[k] ;
5725 Float_t cosPhiK = TMath::Cos(SecPhis[k]) ;
5731 for ( Int_t i = 1 ; i < ROWS-1 ; i++ )
5733 Float_t Radius = IFCRadius + i*GRIDSIZER ;
5734 Float_t local_y_hi = (Radius+GRIDSIZER/2.0) * cosPhiK ;
5735 Float_t local_y_lo = (Radius-GRIDSIZER/2.0) * cosPhiK ;
5750 Float_t SCscale = SpaceChargeRadialDependence(Radius)/SpaceChargeRadialDependence(GAPRADIUS) ;
5752 Float_t top = 0, bottom = 0 ;
5753 Float_t local_charge = 0;
5755 for (Int_t l = 0; l < 4 ; l++ )
5758 if (local_y_hi < GL_charge_y_lo[l] || local_y_lo > GL_charge_y_hi[l])
continue;
5759 top = (local_y_hi > GL_charge_y_hi[l] ? GL_charge_y_hi[l] : local_y_hi);
5760 bottom = (local_y_lo < GL_charge_y_lo[l] ? GL_charge_y_lo[l] : local_y_lo);
5761 local_charge += SCscale * GLWeights[Seclist[k]-1 + l*24] * (top*top - bottom*bottom);
5765 Float_t Weight = 1.0 / (local_y_hi*local_y_hi - local_y_lo*local_y_lo) ;
5769 const Float_t BackwardsCompatibilityRatio = 3.75 ;
5771 for ( Int_t j = 1 ; j < COLUMNS-1 ; j++ )
5774 Charge(i,j) = local_charge * Weight * MiddlGridLeakStrength * BackwardsCompatibilityRatio ;
5783 Poisson3DRelaxation( ArrayofArrayV, ArrayofCharge, ArrayofEroverEz, ArrayofEPhioverEz, PHISLICES,
5784 GRIDSIZEPHI, ITERATIONS, 0) ;
5790 for ( Int_t k = 0 ; k < PHISLICES1 ; k++ )
5792 TMatrix &tiltEr = *ArrayoftiltEr[k] ;
5793 TMatrix &tiltEphi = *ArrayoftiltEphi[k] ;
5794 phi = Philist[k == PHISLICES ? 0 : k] ;
5795 for ( Int_t i = 0 ; i < EMap_nZ ; i++ )
5797 z = TMath::Abs(eZList[i]) ;
5798 TMatrix** ArrayofEroverEz = ( eZList[i] > 0 ? ArrayofEroverEzW : ArrayofEroverEzE ) ;
5799 TMatrix** ArrayofEPhioverEz = ( eZList[i] > 0 ? ArrayofEPhioverEzW : ArrayofEPhioverEzE ) ;
5800 for ( Int_t j = 0 ; j < neR3D ; j++ )
5803 tiltEr(j,i) = Interpolate3DTable( ORDER,r,z,phi,ROWS,COLUMNS,PHISLICES,Rlist,Zedlist,Philist,ArrayofEroverEz ) ;
5804 tiltEphi(j,i) = Interpolate3DTable( ORDER,r,z,phi,ROWS,COLUMNS,PHISLICES,Rlist,Zedlist,Philist,ArrayofEPhioverEz ) ;
5809 for ( Int_t k = 0 ; k < PHISLICES ; k++ )
5811 ArrayofArrayV[k] -> Delete() ;
5812 ArrayofCharge[k] -> Delete() ;
5813 ArrayofEroverEzW[k] -> Delete() ;
5814 ArrayofEPhioverEzW[k] -> Delete() ;
5815 ArrayofEroverEzE[k] -> Delete() ;
5816 ArrayofEPhioverEzE[k] -> Delete() ;
5821 if (usingCartesian) Cart2Polar(x,r,phi);
5822 else { r = x[0]; phi = x[1]; }
5823 if ( phi < 0 ) phi += TMath::TwoPi() ;
5824 z = LimitZ( Sector, x ) ;
5826 Float_t cos_phi_prime, local_y, r_eff ;
5827 cos_phi_prime = TMath::Cos(phi - PiOver6 * ((
int) ((phi/PiOver6)+0.5))) ;
5830 local_y = r * cos_phi_prime ;
5833 const Float_t margin = 0.7975 ;
5835 const Float_t middle_of_inner_sheet = 0.5 * (GL_charge_y_lo[0] + GL_charge_y_hi[0]);
5836 const Float_t middle_of_outer_sheet = 0.5 * (GL_charge_y_lo[3] + GL_charge_y_hi[3]);
5837 if ( local_y > GL_charge_y_lo[0] - margin && local_y < middle_of_inner_sheet) r_eff = (GL_charge_y_lo[0] - margin) / cos_phi_prime;
5838 else if ( local_y < GL_charge_y_hi[0] + margin && local_y > middle_of_inner_sheet) r_eff = (GL_charge_y_hi[0] + margin) / cos_phi_prime;
5839 else if ( local_y > GL_charge_y_lo[1] - margin && local_y < GL_charge_y_hi[1] ) r_eff = (GL_charge_y_lo[1] - margin) / cos_phi_prime;
5840 else if ( local_y < GL_charge_y_hi[2] + margin && local_y > GL_charge_y_lo[2] ) r_eff = (GL_charge_y_hi[2] + margin) / cos_phi_prime;
5841 else if ( local_y > GL_charge_y_lo[3] - margin && local_y < middle_of_outer_sheet) r_eff = (GL_charge_y_lo[3] - margin) / cos_phi_prime;
5842 else if ( local_y < GL_charge_y_hi[3] + margin && local_y > middle_of_outer_sheet) r_eff = (GL_charge_y_hi[3] + margin) / cos_phi_prime;
5844 Er_integral = Interpolate3DTable( ORDER, r_eff, z, phi, neR3D, EMap_nZ, PHISLICES1, eRadius, eZList, Philist, ArrayoftiltEr ) ;
5845 Ephi_integral = Interpolate3DTable( ORDER, r_eff, z, phi, neR3D, EMap_nZ, PHISLICES1, eRadius, eZList, Philist, ArrayoftiltEphi ) ;
5847 if (fSpaceChargeR2) GetSpaceChargeR2();
5853 Float_t Weight = SpaceChargeR2 * (doingDistortion ? SmearCoefSC*SmearCoefGL : 1.0);
5854 if ( z < 0 ) Weight *= SpaceChargeEWRatio ;
5855 phi = phi - Weight * ( Const_0*Ephi_integral - Const_1*Er_integral ) / r ;
5856 r = r - Weight * ( Const_0*Er_integral + Const_1*Ephi_integral ) ;
5859 if (usingCartesian) Polar2Cart(r,phi,Xprime);
5860 else { Xprime[0] = r; Xprime[1] = phi; }
5882 const Int_t ORDER = 1 ;
5883 const Int_t neR3D = 73 ;
5884 const Int_t ITERATIONS = 100 ;
5885 const Int_t ROWS = 257 ;
5886 const Int_t COLUMNS = 65 ;
5887 const Int_t PHISLICES = 180 ;
5889 const Int_t PHISLICES1 = PHISLICES+1 ;
5890 const Float_t GRIDSIZEPHI = TMath::TwoPi() / PHISLICES ;
5891 const Float_t GRIDSIZER = (OFCRadius-IFCRadius) / (ROWS-1) ;
5892 const Float_t GRIDSIZEZ = TPC_Z0 / (COLUMNS-1) ;
5894 const Float_t INNERGGSpan = INNERGGLast - INNERGGFirst ;
5895 const Float_t OUTERGGSpan = OUTERGGLast - OUTERGGFirst ;
5897 static TMatrix *ArrayofArrayV[PHISLICES] , *ArrayofCharge[PHISLICES] ;
5898 static TMatrix *ArrayofEroverEzW[PHISLICES], *ArrayofEPhioverEzW[PHISLICES] ;
5899 static TMatrix *ArrayofEroverEzE[PHISLICES], *ArrayofEPhioverEzE[PHISLICES] ;
5901 static Float_t Rlist[ROWS], Zedlist[COLUMNS] ;
5903 static TMatrix *ArrayoftiltEr[PHISLICES1] ;
5904 static TMatrix *ArrayoftiltEphi[PHISLICES1] ;
5907 static Float_t eRadius[neR3D] = { 50.0, 60.0, 70.0, 80.0, 90.0,
5908 100.0, 104.0, 106.5, 109.0, 111.5,
5909 114.0, 115.0, 116.0, 117.0, 118.0,
5910 118.5, 118.75, 119.0, 119.25, 119.5,
5911 119.75, 120.0, 120.25, 120.5, 120.75,
5912 121.0, 121.25, 121.5, 121.75, 122.0,
5913 122.25, 122.5, 122.75, 123.0, 123.25,
5914 123.5, 123.75, 124.0, 124.25, 124.5,
5915 124.75, 125.0, 125.25, 125.5, 126.0,
5916 126.25, 126.5, 126.75, 127.0, 127.25,
5917 127.5, 128.0, 128.5, 129.0, 129.5,
5918 130.0, 130.5, 131.0, 131.5, 132.0,
5919 133.0, 135.0, 137.5, 140.0, 145.0,
5920 150.0, 160.0, 170.0, 180.0, 190.0,
5921 195.0, 198.0, 200.0 } ;
5923 static Float_t Philist [PHISLICES1] ;
5925 Float_t Er_integral, Ephi_integral ;
5930 memcpy(Xprime,x,threeFloats) ;
5934 cout <<
"StMagUtilities::UndoSectorAlign Please wait for the tables to fill ... ~5 seconds * PHISLICES" << endl ;
5935 Int_t SeclistW[PHISLICES1] ;
5936 Int_t SeclistE[PHISLICES1] ;
5937 Float_t SecPhis [PHISLICES1] ;
5939 for ( Int_t i = 0 ; i < ROWS ; i++ ) Rlist[i] = IFCRadius + i*GRIDSIZER ;
5940 for ( Int_t j = 0 ; j < COLUMNS ; j++ ) Zedlist[j] = j * GRIDSIZEZ ;
5942 for ( Int_t k = 0 ; k < PHISLICES1 ; k++ )
5944 Philist[k] = GRIDSIZEPHI * k ;
5945 Int_t k2 = (12*k+PHISLICES/2)/PHISLICES;
5946 SeclistW[k] = (14-k2)%12+1;
5947 SeclistE[k] = (k2+8)%12+13;
5948 SecPhis[k] = GRIDSIZEPHI * (k2*PHISLICES/12-k) ;
5951 for ( Int_t k = 0 ; k < PHISLICES1 ; k++ )
5953 ArrayoftiltEr[k] =
new TMatrix(neR3D,EMap_nZ) ;
5954 ArrayoftiltEphi[k] =
new TMatrix(neR3D,EMap_nZ) ;
5957 for ( Int_t k = 0 ; k < PHISLICES ; k++ )
5959 ArrayofArrayV[k] =
new TMatrix(ROWS,COLUMNS) ;
5960 ArrayofCharge[k] =
new TMatrix(ROWS,COLUMNS) ;
5961 ArrayofEroverEzW[k] =
new TMatrix(ROWS,COLUMNS) ;
5962 ArrayofEPhioverEzW[k] =
new TMatrix(ROWS,COLUMNS) ;
5963 ArrayofEroverEzE[k] =
new TMatrix(ROWS,COLUMNS) ;
5964 ArrayofEPhioverEzE[k] =
new TMatrix(ROWS,COLUMNS) ;
5967 for ( Int_t m = -1; m < 2 ; m+=2 )
5969 TMatrix** ArrayofEroverEz = ( m < 0 ? ArrayofEroverEzW : ArrayofEroverEzE ) ;
5970 TMatrix** ArrayofEPhioverEz = ( m < 0 ? ArrayofEPhioverEzW : ArrayofEPhioverEzE ) ;
5971 Int_t* Seclist = ( m < 0 ? SeclistW : SeclistE ) ;
5974 for ( Int_t k = 0 ; k < PHISLICES ; k++ )
5976 TMatrix &ArrayV = *ArrayofArrayV[k] ;
5977 TMatrix &Charge = *ArrayofCharge[k] ;
5983 Double_t tanSecPhi = TMath::Tan(SecPhis[k]);
5984 Double_t cosSecPhi = TMath::Cos(SecPhis[k]);
5985 Double_t secSecPhi = 1.0/cosSecPhi;
5986 Double_t iOffsetFirst, iOffsetLast, oOffsetFirst, oOffsetLast;
5988 if (StTpcDb::instance())
5990 Double_t local[3] = {0,0,0};
5999 const TGeoHMatrix& iAlignMatrix = StTpcDb::instance()->SubSInner2Tpc(Seclist[k]);
6000 const TGeoHMatrix& oAlignMatrix = StTpcDb::instance()->SubSOuter2Tpc(Seclist[k]);
6014 Double_t *master = locP.position().xyz();
6020 local[0] = m * INNERGGFirst * tanSecPhi;
6021 local[1] = INNERGGFirst;
6023 iAlignMatrix.LocalToMaster(local,master);
6028 iOffsetFirst = (TPC_Z0 + m * master[2]) * StarMagE;
6030 local[0] = m * INNERGGLast * tanSecPhi;
6031 local[1] = INNERGGLast;
6033 iAlignMatrix.LocalToMaster(local,master);
6038 iOffsetLast = (TPC_Z0 + m * master[2]) * StarMagE;
6040 local[0] = m * OUTERGGFirst * tanSecPhi;
6041 local[1] = OUTERGGFirst;
6043 oAlignMatrix.LocalToMaster(local,master);
6048 oOffsetFirst = (TPC_Z0 + m * master[2]) * StarMagE;
6050 local[0] = m * OUTERGGLast * tanSecPhi;
6051 local[1] = OUTERGGLast;
6053 oAlignMatrix.LocalToMaster(local,master);
6058 oOffsetLast = (TPC_Z0 + m * master[2]) * StarMagE;
6066 oOffsetFirst = (Seclist[k] == 12 || Seclist[k] == 24 ?
6067 0.1 * StarMagE * (1.5 - Seclist[k]/24.) : 0);
6074 for ( Int_t i = 1 ; i < ROWS-1 ; i++ )
6076 Double_t Radius = IFCRadius + i*GRIDSIZER ;
6077 Double_t local_y = Radius * cosSecPhi ;
6079 if (local_y <= INNERGGFirst)
6080 tempV = iOffsetFirst*(Radius - IFCRadius)/(INNERGGFirst*secSecPhi - IFCRadius);
6081 else if (local_y <= INNERGGLast)
6082 tempV = iOffsetFirst + (iOffsetLast -iOffsetFirst)*(local_y-INNERGGFirst)/INNERGGSpan;
6083 else if (local_y <= OUTERGGFirst)
6084 tempV = iOffsetLast + (oOffsetFirst-iOffsetLast) *(local_y-INNERGGLast) /WIREGAP;
6085 else if (local_y <= OUTERGGLast)
6086 tempV = oOffsetFirst + (oOffsetLast -oOffsetFirst)*(local_y-OUTERGGFirst)/OUTERGGSpan;
6088 tempV = oOffsetLast*(OFCRadius - Radius)/(OFCRadius - OUTERGGLast*secSecPhi);
6089 ArrayV(i,COLUMNS-1) = tempV;
6096 Poisson3DRelaxation( ArrayofArrayV, ArrayofCharge, ArrayofEroverEz, ArrayofEPhioverEz, PHISLICES,
6097 GRIDSIZEPHI, ITERATIONS, 0) ;
6103 for ( Int_t k = 0 ; k < PHISLICES1 ; k++ )
6105 TMatrix &tiltEr = *ArrayoftiltEr[k] ;
6106 TMatrix &tiltEphi = *ArrayoftiltEphi[k] ;
6107 phi = Philist[k == PHISLICES ? 0 : k] ;
6108 for ( Int_t i = 0 ; i < EMap_nZ ; i++ )
6110 z = TMath::Abs(eZList[i]) ;
6111 TMatrix** ArrayofEroverEz = ( eZList[i] > 0 ? ArrayofEroverEzW : ArrayofEroverEzE ) ;
6112 TMatrix** ArrayofEPhioverEz = ( eZList[i] > 0 ? ArrayofEPhioverEzW : ArrayofEPhioverEzE ) ;
6113 for ( Int_t j = 0 ; j < neR3D ; j++ )
6116 tiltEr(j,i) = Interpolate3DTable( ORDER,r,z,phi,ROWS,COLUMNS,PHISLICES,Rlist,Zedlist,Philist,ArrayofEroverEz ) ;
6117 tiltEphi(j,i) = Interpolate3DTable( ORDER,r,z,phi,ROWS,COLUMNS,PHISLICES,Rlist,Zedlist,Philist,ArrayofEPhioverEz ) ;
6122 for ( Int_t k = 0 ; k < PHISLICES ; k++ )
6124 ArrayofArrayV[k] -> Delete() ;
6125 ArrayofCharge[k] -> Delete() ;
6126 ArrayofEroverEzW[k] -> Delete() ;
6127 ArrayofEPhioverEzW[k] -> Delete() ;
6128 ArrayofEroverEzE[k] -> Delete() ;
6129 ArrayofEPhioverEzE[k] -> Delete() ;
6134 if (usingCartesian) Cart2Polar(x,r,phi);
6135 else { r = x[0]; phi = x[1]; }
6136 if ( phi < 0 ) phi += TMath::TwoPi() ;
6137 z = LimitZ( Sector, x ) ;
6139 Er_integral = Interpolate3DTable( ORDER, r, z, phi, neR3D, EMap_nZ, PHISLICES1, eRadius, eZList, Philist, ArrayoftiltEr ) ;
6140 Ephi_integral = Interpolate3DTable( ORDER, r, z, phi, neR3D, EMap_nZ, PHISLICES1, eRadius, eZList, Philist, ArrayoftiltEphi ) ;
6144 phi = phi - ( Const_0*Ephi_integral - Const_1*Er_integral ) / r ;
6145 r = r - ( Const_0*Er_integral + Const_1*Ephi_integral ) ;
6148 if (usingCartesian) Polar2Cart(r,phi,Xprime);
6149 else { Xprime[0] = r; Xprime[1] = phi; }
6161 Int_t StMagUtilities::IterationFailCount()
6164 Int_t temp = iterationFailCounter;
6165 iterationFailCounter = 0;
6170 void StMagUtilities::BFieldTpc (
const Float_t xTpc[], Float_t BTpc[], Int_t Sector ) {
6171 if (StTpcDb::IsOldScheme()) {
6172 BField( xTpc, BTpc) ;
6175 Double_t Tpc[3] = {xTpc[0], xTpc[1], xTpc[2]};
6177 StTpcDb::instance()->Tpc2GlobalMatrix().LocalToMaster(Tpc,coorG);
6178 Float_t xyzG[3] = {(Float_t) coorG[0], (Float_t) coorG[1], (Float_t) coorG[2]};
6181 Double_t BGD[3] = {BG[0], BG[1], BG[2]};
6183 StTpcDb::instance()->Tpc2GlobalMatrix().MasterToLocalVect(BGD,BTpcL);
6190 void StMagUtilities::B3DFieldTpc (
const Float_t xTpc[], Float_t BTpc[], Int_t Sector ) {
6191 if (StTpcDb::IsOldScheme()) {
6192 B3DField( xTpc, BTpc) ;
6194 BFieldTpc(xTpc, BTpc, Sector);
virtual void UndoTwistDistortion(const Float_t x[], Float_t Xprime[], Int_t Sector=-1)
Twist distortion.
virtual void UndoFullGridLeakDistortion(const Float_t x[], Float_t Xprime[], Int_t Sector=-1)
Full GridLeak Distortion Calculation.
virtual void UndoIFCShiftDistortion(const Float_t x[], Float_t Xprime[], Int_t Sector=-1)
IFC Shift Distortion.
virtual void UndoClockDistortion(const Float_t x[], Float_t Xprime[], Int_t Sector=-1)
Clock distortion.
virtual void ApplySpaceChargeDistortion(const Double_t sc, const Int_t Charge, const Float_t x[3], const Float_t p[3], const Prime PrimaryOrGlobal, Int_t &new_Charge, Float_t x_new[3], Float_t p_new[3], const unsigned int RowMask1=0xFFFFFF00, const unsigned int RowMask2=0x1FFFFF, const Float_t VertexError=0.0200)
Apply the (1/R**2) space charge correction to selected data from the microDSTs.
virtual void UndoGridLeakDistortion(const Float_t x[], Float_t Xprime[], Int_t Sector=-1)
Grid Leakage entry function.
virtual void FixSpaceChargeDistortion(const Int_t Charge, const Float_t x[3], const Float_t p[3], const Prime PrimaryOrGlobal, Float_t x_new[3], Float_t p_new[3], const unsigned int RowMask1=0xFFFFFF00, const unsigned int RowMask2=0x1FFFFF, const Float_t VertexError=0.0200)
Convert from the old (Uniform) space charge correction to the new (1/R**2) space charge correction...
virtual void UndoBDistortion(const Float_t x[], Float_t Xprime[], Int_t Sector=-1)
B field distortions in 3D ( no Table ) - calculate the distortions due to the shape of the B field...
virtual void Undo2DBDistortion(const Float_t x[], Float_t Xprime[], Int_t Sector=-1)
2D - faster - B field distortions ( no Table ) - calculate the distortions due to the shape of the B ...
virtual void Undo3DGridLeakDistortion(const Float_t x[], Float_t Xprime[], Int_t Sector=-1)
3D GridLeak Distortion Calculation
virtual void FastUndoBDistortion(const Float_t x[], Float_t Xprime[], Int_t Sector=-1)
3D - B field distortions (Table) - calculate the distortions due to the shape of the B field ...
virtual void UndoSectorAlignDistortion(const Float_t x[], Float_t Xprime[], Int_t Sector=-1)
3D Sector Alignment Distortion Calculation
virtual void UndoGGVoltErrorDistortion(const Float_t x[], Float_t Xprime[], Int_t Sector=-1)
Gated Grid Voltage Error.
virtual void Undo2DGridLeakDistortion(const Float_t x[], Float_t Xprime[], Int_t Sector=-1)
Grid Leakage Calculation.
StMagUtilities(StTpcDb *dbin=0, Int_t mode=0)
StMagUtilities constructor using the DataBase.
virtual void UndoSpaceChargeDistortion(const Float_t x[], Float_t Xprime[], Int_t Sector=-1)
Space Charge entry function.
virtual Int_t PredictSpaceChargeDistortion(Int_t sec, Int_t Charge, Float_t Pt, Float_t VertexZ, Float_t PseudoRapidity, Float_t DCA, const unsigned int RowMask1, const unsigned int RowMask2, Float_t &pSpace)
PredictSpaceCharge - Input Physical-Signed DCA and get back spacecharge parameter plus a success or f...
virtual void UndoMembraneDistortion(const Float_t x[], Float_t Xprime[], Int_t Sector=-1)
Membrane distortion.
virtual void UndoEndcapDistortion(const Float_t x[], Float_t Xprime[], Int_t Sector=-1)
Endcap distortion.
virtual void DoDistortion(const Float_t x[], Float_t Xprime[], Int_t Sector=-1)
Main Entry Point for requests to DO the E and B field distortions (for simulations) ...
virtual void UndoAbortGapDistortion(const Float_t x[], Float_t Xprime[], Int_t Sector=-1, Float_t TimeSinceDeposition=-1.0)
Abort Gap Cleaning Cycle space charge correction.
virtual void UndoPad13Distortion(const Float_t x[], Float_t Xprime[], Int_t Sector=-1)
Pad row 13 distortion.
virtual Int_t GetSpaceChargeMode()
Space Charge Correction Mode.
virtual void UndoPad40Distortion(const Float_t x[], Float_t Xprime[], Int_t Sector=-1)
PadRow 40 and/or PadRow 13 distortion correction.
virtual void UndoShortedRingDistortion(const Float_t x[], Float_t Xprime[], Int_t Sector=-1)
Shorted Ring Distortion.
virtual void ManualShortedRing(Int_t EastWest, Int_t InnerOuter, Float_t RingNumber, Float_t MissingRValue, Float_t ExtraRValue)
Manually setup a shorted ring in the TPC.
virtual void UndoSpaceChargeR0Distortion(const Float_t x[], Float_t Xprime[], Int_t Sector=-1)
Space Charge Correction.
virtual void FastUndo2DBDistortion(const Float_t x[], Float_t Xprime[], Int_t Sector=-1)
2D - faster - B field distortions (Table) - calculate the distortions due to the shape of the B field...
virtual void UndoSpaceChargeR2Distortion(const Float_t x[], Float_t Xprime[], Int_t Sector=-1)
1/R**2 SpaceCharge Distortion
virtual void UndoSpaceChargeFXTDistortion(const Float_t x[], Float_t Xprime[], Int_t Sector=-1)
FixedTarget 1/r SpaceCharge Distortion.