diff --git a/src/THcHodoscope.cxx b/src/THcHodoscope.cxx index e9a8b1d032a92e1c64f717e576395b2b6a00dbb0..d00b98d95f8366a61da9200917d310b19f4fbdbe 100644 --- a/src/THcHodoscope.cxx +++ b/src/THcHodoscope.cxx @@ -15,6 +15,13 @@ /////////////////////////////////////////////////////////////////////////////// #include "THcSignalHit.h" +#include "THcShower.h" + +#include "THcHitList.h" +#include "THcRawShowerHit.h" +#include "TClass.h" +#include "math.h" +#include "THaSubDetector.h" #include "THcHodoscope.h" #include "THaEvData.h" @@ -104,6 +111,9 @@ void THcHodoscope::Setup(const char* name, const char* description) fPlaneNames[i] = new char[plane_names[i].length()]; strcpy(fPlaneNames[i], plane_names[i].c_str()); } + + //myShower = new THcShower("cal", "Shower" ); + /* fPlaneNames = new char* [fNPlanes]; for(Int_t i=0;i<fNPlanes;i++) {fPlaneNames[i] = new char[3];} // Should get the plane names from parameters. @@ -136,6 +146,23 @@ THaAnalysisObject::EStatus THcHodoscope::Init( const TDatime& date ) // But it needs to happen before the sub detectors are initialized // so that they can get the pointer to the hitlist. + // --------------- To get energy from THcShower ---------------------- + + const char* shower_detector_name = "cal"; + THaApparatus* app = GetApparatus(); + THaDetector* det = app->GetDetector( shower_detector_name ); + + if( !dynamic_cast<THcShower*>(det) ) { + Error("THcHodoscope", "Cannot find shower detector %s", + shower_detector_name ); + return fStatus = kInitError; + } + + fShower = static_cast<THcShower*>(det); // fShower is a membervariable + + // --------------- To get energy from THcShower ---------------------- + + InitHitList(fDetMap, "THcHodoscopeHit", 100); EStatus status; @@ -165,6 +192,50 @@ THaAnalysisObject::EStatus THcHodoscope::Init( const TDatime& date ) return kInitError; } + fKeepPos = new Bool_t [MAXHODHITS]; + fKeepNeg = new Bool_t [MAXHODHITS]; + fGoodTDCPos = new Bool_t [MAXHODHITS]; + fGoodTDCNeg = new Bool_t [MAXHODHITS]; + + fGoodRawPad = new Int_t [MAXHODHITS]; + fHitPaddle = new Int_t [MAXHODHITS]; + fNScinHit = new Int_t [MAXHODHITS]; + fNPmtHit = new Int_t [MAXHODHITS]; + fTimeHist = new Int_t [200]; + + fTimeAtFP = new Double_t [MAXHODHITS]; + fScinSigma = new Double_t [MAXHODHITS]; + fGoodScinTime = new Double_t [MAXHODHITS]; + fScinTime = new Double_t [MAXHODHITS]; + fTime = new Double_t [MAXHODHITS]; + adcPh = new Double_t [MAXHODHITS]; + fPath = new Double_t [MAXHODHITS]; + fTimePos = new Double_t [MAXHODHITS]; + fTimeNeg = new Double_t [MAXHODHITS]; + fScinTimefp = new Double_t [MAXHODHITS]; + fScinPosTime = new Double_t [MAXHODHITS]; + fScinNegTime = new Double_t [MAXHODHITS]; + + fNScinHits = new Int_t [fNPlanes]; + fGoodPlaneTime = new Bool_t [fNPlanes]; + fNPlaneTime = new Int_t [fNPlanes]; + fSumPlaneTime = new Double_t [fNPlanes]; + + // Double_t fHitCnt4 = 0., fHitCnt3 = 0.; + + Int_t m = 0; + + fdEdX = new Double_t*[MAXHODHITS]; + for ( m = 0; m < MAXHODHITS; m++ ){ + fdEdX[m] = new Double_t[MAXHODHITS]; + } + + fScinHit = new Double_t*[fNPlanes]; + for ( m = 0; m < fNPlanes; m++ ){ + fScinHit[m] = new Double_t[fNPaddle[0]]; + } + + return fStatus = kOK; } //_____________________________________________________________________________ @@ -290,7 +361,7 @@ void THcHodoscope::DefineArray(const char* fName, char** Suffix, const Int_t ind Int_t THcHodoscope::ReadDatabase( const TDatime& date ) { - MAXHODHITS = 53; + MAXHODHITS = 30; fBeta = new Double_t[ MAXHODHITS ]; fBetaChisq = new Double_t[ MAXHODHITS ]; @@ -369,32 +440,58 @@ Int_t THcHodoscope::ReadDatabase( const TDatime& date ) fHodoPosInvAdcAdc=new Double_t [fMaxHodoScin]; fHodoNegInvAdcAdc=new Double_t [fMaxHodoScin]; + + prefix[1]='\0'; DBRequest list[]={ - {"start_time_center", &fStartTimeCenter, kDouble}, - {"start_time_slop", &fStartTimeSlop, kDouble}, - {"scin_tdc_to_time", &fScinTdcToTime, kDouble}, - {"scin_tdc_min", &fScinTdcMin, kDouble}, - {"scin_tdc_max", &fScinTdcMax, kDouble}, - {"tof_tolerance", &fTofTolerance, kDouble,0,1}, - {"pathlength_central", &fPathLengthCentral, kDouble}, - {"hodo_vel_light",&fHodoVelLight[0],kDouble,fMaxHodoScin}, - {"hodo_pos_sigma",&fHodoPosSigma[0],kDouble,fMaxHodoScin}, - {"hodo_neg_sigma",&fHodoNegSigma[0],kDouble,fMaxHodoScin}, - {"hodo_pos_minph",&fHodoPosMinPh[0],kDouble,fMaxHodoScin}, - {"hodo_neg_minph",&fHodoNegMinPh[0],kDouble,fMaxHodoScin}, - {"hodo_pos_phc_coeff",&fHodoPosPhcCoeff[0],kDouble,fMaxHodoScin}, - {"hodo_neg_phc_coeff",&fHodoNegPhcCoeff[0],kDouble,fMaxHodoScin}, - {"hodo_pos_time_offset",&fHodoPosTimeOffset[0],kDouble,fMaxHodoScin}, - {"hodo_neg_time_offset",&fHodoNegTimeOffset[0],kDouble,fMaxHodoScin}, - {"hodo_pos_ped_limit",&fHodoPosPedLimit[0],kInt,fMaxHodoScin}, - {"hodo_neg_ped_limit",&fHodoNegPedLimit[0],kInt,fMaxHodoScin}, - {"tofusinginvadc",&fTofUsingInvAdc,kInt,0,1}, + // {"scin_2x_zpos", &fScin2XZpos, kDouble, 0, 1}, + // {"scin_2x_dzpos", &fScin2XdZpos, kDouble, 0, 1}, + // {"scin_2y_zpos", &fScin2YZpos, kDouble, 0, 1}, + // {"scin_2y_dzpos", &fScin2YdZpos, kDouble, 0, 1}, + // {"sel_betamin", &fSelBetaMin, kDouble, 0, 1}, + // {"sel_dedx1min", &fSeldEdX1Min, kDouble, 0, 1}, + // {"sel_dedx1max", &fSeldEdX1Max, kDouble, 0, 1}, + // {"sel_using_scin", &fSelUsingScin, kInt, 0, 1}, + // {"sel_ndegreesmin", &fSelNDegreesMin, kDouble, 0, 1}, + {"start_time_center", &fStartTimeCenter, kDouble}, + {"start_time_slop", &fStartTimeSlop, kDouble}, + {"scin_tdc_to_time", &fScinTdcToTime, kDouble}, + {"scin_tdc_min", &fScinTdcMin, kDouble}, + {"scin_tdc_max", &fScinTdcMax, kDouble}, + {"tof_tolerance", &fTofTolerance, kDouble, 0, 1}, + {"pathlength_central", &fPathLengthCentral, kDouble}, + {"hodo_vel_light", &fHodoVelLight[0], kDouble, fMaxHodoScin}, + {"hodo_pos_sigma", &fHodoPosSigma[0], kDouble, fMaxHodoScin}, + {"hodo_neg_sigma", &fHodoNegSigma[0], kDouble, fMaxHodoScin}, + {"hodo_pos_minph", &fHodoPosMinPh[0], kDouble, fMaxHodoScin}, + {"hodo_neg_minph", &fHodoNegMinPh[0], kDouble, fMaxHodoScin}, + {"hodo_pos_phc_coeff", &fHodoPosPhcCoeff[0], kDouble, fMaxHodoScin}, + {"hodo_neg_phc_coeff", &fHodoNegPhcCoeff[0], kDouble, fMaxHodoScin}, + {"hodo_pos_time_offset", &fHodoPosTimeOffset[0], kDouble, fMaxHodoScin}, + {"hodo_neg_time_offset", &fHodoNegTimeOffset[0], kDouble, fMaxHodoScin}, + {"hodo_pos_ped_limit", &fHodoPosPedLimit[0], kInt, fMaxHodoScin}, + {"hodo_neg_ped_limit", &fHodoNegPedLimit[0], kInt, fMaxHodoScin}, + {"tofusinginvadc", &fTofUsingInvAdc, kInt, 0, 1}, {0} }; fTofUsingInvAdc = 0; // Default if not defined fTofTolerance = 3.0; // Default if not defined gHcParms->LoadParmValues((DBRequest*)&list,prefix); + + + cout << "\n\n\n\n\n\nPaddles1x = " << fNPaddle[0] + << "\nscin_2y_zpos = " << fScin2YZpos + << "\nscin_2y_dzpos = " << fScin2YdZpos + << endl; + // << "\ndedx max = " << fSeldEdX1Max + // << "\nbeta min = " << fSelBetaMin + // << "\nbeta max = " << fSelBetaMax + // << "\net min = " << fSelEtMin + // << "\net max = " << fSelEtMax + + // << endl; + + if (fTofUsingInvAdc) { DBRequest list2[]={ {"hodo_pos_invadc_offset",&fHodoPosInvAdcOffset[0],kDouble,fMaxHodoScin}, @@ -480,7 +577,7 @@ Int_t THcHodoscope::DefineVariables( EMode mode ) // { "y_adc", "y-position from amplitudes (m)", "fYa" }, // { "time", "Time of hit at plane (s)", "fTime" }, // { "dtime", "Est. uncertainty of time (s)", "fdTime" }, - // { "dedx", "dEdX-like deposited in paddle", "fAmpl" }, + { "dedx", "dEdX-like deposited in paddle", "fdEdX" }, // { "troff", "Trigger offset for paddles", "fTrigOff"}, // { "trn", "Number of tracks for hits", "GetNTracks()" }, // { "trx", "x-position of track in det plane", "fTrackProj.THaTrackProj.fX" }, @@ -517,7 +614,17 @@ THcHodoscope::~THcHodoscope() void THcHodoscope::DeleteArrays() { // Delete member arrays. Used by destructor. - + Int_t k; + for( k = 0; k < MAXHODHITS; k++){ + delete [] fdEdX[k]; + } + delete [] fdEdX; + + for( k = 0; k < fNPlanes; k++){ + delete [] fScinHit[k]; + } + delete [] fScinHit; + delete [] fNPaddle; fNPaddle = NULL; delete [] fHodoVelLight; fHodoVelLight = NULL; delete [] fHodoPosSigma; fHodoPosSigma = NULL; @@ -536,12 +643,14 @@ void THcHodoscope::DeleteArrays() delete [] fHodoNegInvAdcLinear; fHodoNegInvAdcLinear = NULL; delete [] fHodoPosInvAdcAdc; fHodoPosInvAdcAdc = NULL; + delete [] fGoodRawPad; fGoodRawPad = NULL; // Ahmed delete [] fHitPaddle; fHitPaddle = NULL; // Ahmed delete [] fNScinHit; fNScinHit = NULL; // Ahmed delete [] fNPmtHit; fNPmtHit = NULL; // Ahmed delete [] fTimeHist; fTimeHist = NULL; // Ahmed delete [] fTimeAtFP; fTimeAtFP = NULL; // Ahmed + delete [] fScinSigma; fScinSigma = NULL; // Ahmed delete [] fGoodScinTime; fGoodScinTime = NULL; // Ahmed delete [] fScinTime; fScinTime = NULL; // Ahmed @@ -561,6 +670,8 @@ void THcHodoscope::DeleteArrays() delete [] fNPlaneTime; fNPlaneTime = NULL; // Ahmed delete [] fSumPlaneTime; fSumPlaneTime = NULL; // Ahmed + delete [] fNScinHits; fNScinHits = NULL; // Ahmed + // delete [] fSpacing; fSpacing = NULL; //delete [] fCenter; fCenter = NULL; // This 2D. What is correct way to delete? @@ -616,8 +727,10 @@ Int_t THcHodoscope::Decode( const THaEvData& evdata ) // // GN: print event number so we can cross-check with engine // if (evdata.GetEvNum()>1000) - // cout <<"\nhcana event = " << evdata.GetEvNum()<<endl; + // cout <<"\nhcana_event " << evdata.GetEvNum()<<endl; + fCheckEvent = evdata.GetEvNum(); + if(gHaCuts->Result("Pedestal_event")) { Int_t nexthit = 0; for(Int_t ip=0;ip<fNPlanes;ip++) { @@ -648,6 +761,15 @@ Int_t THcHodoscope::Decode( const THaEvData& evdata ) fNfptimes=0; for(Int_t ip=0;ip<fNPlanes;ip++) { + if ( ip == 2 ){ + fHodoCenter3 = fPlanes[ip]->GetPosCenter(0) + fPlanes[ip]->GetPosOffset(); + fScin2XSpacing = fPlanes[ip]->GetSpacing(); + } + if ( ip == 3 ){ + fHodoCenter4 = fPlanes[ip]->GetPosCenter(0) + fPlanes[ip]->GetPosOffset(); + fScin2YSpacing = fPlanes[ip]->GetSpacing(); + } + // if ( !fPlanes[ip] ) // Ahmed // return -1; // Ahmed @@ -703,47 +825,34 @@ Double_t THcHodoscope::TimeWalkCorrection(const Int_t& paddle, //_____________________________________________________________________________ Int_t THcHodoscope::CoarseProcess( TClonesArray& tracks ) +{ + + ApplyCorrections(); + + return 0; +} + +//_____________________________________________________________________________ +Int_t THcHodoscope::FineProcess( TClonesArray& tracks ) { Int_t Ntracks = tracks.GetLast()+1; // Number of reconstructed tracks - Int_t fPaddle = 0, fIndex, k, fJMax, fMaxHit, ip, ihit; - Int_t fNfpTime, fSumfpTime, fNScinHits; + Int_t fPaddle = 0, fIndex, k, fJMax, fMaxHit, ip, ihit, itrack; + Int_t fNfpTime, fSumfpTime, fRawIndex = -1; //, fNScinHits[fNPlanes]; Double_t fScinTrnsCoord, fScinLongCoord, fScinCenter; Double_t fP, fBetaP, fXcoord, fYcoord, fTMin; - // Int_t fNtof, fNtofPairs; - // ------------------------------------------------- - - fKeepPos = new Bool_t [MAXHODHITS]; - fKeepNeg = new Bool_t [MAXHODHITS]; - fGoodTDCPos = new Bool_t [MAXHODHITS]; - fGoodTDCNeg = new Bool_t [MAXHODHITS]; - - fHitPaddle = new Int_t [MAXHODHITS]; - fNScinHit = new Int_t [MAXHODHITS]; - fNPmtHit = new Int_t [MAXHODHITS]; - fTimeHist = new Int_t [MAXHODHITS]; - fTimeAtFP = new Double_t [MAXHODHITS]; - - fScinSigma = new Double_t [MAXHODHITS]; - fGoodScinTime = new Double_t [MAXHODHITS]; - fScinTime = new Double_t [MAXHODHITS]; - fTime = new Double_t [MAXHODHITS]; - adcPh = new Double_t [MAXHODHITS]; - - fPath = new Double_t [MAXHODHITS]; - fTimePos = new Double_t [MAXHODHITS]; - fTimeNeg = new Double_t [MAXHODHITS]; - fScinTimefp = new Double_t [MAXHODHITS]; - - fTimeHist = new Int_t [200]; - // ------------------------------------------------- Double_t hpartmass=0.00051099; // Fix it for ( Int_t m = 0; m < MAXHODHITS; m++ ){ + + for(k = 0; k < MAXHODHITS; k++){ + fdEdX[m][k] = 0.; + } + fGoodRawPad[m] = 0; fScinSigma[m] = 0.; fHitPaddle[m] = 0.; fScinTime[m] = 0.; @@ -764,24 +873,25 @@ Int_t THcHodoscope::CoarseProcess( TClonesArray& tracks ) } - if (tracks.GetLast()+1 > 0 ) { // **MAIN LOOP: Loop over all tracks and get corrected time, tof, beta... - for ( Int_t itrack = 0; itrack < Ntracks; itrack++ ) { // Line 133 + for ( itrack = 0; itrack < Ntracks; itrack++ ) { // Line 133 THaTrack* theTrack = dynamic_cast<THaTrack*>( tracks.At(itrack) ); if (!theTrack) return -1; - fGoodPlaneTime = new Bool_t[4]; - for ( ip = 0; ip < fNPlanes; ip++ ){ fGoodPlaneTime[ip] = kFALSE; } - + for ( ip = 0; ip < fNPlanes; ip++ ){ + fGoodPlaneTime[ip] = kFALSE; + fNScinHits[ip] = 0; + } + fNfpTime = 0; fBetaChisq[itrack] = -3; fTimeAtFP[itrack] = 0.; fSumfpTime = 0.; // Line 138 fNScinHit[itrack] = 0; // Line 140 - fP = theTrack->GetP(); // Line 142 Fix it + fP = theTrack->GetP(); // Line 142 fBetaP = fP/( TMath::Sqrt( fP * fP + hpartmass * hpartmass) ); //! Calculate all corrected hit times and histogram @@ -799,12 +909,20 @@ Int_t THcHodoscope::CoarseProcess( TClonesArray& tracks ) for (Int_t j=0; j<200; j++) { fTimeHist[j]=0; } // Line 176 - fNPlaneTime = new Int_t [4]; - fSumPlaneTime = new Double_t [4]; + // fNPlaneTime = new Int_t [4]; + // fSumPlaneTime = new Double_t [4]; for ( ip = 0; ip < fNPlanes; ip++ ){ fNPlaneTime[ip] = 0; fSumPlaneTime[ip] = 0.; + // if ( ip == 2 ){ + // fHodoCenter3 = fPlanes[ip]->GetPosCenter(0) + fPlanes[ip]->GetPosOffset(); + // fScin2XSpacing = fPlanes[ip]->GetSpacing(); + // } + // if ( ip == 3 ){ + // fHodoCenter4 = fPlanes[ip]->GetPosCenter(0) + fPlanes[ip]->GetPosOffset(); + // fScin2YSpacing = fPlanes[ip]->GetSpacing(); + // } } // Loop over scintillator planes. // In ENGINE, its loop over good scintillator hits. @@ -812,10 +930,10 @@ Int_t THcHodoscope::CoarseProcess( TClonesArray& tracks ) fGoodTimeIndex = -1; for( ip = 0; ip < fNPlanes; ip++ ) { - fNScinHits = fPlanes[ip]->GetNScinHits(); + fNScinHits[ip] = fPlanes[ip]->GetNScinHits(); // first loop over hits with in a single plane - for ( ihit = 0; ihit < fNScinHits; ihit++ ){ + for ( ihit = 0; ihit < fNScinHits[ip]; ihit++ ){ fTimePos[ihit] = -99.; fTimeNeg[ihit] = -99.; @@ -917,7 +1035,7 @@ Int_t THcHodoscope::CoarseProcess( TClonesArray& tracks ) if ( fJMax >= 0 ){ // Line 248. Here I followed the code of THcSCintilaltorPlane::PulseHeightCorrection fTMin = 0.5 * fJMax; - for( ihit = 0; ihit < fNScinHits; ihit++) { // Loop over sinc. hits. in plane + for( ihit = 0; ihit < fNScinHits[ip]; ihit++) { // Loop over sinc. hits. in plane if ( ( fTimePos[ihit] > fTMin ) && ( fTimePos[ihit] < ( fTMin + fTofTolerance ) ) ) { fKeepPos[ihit]=kTRUE; } @@ -927,8 +1045,8 @@ Int_t THcHodoscope::CoarseProcess( TClonesArray& tracks ) } } // fJMax > 0 condition - fScinPosTime = new Double_t [MAXHODHITS]; - fScinNegTime = new Double_t [MAXHODHITS]; + // fScinPosTime = new Double_t [MAXHODHITS]; + // fScinNegTime = new Double_t [MAXHODHITS]; for ( k = 0; k < MAXHODHITS; k++ ){ fScinPosTime[k] = 0; @@ -939,12 +1057,15 @@ Int_t THcHodoscope::CoarseProcess( TClonesArray& tracks ) // ---------------------- Scond loop over scint. hits in a plane ------------------------------ //--------------------------------------------------------------------------------------------- - for ( ihit = 0; ihit < fNScinHits; ihit++ ){ + for ( ihit = 0; ihit < fNScinHits[ip]; ihit++ ){ fGoodTimeIndex ++; + fRawIndex ++; fPaddle = ((THcSignalHit*)scinPosTDC->At(ihit))->GetPaddleNumber()-1; fHitPaddle[fGoodTimeIndex] = fPaddle; + + fGoodRawPad[fRawIndex] = fPlanes[ip]->GetGoodRawPadNum(ihit); fXcoord = theTrack->GetX() + theTrack->GetTheta() * ( fPlanes[ip]->GetZpos() + ( fPaddle % 2 ) * fPlanes[ip]->GetDzpos() ); // Line 277 @@ -1069,7 +1190,6 @@ Int_t THcHodoscope::CoarseProcess( TClonesArray& tracks ) // scinHit[itrack][fNScinHit[itrack]] = ihit; // scinfFPTime[itrack][fNScinHit[itrack]] = fScinTimefp[ihit]; - // --------------------------------------------------------------------------- // Date: July 8 2014 // This counts the pmt hits. Right now we don't need it so it is commentd off @@ -1083,29 +1203,34 @@ Int_t THcHodoscope::CoarseProcess( TClonesArray& tracks ) // --------------------------------------------------------------------------- + // fdEdX[itrack][ihit] = 5.0; + + // -------------------------------------------------------------------------------------------- // Date: July 8 201 May be we need this, not sure. // - // if ( fGoodTDCPos[itrack][ihit] ){ - // if ( fGoodTDCNeg[itrack][ihit] ){ - // dedX[itrack][fNScinHit[itrack]] = - // TMath::Sqrt( TMath::Max( 0., ((THcSignalHit*)scinPosADC->At(ihit))->GetData() * - // ((THcSignalHit*)scinNegADC->At(ihit))->GetData() ) ); - // } - // else{ - // dedX[itrack][fNScinHit[itrack]] = - // TMath::Max( 0., ((THcSignalHit*)scinPosADC->At(ihit))->GetData() ); - // } - // } - // else{ - // if ( fGoodTDCNeg[itrack][ihit] ){ - // dedX[itrack][fNScinHit[itrack]] = - // TMath::Max( 0., ((THcSignalHit*)scinNegADC->At(ihit))->GetData() ); - // } - // else{ - // dedX[itrack][fNScinHit[itrack]] = 0.; - // } - // } + + if ( fGoodTDCPos[fGoodTimeIndex] ){ + if ( fGoodTDCNeg[fGoodTimeIndex] ){ + + fdEdX[itrack][fNScinHit[itrack]-1] = + TMath::Sqrt( TMath::Max( 0., ((THcSignalHit*)scinPosADC->At(ihit))->GetData() * + ((THcSignalHit*)scinNegADC->At(ihit))->GetData() ) ); + } + else{ + fdEdX[itrack][fNScinHit[itrack]-1] = + TMath::Max( 0., ((THcSignalHit*)scinPosADC->At(ihit))->GetData() ); + } + } + else{ + if ( fGoodTDCNeg[fGoodTimeIndex] ){ + fdEdX[itrack][fNScinHit[itrack]-1] = + TMath::Max( 0., ((THcSignalHit*)scinNegADC->At(ihit))->GetData() ); + } + else{ + fdEdX[itrack][fNScinHit[itrack]-1] = 0.; + } + } // -------------------------------------------------------------------------------------------- @@ -1144,8 +1269,8 @@ Int_t THcHodoscope::CoarseProcess( TClonesArray& tracks ) if (!fPlanes[ip]) return -1; - fNScinHits = fPlanes[ip]->GetNScinHits(); - for (ihit = 0; ihit < fNScinHits; ihit++ ){ + fNScinHits[ip] = fPlanes[ip]->GetNScinHits(); + for (ihit = 0; ihit < fNScinHits[ip]; ihit++ ){ fGoodTimeIndex ++; if ( fGoodScinTime[fGoodTimeIndex] ) { @@ -1178,8 +1303,8 @@ Int_t THcHodoscope::CoarseProcess( TClonesArray& tracks ) if (!fPlanes[ip]) return -1; - fNScinHits = fPlanes[ip]->GetNScinHits(); - for (ihit = 0; ihit < fNScinHits; ihit++ ){ // Loop over hits of a plane + fNScinHits[ip] = fPlanes[ip]->GetNScinHits(); + for (ihit = 0; ihit < fNScinHits[ip]; ihit++ ){ // Loop over hits of a plane fGoodTimeIndex ++; if ( fGoodScinTime[fGoodTimeIndex] ){ @@ -1246,19 +1371,20 @@ Int_t THcHodoscope::CoarseProcess( TClonesArray& tracks ) } } - } // Main loop over tracks ends here. - } // If condition for at least one track - ApplyCorrections(); - - return 0; -} + // fBetaChisq[itrack] + // fFPTime[ind] + + theTrack->SetDedx(fdEdX[itrack][0]); + theTrack->SetBeta(fBeta[itrack]); + theTrack->SetBetaChi2( fBetaChisq[itrack] ); + + } // Main loop over tracks ends here. + + } // If condition for at least one track -//_____________________________________________________________________________ -Int_t THcHodoscope::FineProcess( TClonesArray& tracks ) -{ - return 0; + } //_____________________________________________________________________________ Int_t THcHodoscope::GetScinIndex( Int_t nPlane, Int_t nPaddle ) { diff --git a/src/THcHodoscope.h b/src/THcHodoscope.h index aba78e9c642d193539fdcfe77065d91e17734103..87c9f5a32e82c0b45f0693687b8bc0f95f13fa85 100644 --- a/src/THcHodoscope.h +++ b/src/THcHodoscope.h @@ -14,6 +14,21 @@ #include "THcHitList.h" #include "THcHodoscopeHit.h" #include "THcScintillatorPlane.h" +#include "THcShower.h" + +#include "THaTrackingDetector.h" +#include "THcHitList.h" +#include "THcRawDCHit.h" +#include "THcSpacePoint.h" +#include "THcDriftChamberPlane.h" +#include "THcDriftChamber.h" +#include "TMath.h" + +#include "THaSubDetector.h" +#include "TClonesArray.h" +#include <iostream> +#include <fstream> + class THaScCalib; @@ -54,13 +69,22 @@ public: Double_t GetStartTimeSlop() const {return fStartTimeSlop;} Double_t GetBetaNotrk() const {return fBetaNotrk;} + Int_t GetGoodRawPad(Int_t iii){return fGoodRawPad[iii];} + Double_t GetNScinHits(Int_t iii){return fNScinHits[iii];} + Double_t GetHodoCenter3() { return fHodoCenter3;} + Double_t GetHodoCenter4() { return fHodoCenter4;} + Double_t GetScin2XSpacing() { return fScin2XSpacing;} + Double_t GetScin2YSpacing() { return fScin2YSpacing;} + // Double_t GetBeta() const {return fBeta[];} Double_t GetBeta(Int_t iii) const {return fBeta[iii];} // Ahmed + Int_t GetEvent(){ return fCheckEvent;} Double_t GetHodoPosSigma(Int_t iii) const {return fHodoPosSigma[iii];} Double_t GetHodoNegSigma(Int_t iii) const {return fHodoNegSigma[iii];} + const TClonesArray* GetTrackHits() const { return fTrackProj; } friend class THaScCalib; @@ -79,7 +103,7 @@ protected: Double_t fStartTime; Int_t fNfptimes; - Double_t fdEdX; + // Double_t fdEdX; Double_t fBetaNotrk; // Double_t fBeta; @@ -121,19 +145,43 @@ protected: //-------------------------- Ahmed ----------------------------- - Int_t MAXHODHITS; - - Double_t* fTestArr; // [MAXHODHITS] Array + THcShower* fShower; + + Int_t fCheckEvent; + + Int_t fGoodTrack; + Int_t MAXHODHITS; + // Int_t fSelUsingScin; + Int_t fSelNDegreesMin; + Double_t fSeldEdX1Min; + Double_t fSeldEdX1Max; + Double_t fSelBetaMin; + Double_t fSelBetaMax; + Double_t fSelEtMin; + Double_t fSelEtMax; + Double_t fScin2XZpos; + Double_t fScin2XdZpos; + Double_t fScin2YZpos; + Double_t fScin2YdZpos; + + Double_t fChi2Min; + Double_t fHodoCenter4, fHodoCenter3; + Double_t fScin2YSpacing, fScin2XSpacing; + + Double_t** fdEdX; // [MAXHODHITS] Array + Double_t** fScinHit; // [fNPlanes] Array + Int_t* fGoodRawPad; Double_t* fBeta; // [MAXHODHITS] Array Double_t* fBetaChisq; // [MAXHODHITS] Array Double_t* fFPTime; // [fNPlanes] Array + Double_t* fScinSigma; Double_t* fGoodScinTime; Double_t* fScinTime; Double_t* fTime; - Double_t* adcPh; + Double_t* adcPh; // Correct it Double_t* fTimeAtFP; Double_t* fPath; Double_t* fTimePos; @@ -145,6 +193,7 @@ protected: Int_t* fHitPaddle; Int_t* fNScinHit; + Int_t* fNScinHits; Int_t* fNPmtHit; Int_t* fTimeHist; Int_t* fNPlaneTime; @@ -163,6 +212,10 @@ protected: TClonesArray* scinPosTDC; TClonesArray* scinNegTDC; + //test array + Double_t test_arr[53]; + Double_t test_arr1[2]; + //---------------------------------------------------------------- // Useful derived quantities diff --git a/src/THcScintillatorPlane.cxx b/src/THcScintillatorPlane.cxx index e945e5105a40c10248255df7f7163b072740ff23..f85baadd42f957141c9a7450c7bd8751840f84f0 100644 --- a/src/THcScintillatorPlane.cxx +++ b/src/THcScintillatorPlane.cxx @@ -45,6 +45,9 @@ THcScintillatorPlane::THcScintillatorPlane( const char* name, fNScinHits = 0; // fMaxHits=53; + + fGoodRawPadNum = new Int_t [fMaxHits]; // Ahmed + fpTimes = new Double_t [fMaxHits]; fScinTime = new Double_t [fMaxHits]; fScinSigma = new Double_t [fMaxHits]; @@ -72,6 +75,9 @@ THcScintillatorPlane::THcScintillatorPlane( const char* name, fNScinHits = 0; // fMaxHits=53; + + fGoodRawPadNum = new Int_t [fMaxHits]; + fpTimes = new Double_t [fMaxHits]; fScinTime = new Double_t [fMaxHits]; fScinSigma = new Double_t [fMaxHits]; @@ -94,6 +100,10 @@ THcScintillatorPlane::~THcScintillatorPlane() delete fScinTime; delete fScinSigma; delete fScinZpos; + + delete fGoodRawPadNum; + + } //______________________________________________________________________________ @@ -303,20 +313,35 @@ Int_t THcScintillatorPlane::ProcessHits(TClonesArray* rawhits, Int_t nexthit) mintdc=((THcHodoscope *)GetParent())->GetTdcMin(); maxtdc=((THcHodoscope *)GetParent())->GetTdcMax(); Int_t ihit = nexthit; + + // cout << "THcScintillatorPlane: raw htis = " << nrawhits << endl; + while(ihit < nrawhits) { THcHodoscopeHit* hit = (THcHodoscopeHit *) rawhits->At(ihit); if(hit->fPlane > fPlaneNum) { break; } Int_t padnum=hit->fCounter; + Int_t index=padnum-1; - if (hit->fTDC_pos > 0) ((THcSignalHit*) frPosTDCHits->ConstructedAt(nrPosTDCHits++))->Set(padnum, hit->fTDC_pos); - if (hit->fTDC_neg > 0) ((THcSignalHit*) frNegTDCHits->ConstructedAt(nrNegTDCHits++))->Set(padnum, hit->fTDC_neg); - if ((hit->fADC_pos-fPosPed[index]) >= 50) ((THcSignalHit*) frPosADCHits->ConstructedAt(nrPosADCHits++))->Set(padnum, hit->fADC_pos-fPosPed[index]); - if ((hit->fADC_neg-fNegPed[index]) >= 50) ((THcSignalHit*) frNegADCHits->ConstructedAt(nrNegADCHits++))->Set(padnum, hit->fADC_neg-fNegPed[index]); + if (hit->fTDC_pos > 0) + ((THcSignalHit*) frPosTDCHits->ConstructedAt(nrPosTDCHits++))->Set(padnum, hit->fTDC_pos); + if (hit->fTDC_neg > 0) + ((THcSignalHit*) frNegTDCHits->ConstructedAt(nrNegTDCHits++))->Set(padnum, hit->fTDC_neg); + if ((hit->fADC_pos-fPosPed[index]) >= 50) + ((THcSignalHit*) frPosADCHits->ConstructedAt(nrPosADCHits++))->Set(padnum, hit->fADC_pos-fPosPed[index]); + if ((hit->fADC_neg-fNegPed[index]) >= 50) + ((THcSignalHit*) frNegADCHits->ConstructedAt(nrNegADCHits++))->Set(padnum, hit->fADC_neg-fNegPed[index]); // check TDC values if (((hit->fTDC_pos >= mintdc) && (hit->fTDC_pos <= maxtdc)) || ((hit->fTDC_neg >= mintdc) && (hit->fTDC_neg <= maxtdc))) { + + fGoodRawPadNum[fNScinHits] = hit->fCounter; + + // cout << " plane = " << fPlaneNum << " hit = " << fNScinHits + 1 + // << " index = " << fNScinHits + // << " raw pad = " << fGoodRawPadNum[fNScinHits] << endl; + //TDC positive hit THcSignalHit *sighit = (THcSignalHit*) fPosTDCHits->ConstructedAt(nPosTDCHits++); sighit->Set(padnum, hit->fTDC_pos); @@ -339,6 +364,8 @@ Int_t THcScintillatorPlane::ProcessHits(TClonesArray* rawhits, Int_t nexthit) ihit++; } + // cout << "THcScintillatorPlane: ihit = " << ihit << endl; + return(ihit); } //________________________________________________________________________________ @@ -398,6 +425,9 @@ Int_t THcScintillatorPlane::PulseHeightCorrection() postime[i]=((THcSignalHit*) fPosTDCHits->At(i))->GetData()*tdctotime; j=((THcSignalHit*)fPosTDCHits->At(i))->GetPaddleNumber()-1; index=((THcHodoscope *)GetParent())->GetScinIndex(fPlaneNum-1,j); + + // cout << "THcScintillatorPlane: index = " << index << endl; + postime[i]=postime[i]-((THcHodoscope *)GetParent())->GetHodoPosPhcCoeff(index)* TMath::Sqrt(TMath::Max(0.,(pos_ph[i]/((THcHodoscope *)GetParent())->GetHodoPosMinPh(index)-1))); postime[i]=postime[i]-((THcHodoscope *)GetParent())->GetHodoPosTimeOffset(index); diff --git a/src/THcScintillatorPlane.h b/src/THcScintillatorPlane.h index 959038f6f3c926e468ee8e8f832d13e2fdc82076..2ecae143556e00169e1e99d4e307b3998d2edcf9 100644 --- a/src/THcScintillatorPlane.h +++ b/src/THcScintillatorPlane.h @@ -63,15 +63,17 @@ class THcScintillatorPlane : public THaSubDetector { TClonesArray* fParentHitList; - TClonesArray* GetPosADC() { return fPosADCHits;}; // Ahmed - TClonesArray* GetNegADC() { return fNegADCHits;}; // Ahmed - TClonesArray* GetPosTDC() { return fPosTDCHits;}; // Ahmed - TClonesArray* GetNegTDC() { return fNegTDCHits;}; // Ahmed - + TClonesArray* GetPosADC() { return fPosADCHits;}; // Ahmed + TClonesArray* GetNegADC() { return fNegADCHits;}; // Ahmed + TClonesArray* GetPosTDC() { return fPosTDCHits;}; // Ahmed + TClonesArray* GetNegTDC() { return fNegTDCHits;}; // Ahmed + Int_t GetGoodRawPadNum(Int_t index) {return fGoodRawPadNum[index];}; protected: + Int_t *fGoodRawPadNum; // array + TClonesArray* frPosTDCHits; TClonesArray* frNegTDCHits; TClonesArray* frPosADCHits;