00001 #include "GaudiKernel/MsgStream.h"
00002 #include "GaudiKernel/StatusCode.h"
00003 #include "GaudiKernel/INTupleSvc.h"
00004 #include "GaudiKernel/SmartDataPtr.h"
00005
00006 #include "EventModel/EventHeader.h"
00007 #include "EvTimeEvent/RecEsTime.h"
00008 #include "MdcRecEvent/RecMdcHit.h"
00009 #include "MdcRecEvent/RecMdcTrack.h"
00010 #include "MdcRecEvent/RecMdcKalTrack.h"
00011 #include "MdcRecEvent/RecMdcKalHelixSeg.h"
00012 #include "MdcRecEvent/RecMdcDedx.h"
00013 #include "MdcRecEvent/RecMdcDedxHit.h"
00014 #include "Identifier/Identifier.h"
00015 #include "Identifier/MdcID.h"
00016
00017 #include <TMath.h>
00018
00019 #include "DedxCalibAlg/DedxCalibEvent.h"
00020
00021 typedef std::vector<int> Vint;
00022
00023 using namespace std;
00024 using CLHEP::HepVector;
00025
00026
00027 DedxCalibEvent::DedxCalibEvent(const std::string& name, ISvcLocator* pSvcLocator): DedxCalib(name, pSvcLocator) {}
00028
00029
00030 void DedxCalibEvent::initializing()
00031 {
00032 MsgStream log(msgSvc(), name());
00033 log << MSG::INFO << "DedxCalibEvent::initializing()" << endreq;
00034
00035 StatusCode status;
00036 NTuplePtr nt1(ntupleSvc(),"FILE100/n103");
00037 if ( nt1 )
00038 m_nt1 = nt1;
00039 else
00040 {
00041 m_nt1= ntupleSvc()->book("FILE100/n103",CLID_ColumnWiseTuple,"dEdx per track");
00042 if ( m_nt1 )
00043 {
00044 status = m_nt1->addItem("ptrk",m_ptrk);
00045 status = m_nt1->addItem("ptrk_t",m_ptrk_t);
00046 status = m_nt1->addItem("sintheta",m_sintheta);
00047 status = m_nt1->addItem("costheta",m_costheta);
00048 status = m_nt1->addItem("charge",m_charge);
00049 status = m_nt1->addItem("runNO",m_runNO);
00050 status = m_nt1->addItem("runFlag",m_runFlag);
00051 status = m_nt1->addItem("evtNO",m_evtNO);
00052 status = m_nt1->addItem("t0",m_t0);
00053 status = m_nt1->addItem("trackId",m_trackId);
00054 status = m_nt1->addItem("poca_x",m_poca_x);
00055 status = m_nt1->addItem("poca_y",m_poca_y);
00056 status = m_nt1->addItem("poca_z",m_poca_z);
00057 status = m_nt1->addItem("recalg",m_recalg);
00058 status = m_nt1->addItem("nhit",m_nhit);
00059 status = m_nt1->addItem("nhits",m_nhits);
00060 status = m_nt1->addItem("usedhit",m_usedhit);
00061
00062 status = m_nt1->addItem("ndedxhit",m_nphlisthit,0,100);
00063 status = m_nt1->addIndexedItem("dEdx_hit",m_nphlisthit,m_dEdx_hit);
00064 status = m_nt1->addIndexedItem("pathlength_hit",m_nphlisthit,m_pathlength_hit);
00065 status = m_nt1->addIndexedItem("wid_hit",m_nphlisthit,m_wid_hit);
00066 status = m_nt1->addIndexedItem("layid_hit",m_nphlisthit,m_layid_hit);
00067 status = m_nt1->addIndexedItem("dd_in_hit",m_nphlisthit,m_dd_in_hit);
00068 status = m_nt1->addIndexedItem("eangle_hit",m_nphlisthit,m_eangle_hit);
00069 status = m_nt1->addIndexedItem("zhit_hit",m_nphlisthit,m_zhit_hit);
00070
00071
00072 status = m_nt1->addItem("dEdx_meas", m_dEdx_meas);
00073
00074
00075
00076 status = m_nt1->addItem("type",m_parttype);
00077 status = m_nt1->addItem("chidedx_e",m_chidedxe);
00078 status = m_nt1->addItem("chidedx_mu",m_chidedxmu);
00079 status = m_nt1->addItem("chidedx_pi",m_chidedxpi);
00080 status = m_nt1->addItem("chidedx_k",m_chidedxk);
00081 status = m_nt1->addItem("chidedx_p",m_chidedxp);
00082 status = m_nt1->addItem("partid",5,m_probpid);
00083 status = m_nt1->addItem("expectid",5,m_expectid);
00084 status = m_nt1->addItem("sigmaid",5,m_sigmaid);
00085 }
00086 }
00087
00088 NTuplePtr nt2(ntupleSvc(),"FILE100/n102");
00089 if ( nt2 ) m_nt2 = nt2;
00090 else
00091 {
00092 m_nt2= ntupleSvc()->book("FILE100/n102",CLID_RowWiseTuple,"dE/dx per hit");
00093 if ( m_nt2 )
00094 {
00095 status = m_nt2->addItem("charge",m_charge1);
00096 status = m_nt2->addItem("adc_raw",m_phraw);
00097 status = m_nt2->addItem("exraw",m_exraw);
00098 status = m_nt2->addItem("runNO",m_runNO1);
00099 status = m_nt2->addItem("runFlag",m_runFlag1);
00100 status = m_nt2->addItem("wire",m_wire);
00101 status = m_nt2->addItem("doca_in",m_doca_in);
00102 status = m_nt2->addItem("doca_ex",m_doca_ex);
00103 status = m_nt2->addItem("driftdist",m_driftdist);
00104 status = m_nt2->addItem("eangle",m_eangle);
00105 status = m_nt2->addItem("zhit",m_zhit);
00106 status = m_nt2->addItem("costheta1",m_costheta1);
00107 status = m_nt2->addItem("path_rphi",m_pathL);
00108 status = m_nt2->addItem("layer",m_layer);
00109 status = m_nt2->addItem("ptrk1",m_ptrk1);
00110 status = m_nt2->addItem("ptrk_hit",m_ptrk_hit);
00111 status = m_nt2->addItem("t01",m_t01);
00112 status = m_nt2->addItem("tdc_raw",m_tdc_raw);
00113 status = m_nt2->addItem("driftT",m_driftT);
00114 status = m_nt2->addItem("localwid",m_localwid);
00115 status = m_nt2->addItem("trackId1",m_trackId1);
00116 }
00117 }
00118 }
00119
00120 void DedxCalibEvent::genNtuple()
00121 {
00122 MsgStream log(msgSvc(), name());
00123 log << MSG::INFO << "DedxCalibEvent::genNtuple()" << endreq;
00124
00125 SmartDataPtr<Event::EventHeader> eventHeader(eventSvc(),"/Event/EventHeader");
00126 if (!eventHeader)
00127 {
00128 log << MSG::INFO << "Could not find Event Header" << endreq;
00129 return;
00130 }
00131 int eventNO = eventHeader->eventNumber();
00132 int runNO = eventHeader->runNumber();
00133 log << MSG::INFO << "now begin the event: runNO= "<<runNO<<" eventNO= "<< eventNO<< endreq;
00134
00135 int runFlag=0;
00136 if(runNO<RUN0) runFlag =0;
00137 if(runNO>=RUN1 && runNO<RUN2) runFlag =1;
00138 if(runNO>=RUN2 && runNO<RUN3) runFlag =2;
00139 if(runNO>=RUN3 && runNO<RUN4) runFlag =3;
00140 if(runNO>=RUN4 && runNO<RUN5) runFlag =4;
00141 if(runNO>=RUN5 && runNO<RUN6) runFlag =5;
00142 if(runNO>=RUN6 && runNO<RUN7) runFlag =6;
00143 if(runNO>=RUN7) runFlag =7;
00144
00145 double tes = -999.0;
00146 int esTimeflag = -1;
00147 SmartDataPtr<RecEsTimeCol> aevtimeCol(eventSvc(),"/Event/Recon/RecEsTimeCol");
00148 if( (!aevtimeCol) || (aevtimeCol->size()==0) ){
00149 tes = -9999.0;
00150 esTimeflag = -1;
00151 }else{
00152 RecEsTimeCol::iterator iter_evt = aevtimeCol->begin();
00153 for(; iter_evt!=aevtimeCol->end(); iter_evt++){
00154 tes = (*iter_evt)->getTest();
00155 esTimeflag = (*iter_evt)->getStat();
00156 }
00157 }
00158 if(runFlag ==2) {if( tes>1000 ) return;}
00159 else if(runFlag ==3 ){if (tes>700 ) return;}
00160 else {if (tes>1400 ) return;}
00161
00162 SmartDataPtr<RecMdcDedxCol> newexCol(eventSvc(),"/Event/Recon/RecMdcDedxCol");
00163 if (!newexCol)
00164 {
00165 log << MSG::FATAL << "Could not find RecMdcDedxCol" << endreq;
00166 return;
00167 }
00168
00169 Vint iGood;
00170 iGood.clear();
00171 int nCharge = 0;
00172 double db=0,dz=0,pt0=0,charge0=0;
00173 for(RecMdcDedxCol::iterator it = newexCol->begin(); it != newexCol->end(); it++)
00174 {
00175 HepVector a;
00176 if((*it)->getMdcKalTrack())
00177 {
00178 RecMdcKalTrack* trk = (*it)->getMdcKalTrack();
00179 if(ParticleFlag>-1&&ParticleFlag<5)
00180 {
00181 DstMdcKalTrack* dstTrk = trk;
00182 a = dstTrk->getZHelix(ParticleFlag);
00183 }
00184 else
00185 a = trk->getZHelix();
00186 }
00187 else if((*it)->getMdcTrack())
00188 {
00189 RecMdcTrack* trk = (*it)->getMdcTrack();
00190 a = trk->helix();
00191 }
00192 else continue;
00193 db = a[0];
00194 dz = a[3];
00195 pt0 = fabs(1.0/a[2]);
00196 charge0 = ( a[2] > 0 )? 1 : -1;
00197
00198
00199 if(fabs(dz) >= VZ0CUT) continue;
00200 if(db >= VR0CUT) continue;
00201 if(pt0 >= PT0HighCut) continue;
00202 if(pt0 <= PT0LowCut) continue;
00203 iGood.push_back((*it)->trackId());
00204 nCharge += charge0;
00205 }
00206
00207 if((m_eventType == "isBhabha")||(m_eventType == "isRadBhabha"))
00208 {
00209 if(iGood.size()!=2 || nCharge!=0 ) return;
00210 }
00211
00212
00213
00214 double poca_x=0,poca_y=0,poca_z=0;
00215 float sintheta=0,costheta=0,ptrk=0,ptrk_t=0,charge=0,trackId=0;
00216 int Nhit=0,usedhit=0,Nhits=0,Nphlisthit=0;
00217 double dEdx_meas_hit=0, dEdx_meas=0,dEdx_meas_esat=0,dEdx_meas_norun=0;
00218 double dEdx_hit[100]={0},pathlength_hit[100]={0},wid_hit[100]={0},layid_hit[100]={0},dd_in_hit[100]={0},eangle_hit[100]={0},zhit_hit[100]={0};
00219 int trk_recalg = -1;
00220 Identifier mdcid;
00221
00222 for(RecMdcDedxCol::iterator it = newexCol->begin(); it != newexCol->end(); it++)
00223 {
00224
00225 bool flag = false;
00226 for(int i = 0; i < iGood.size(); i++)
00227 {
00228 if((*it)->trackId()==iGood[i]) flag=true;
00229 }
00230 if(flag==false) continue;
00231
00232 HepVector a;
00233 HepSymMatrix tkErrM;
00234 if((*it)->getMdcKalTrack())
00235 {
00236
00237 poca_x = (*it)->getMdcKalTrack()->x();
00238 poca_y = (*it)->getMdcKalTrack()->y();
00239 poca_z = (*it)->getMdcKalTrack()->z();
00240
00241 RecMdcKalTrack* trk = (*it)->getMdcKalTrack();
00242
00243 if(ParticleFlag>-1&&ParticleFlag<5)
00244 {
00245 DstMdcKalTrack* dstTrk = trk;
00246 a = dstTrk->getFHelix(ParticleFlag);
00247 tkErrM = dstTrk->getFError(ParticleFlag);
00248 }
00249 else
00250 {
00251 a = trk->getFHelix();
00252 tkErrM = trk->getFError();
00253 }
00254 }
00255 else if((*it)->getMdcTrack())
00256 {
00257 poca_x = (*it)->getMdcTrack()->x();
00258 poca_y = (*it)->getMdcTrack()->y();
00259 poca_z = (*it)->getMdcTrack()->z();
00260
00261 RecMdcTrack* trk = (*it)->getMdcTrack();
00262 a = trk->helix();
00263 tkErrM = trk->err();
00264 }
00265 else continue;
00266
00267 sintheta = sin(M_PI_2 - atan(a[4]));
00268 costheta = cos(M_PI_2 - atan(a[4]));
00269 ptrk_t = 1.0/fabs( a[2] );
00270 ptrk = ptrk_t*sqrt(1 + a[4]*a[4]);
00271 charge = ( a[2] > 0 )? 1 : -1;
00272
00273 Nhit = (*it)->numTotalHits();
00274 Nhits = ((*it)->getVecDedxHits()).size();
00275 usedhit = (*it)->numGoodHits();
00276 trk_recalg = (*it)->status();
00277 trackId = (*it)->trackId();
00278
00279 if(m_eventType == "isBhabha")
00280 {
00281 if(runFlag ==3 &&(ptrk>1.88 || ptrk<1.80)) continue;
00282 if(runFlag ==4 &&(ptrk>1.72 || ptrk<1.45)) continue;
00283 if(runFlag ==5 &&(ptrk>2.00 || ptrk<1.70)) continue;
00284 if(runFlag ==6 &&(ptrk>1.90 || ptrk<1.00)) continue;
00285 if(runFlag ==7 &&(ptrk>1.90 || ptrk<1.40)) continue;
00286 if(abs(costheta)>0.9) continue;
00287
00288 if(Nhit<20) continue;
00289 if(usedhit<6) continue;
00290 }
00291
00292
00293 int layid=0,localwid=0,w0id=0,wid=0,lr=0;
00294 double p_hit=0,adc_raw=0,tdc_raw=0,zhit=0,driftd=0,driftD=0,driftT=0,dd_in=0,dd_ex=0,eangle=0;
00295 double ph=0,pathlength=0,rphi_path=0;
00296 long k=0;
00297
00298 DedxHitRefVec gothits = (*it)->getVecDedxHits();
00299 DedxHitRefVec::iterator it_gothit = gothits.begin();
00300 for(;it_gothit!=gothits.end(); it_gothit++)
00301 {
00302
00303 if((*it_gothit)->isMdcHitValid())
00304 {
00305 RecMdcHit* itor = (*it_gothit)->getMdcHit();
00306 mdcid = itor->getMdcId();
00307 layid = MdcID::layer(mdcid);
00308 localwid = MdcID::wire(mdcid);
00309 w0id = geosvc->Layer(layid)->Wirst();
00310 wid = w0id + localwid;
00311 adc_raw = itor->getAdc();
00312 tdc_raw = itor->getTdc();
00313 zhit = itor->getZhit();
00314
00315 lr = itor->getFlagLR();
00316 if(lr == 2) cout<<"lr = "<<lr<<endl;
00317 if(lr == 0 || lr == 2) driftD = itor->getDriftDistLeft();
00318 else driftD = itor->getDriftDistRight();
00319 driftd = abs(driftD);
00320 dd_in = itor->getDoca();
00321 dd_ex = itor->getDoca();
00322 if(lr == 0 || lr == 2 ) {dd_in = -abs(dd_in);dd_ex = -abs(dd_ex);}
00323 if(lr == 1 ) {dd_in = abs(dd_in);dd_ex = abs(dd_ex);}
00324 driftT = itor->getDriftT();
00325 eangle = itor->getEntra();
00326 eangle = eangle/M_PI;
00327 }
00328 else if((*it_gothit)->isMdcKalHelixSegValid())
00329 {
00330 RecMdcKalHelixSeg* itor = (*it_gothit)->getMdcKalHelixSeg();
00331 HepVector a_hit_in = itor->getHelixIncl();
00332 p_hit = 1.0/fabs(a_hit_in(3))*sqrt(1+a_hit_in(5)*a_hit_in(5));
00333
00334 mdcid = itor->getMdcId();
00335 layid = MdcID::layer(mdcid);
00336 localwid = MdcID::wire(mdcid);
00337 w0id = geosvc->Layer(layid)->Wirst();
00338 wid = w0id + localwid;
00339 adc_raw = itor->getAdc();
00340 tdc_raw = itor->getTdc();
00341 zhit = itor->getZhit();
00342
00343 lr = itor->getFlagLR();
00344 if(lr == 2) cout<<"lr = "<<lr<<endl;
00345 driftD = itor->getDD();
00346 driftd = abs(driftD);
00347 driftT = itor->getDT();
00348 dd_in = itor->getDocaIncl();
00349 dd_ex = itor->getDocaExcl();
00350 if(lr==-1 || lr == 2) {dd_in = dd_in; dd_ex = dd_ex;}
00351 else if(lr ==1) {dd_in = -dd_in; dd_ex = -dd_ex;}
00352 eangle = itor->getEntra();
00353 eangle = eangle/M_PI;
00354 }
00355 else continue;
00356
00357 pathlength=(*it_gothit)->getPathLength();
00358 rphi_path=pathlength*sintheta;
00359 ph = (*it_gothit)->getDedx();
00360 if(layid>3)
00361 {
00362 dEdx_hit[k]=adc_raw;
00363 pathlength_hit[k]=pathlength;
00364 wid_hit[k]=wid;
00365 layid_hit[k]=layid;
00366 dd_in_hit[k]=dd_in;
00367 eangle_hit[k]=eangle;
00368 zhit_hit[k]=zhit;
00369
00370 k++;
00371 }
00372
00373
00374 m_charge1 = charge;
00375 m_t01 = tes;
00376 m_driftT = driftT;
00377 m_tdc_raw = tdc_raw;
00378 m_phraw = adc_raw;
00379 m_exraw = ph;
00380 m_localwid = localwid;
00381 m_wire = wid;
00382 m_runNO1 = runNO;
00383 m_runFlag1 = runFlag;
00384 m_doca_in = dd_in;
00385 m_doca_ex = dd_ex;
00386 m_driftdist = driftD;
00387 m_eangle = eangle;
00388 m_zhit = zhit;
00389 m_costheta1 = costheta;
00390 m_pathL = pathlength;
00391 m_layer = layid;
00392 m_ptrk1 = ptrk;
00393 m_ptrk_hit = p_hit;
00394 m_trackId1 = trackId;
00395
00396 StatusCode status = m_nt2->write();
00397 if ( status.isFailure() )
00398 log << MSG::ERROR << "Cannot fill Ntuple n102" << endreq;
00399 }
00400
00401 Nphlisthit = k;
00402 dEdx_meas = (*it)->probPH();
00403
00404 dEdx_meas_esat = (*it)->getDedxEsat();
00405 dEdx_meas_norun = (*it)->getDedxNoRun();
00406 dEdx_meas_hit = (*it)->getDedxHit();
00407
00408
00409
00410 m_poca_x = poca_x;
00411 m_poca_y = poca_y;
00412 m_poca_z = poca_z;
00413 m_ptrk_t=ptrk_t;
00414 m_ptrk=ptrk;
00415 m_sintheta=sintheta;
00416 m_costheta=costheta;
00417 m_charge=charge;
00418 m_runNO = runNO;
00419 m_runFlag = runFlag;
00420 m_evtNO = eventNO;
00421 m_t0 = tes;
00422 m_trackId = trackId;
00423 m_recalg = trk_recalg;
00424
00425 m_nhit=Nhit;
00426 m_nhits=Nhits;
00427 m_nphlisthit=Nphlisthit;
00428 m_usedhit=usedhit;
00429 for(int j=0; j<Nphlisthit; j++)
00430 {
00431 m_dEdx_hit[j]=dEdx_hit[j];
00432 m_pathlength_hit[j]=pathlength_hit[j];
00433 m_wid_hit[j]=wid_hit[j];
00434 m_layid_hit[j]=layid_hit[j];
00435 m_dd_in_hit[j]=dd_in_hit[j];
00436 m_eangle_hit[j]=eangle_hit[j];
00437 m_zhit_hit[j]=zhit_hit[j];
00438 }
00439
00440
00441 m_dEdx_meas = dEdx_meas;
00442
00443
00444
00445 m_parttype = (*it)->particleId();
00446 m_chidedxe=(*it)->chiE();
00447 m_chidedxmu=(*it)->chiMu();
00448 m_chidedxpi=(*it)->chiPi();
00449 m_chidedxk=(*it)->chiK();
00450 m_chidedxp=(*it)->chiP();
00451 for(int i=0;i<5;i++)
00452 {
00453 m_probpid[i]= (*it)->getPidProb(i);
00454 m_expectid[i]= (*it)->getDedxExpect(i);
00455 m_sigmaid[i]= (*it)->getSigmaDedx(i);
00456 }
00457 StatusCode status = m_nt1->write();
00458 if ( status.isFailure() )
00459 {
00460 log << MSG::ERROR << "Cannot fill Ntuple n103" << endreq;
00461 }
00462 }
00463
00464 }