WCSimRecoNtuple.cc

Go to the documentation of this file.
00001 #include "WCSimRecoNtuple.hh"
00002 
00003 #include "WCSimVertexFinder.hh"
00004 
00005 #include "WCSimGeometry.hh"
00006 
00007 #include "WCSimRecoVertex.hh"
00008 #include "WCSimRecoRing.hh"
00009 
00010 #include "TDirectory.h"
00011 #include "TFile.h"
00012 #include "TTree.h"
00013 #include "TMath.h"
00014 
00015 #include <iostream>
00016 #include <cmath>
00017 
00018 ClassImp(WCSimRecoNtuple)
00019 
00020 WCSimRecoNtuple::WCSimRecoNtuple() : WCSimNtuple(),
00021   fRecoFile(0),
00022   fRecoTree(0)
00023 {
00024   
00025 }
00026   
00027 WCSimRecoNtuple::~WCSimRecoNtuple()
00028 {
00029   
00030 }
00031 
00032 void WCSimRecoNtuple::Reset() 
00033 { 
00034   std::cout << " *** WCSimRecoNtuple::Reset() *** " << std::endl;
00035 
00036   fRecoFile = 0;
00037   fRecoTree = 0;
00038 
00039   return;
00040 }
00041 
00042 void WCSimRecoNtuple::Fill( WCSimTrueEvent* trueEvent, WCSimRecoEvent* recoEvent )
00043 {
00044   this->ResetVariables();
00045   this->WriteVariables( trueEvent, recoEvent );
00046 }
00047   
00048 void WCSimRecoNtuple::ResetVariables()
00049 {
00050   // reset variables
00051   // ===============
00052   fRunNum = -1;
00053   fEventNum = -1;
00054   fTriggerNum = -1;
00055 
00056   fTruePDG = 0;
00057   fTrueP = 0.0;
00058   fTrueE = 0.0; 
00059   fTrueG4VtxX = 0.0;
00060   fTrueG4VtxY = 0.0;
00061   fTrueG4VtxZ = 0.0;
00062   fTrueVtxX = 0.0;
00063   fTrueVtxY = 0.0;
00064   fTrueVtxZ = 0.0;
00065   fTrueDirX = 0.0;
00066   fTrueDirY = 0.0;
00067   fTrueDirZ = 0.0;
00068   fTrueVtxTraceFwd = -99999.9;
00069   fTrueVtxTraceBwd = -99999.9;
00070   fTrueVtxDistToEdge = -99999.9;
00071 
00072   fTrueG4EndX = 0.0;
00073   fTrueG4EndY = 0.0;
00074   fTrueG4EndZ = 0.0;
00075   fTrueEndX = 0.0;
00076   fTrueEndY = 0.0;
00077   fTrueEndZ = 0.0;
00078   fTrueEndDistToEdge = -99999.9;
00079   fTrueLength = 0.0;
00080 
00081   fNDigits = 0;
00082   fNFilterDigits = 0;
00083   fNRings = 0;
00084 
00085   fFoundVertex = 0;
00086   fFoundDirection = 0;
00087   fFoundRings = 0;
00088 
00089   fRecoVtxX = 0.0;
00090   fRecoVtxY = 0.0;
00091   fRecoVtxZ = 0.0;
00092   fRecoVtxTime = 0.0;
00093   fRecoDirX = 0.0;
00094   fRecoDirY = 0.0;
00095   fRecoDirZ = 0.0;
00096   fRecoConeAngle = 0.0;
00097   fRecoTrackLength = 0.0;
00098   fRecoVtxFOM = 0.0;
00099   fRecoVtxSteps = 0;
00100   fRecoVtxPass = 0;
00101   fRecoVtxStatus = 0;
00102 
00103   fRecoVtxTraceFwd = -99999.9;
00104   fRecoVtxTraceBwd = -99999.9;
00105   fRecoVtxDistToEdge = -99999.9;
00106 
00107   fRecoDelta = -99999.9;
00108   fRecoDeltaTrans = -99999.9;
00109   fRecoDeltaLong = -99999.9;
00110   fRecoDeltaAngle = -99999.9;
00111 
00112   fRecoRingVtxX = 0.0;
00113   fRecoRingVtxY = 0.0;
00114   fRecoRingVtxZ = 0.0;  
00115   fRecoRingDirX = 0.0;
00116   fRecoRingDirY = 0.0;
00117   fRecoRingDirZ = 0.0;
00118   fRecoRingAngle = 0.0;
00119   fRecoRingHeight = 0.0;
00120   fRecoRingHeightFrac = 0.0;
00121 
00122   fRecoRingDeltaAngle = -99999.9;
00123 
00124   // vertex reconstruction
00125   // =====================
00126   fSimplePositionX = 0.0;
00127   fSimplePositionY = 0.0;
00128   fSimplePositionZ = 0.0;
00129   fSimplePositionDelta = -99999.9;
00130   fSimplePositionDeltaTrans = -99999.9;
00131   fSimplePositionDeltaLong = -99999.9;
00132 
00133   fPointPositionX = 0.0;
00134   fPointPositionY = 0.0;
00135   fPointPositionZ = 0.0;
00136   fPointPositionFOM = 0.0;
00137   fPointPositionSteps = 0;
00138   fPointPositionPass = 0;
00139   fPointPositionDelta = -99999.9;
00140   fPointPositionDeltaTrans = -99999.9;
00141   fPointPositionDeltaLong = -99999.9;
00142 
00143   fSimpleDirectionX = 0.0;
00144   fSimpleDirectionY = 0.0;
00145   fSimpleDirectionZ = 0.0;
00146   fSimpleDirectionDeltaAngle = -99999.9;
00147 
00148   fPointDirectionX = 0.0;
00149   fPointDirectionY = 0.0;
00150   fPointDirectionZ = 0.0;
00151   fPointDirectionFOM = 0.0;
00152   fPointDirectionSteps = 0;
00153   fPointDirectionPass = 0;  
00154   fPointDirectionDeltaAngle = -99999.9;
00155 
00156   fPointVtxX = 0.0;
00157   fPointVtxY = 0.0;
00158   fPointVtxZ = 0.0;
00159   fPointDirX = 0.0;
00160   fPointDirY = 0.0;
00161   fPointDirZ = 0.0;
00162   fPointConeAngle = 0.0;
00163   fPointTrackLength = 0.0;
00164   fPointVtxFOM = 0.0;
00165   fPointVtxSteps = 0;
00166   fPointVtxPass = 0;  
00167   fPointVtxDelta = -99999.9;
00168   fPointVtxDeltaTrans = -99999.9;
00169   fPointVtxDeltaLong = -99999.9;
00170   fPointVtxDeltaAngle = -99999.9;
00171 
00172   fExtendedVtxX = 0.0;
00173   fExtendedVtxY = 0.0;
00174   fExtendedVtxZ = 0.0;
00175   fExtendedDirX = 0.0;
00176   fExtendedDirY = 0.0;
00177   fExtendedDirZ = 0.0;  
00178   fExtendedConeAngle = 0.0;
00179   fExtendedTrackLength = 0.0;
00180   fExtendedVtxFOM = 0.0;
00181   fExtendedVtxSteps = 0;
00182   fExtendedVtxPass = 0;  
00183   fExtendedVtxDelta = -99999.9;
00184   fExtendedVtxDeltaTrans = -99999.9;
00185   fExtendedVtxDeltaLong = -99999.9;
00186   fExtendedVtxDeltaAngle = -99999.9;
00187 
00188   // --- debug ---
00189   fExtendedVtxTimeParam0 = 0.0;
00190   fExtendedVtxConeParam0 = 0.0;
00191   fExtendedVtxConeParam1 = 0.0;
00192   fExtendedVtxConeParam2 = 0.0;
00193 
00194   return;
00195 }
00196 
00197 void WCSimRecoNtuple::WriteVariables( WCSimTrueEvent* fTrueEvent, WCSimRecoEvent* fRecoEvent )
00198 {
00199   std::cout << " *** WCSimRecoNtuple::WriteVariables() *** " << std::endl;
00200 
00201   // header information
00202   // ==================
00203   fRunNum = fRecoEvent->GetRun();
00204   fEventNum = fRecoEvent->GetEvent();
00205   fTriggerNum = fRecoEvent->GetTrigger();
00206 
00207   // truth information
00208   // =================
00209   fTruePDG  = fTrueEvent->GetPDG();
00210   fTrueP     = fTrueEvent->GetMomentum();
00211   fTrueE     = fTrueEvent->GetEnergy();  
00212 
00213   fTrueG4VtxX  = fTrueEvent->GetG4VtxX();
00214   fTrueG4VtxY  = fTrueEvent->GetG4VtxY();
00215   fTrueG4VtxZ  = fTrueEvent->GetG4VtxZ();
00216   fTrueG4EndX  = fTrueEvent->GetG4EndX();
00217   fTrueG4EndY  = fTrueEvent->GetG4EndY();
00218   fTrueG4EndZ  = fTrueEvent->GetG4EndZ();
00219 
00220   fTrueVtxX  = fTrueEvent->GetVtxX();
00221   fTrueVtxY  = fTrueEvent->GetVtxY();
00222   fTrueVtxZ  = fTrueEvent->GetVtxZ();
00223   fTrueEndX  = fTrueEvent->GetEndX();
00224   fTrueEndY  = fTrueEvent->GetEndY();
00225   fTrueEndZ  = fTrueEvent->GetEndZ();
00226   fTrueDirX  = fTrueEvent->GetDirX();
00227   fTrueDirY  = fTrueEvent->GetDirY();
00228   fTrueDirZ  = fTrueEvent->GetDirZ();
00229 
00230   fTrueLength = fTrueEvent->GetLength();
00231   
00232   if( fTruePDG!=0 ){
00233     fTrueVtxTraceFwd = WCSimGeometry::Instance()->ForwardProjectionToEdge(fTrueVtxX,fTrueVtxY,fTrueVtxZ,
00234                                                                           fTrueDirX,fTrueDirY,fTrueDirZ);
00235     fTrueVtxTraceBwd = WCSimGeometry::Instance()->BackwardProjectionToEdge(fTrueVtxX,fTrueVtxY,fTrueVtxZ,
00236                                                                            fTrueDirX,fTrueDirY,fTrueDirZ);
00237     fTrueVtxDistToEdge = WCSimGeometry::Instance()->DistanceToEdge(fTrueVtxX,fTrueVtxY,fTrueVtxZ);
00238 
00239     fTrueEndDistToEdge = WCSimGeometry::Instance()->DistanceToEdge(fTrueEndX,fTrueEndY,fTrueEndZ);
00240   }
00241 
00242   // reconstructed variables
00243   // =======================
00244   fNDigits = fRecoEvent->GetNDigits();
00245   fNFilterDigits = fRecoEvent->GetNFilterDigits();
00246   fNRings = fRecoEvent->GetNRings();
00247 
00248   fFoundVertex = fRecoEvent->FoundVertex();
00249   fRecoVtxX = fRecoEvent->GetVtxX();
00250   fRecoVtxY = fRecoEvent->GetVtxY();
00251   fRecoVtxZ = fRecoEvent->GetVtxZ();
00252   fRecoVtxTime = fRecoEvent->GetVtxTime();
00253   
00254   fFoundDirection = fRecoEvent->FoundDirection();
00255   fRecoDirX = fRecoEvent->GetDirX();
00256   fRecoDirY = fRecoEvent->GetDirY();
00257   fRecoDirZ = fRecoEvent->GetDirZ();
00258 
00259   fRecoConeAngle = fRecoEvent->GetConeAngle();
00260   fRecoTrackLength = fRecoEvent->GetTrackLength();
00261 
00262   fRecoVtxFOM = fRecoEvent->GetVtxFOM(); 
00263   fRecoVtxSteps = fRecoEvent->GetVtxIterations(); 
00264   fRecoVtxPass = fRecoEvent->GetVtxPass();
00265   fRecoVtxStatus = fRecoEvent->GetVtxStatus();
00266 
00267   if( fFoundVertex && fFoundDirection ){
00268     fRecoVtxTraceFwd = WCSimGeometry::Instance()->ForwardProjectionToEdge(fRecoVtxX,fRecoVtxY,fRecoVtxZ,
00269                                                                           fRecoDirX,fRecoDirY,fRecoDirZ);
00270     fRecoVtxTraceBwd = WCSimGeometry::Instance()->BackwardProjectionToEdge(fRecoVtxX,fRecoVtxY,fRecoVtxZ,
00271                                                                            fRecoDirX,fRecoDirY,fRecoDirZ);
00272     fRecoVtxDistToEdge = WCSimGeometry::Instance()->DistanceToEdge(fRecoVtxX,fRecoVtxY,fRecoVtxZ);
00273   }
00274 
00275   fFoundRings = fRecoEvent->FoundRings();
00276   
00277   if( fFoundRings ){
00278     WCSimRecoRing* myPrimaryRing = (WCSimRecoRing*)(fRecoEvent->GetPrimaryRing());
00279     fRecoRingVtxX = myPrimaryRing->GetVtxX();
00280     fRecoRingVtxY = myPrimaryRing->GetVtxY();
00281     fRecoRingVtxZ = myPrimaryRing->GetVtxZ();  
00282     fRecoRingDirX = myPrimaryRing->GetDirX();
00283     fRecoRingDirY = myPrimaryRing->GetDirY();
00284     fRecoRingDirZ = myPrimaryRing->GetDirZ();
00285     fRecoRingAngle = myPrimaryRing->GetAngle();
00286     fRecoRingHeight =  myPrimaryRing->GetHeight();
00287     fRecoRingHeightFrac = 0.0;
00288 
00289     if( fRecoRingHeight>0.0 ){
00290       fRecoRingHeightFrac = fRecoRingHeight/(double)fNFilterDigits;
00291     }
00292   }
00293 
00294   // compare with truth
00295   // ==================
00296   if( fTruePDG!=0 
00297    && fFoundVertex && fFoundDirection ){
00298     Double_t dx = fRecoVtxX - fTrueVtxX;
00299     Double_t dy = fRecoVtxY - fTrueVtxY;
00300     Double_t dz = fRecoVtxZ - fTrueVtxZ;
00301     Double_t ds = sqrt( dx*dx + dy*dy + dz*dz );
00302     Double_t px = fTrueDirX;
00303     Double_t py = fTrueDirY;
00304     Double_t pz = fTrueDirZ;
00305 
00306     Double_t epsilon = 1.0e-7;
00307 
00308     if( ds>epsilon ){
00309       px = dx/ds;
00310       py = dy/ds;
00311       pz = dz/ds;
00312     }
00313 
00314     Double_t costh = px*fTrueDirX+py*fTrueDirY+pz*fTrueDirZ; 
00315     Double_t sinth = 0.0;
00316     if( costh<1.0-epsilon ) sinth = sqrt(1.0-costh*costh);
00317     else costh = 1.0;
00318 
00319     Double_t costh2 = fTrueDirX*fRecoDirX + fTrueDirY*fRecoDirY + fTrueDirZ*fRecoDirZ;
00320     Double_t angle = 0.0;
00321     if( costh2<1.0-epsilon ) angle = (180.0/TMath::Pi())*acos(costh2);
00322     else angle = 0.0;
00323 
00324     fRecoDelta = ds;
00325     fRecoDeltaTrans = ds*sinth;
00326     fRecoDeltaLong = ds*costh;
00327     fRecoDeltaAngle = angle;
00328   }
00329 
00330   if( fTruePDG!=0 
00331    && fFoundRings ){
00332     Double_t epsilon = 1.0e-7;
00333     Double_t costh2 = fTrueDirX*fRecoRingDirX+fTrueDirY*fRecoRingDirY+fTrueDirZ*fRecoRingDirZ;  
00334     Double_t angle = 0.0;
00335     if( costh2<1.0-epsilon ) angle = (180.0/TMath::Pi())*acos(costh2);
00336     else angle = 0.0;
00337 
00338     fRecoRingDeltaAngle = angle;  
00339   }
00340 
00341   // reconstruction chain
00342   // =====================
00343   if( fFoundVertex && fFoundDirection ){
00344     Double_t dx = 0.0;
00345     Double_t dy = 0.0;
00346     Double_t dz = 0.0;
00347     Double_t ds = 0.0;
00348     Double_t px = 0.0;
00349     Double_t py = 0.0;
00350     Double_t pz = 0.0;  
00351     Double_t qx = 0.0;
00352     Double_t qy = 0.0;
00353     Double_t qz = 0.0;
00354 
00355     Double_t epsilon = 1.0e-7;
00356 
00357     Double_t costh = 0.0;
00358     Double_t sinth = 0.0;
00359     Double_t costh2 = 0.0;
00360     Double_t angle = 0.0;
00361 
00362     // simple position
00363     WCSimRecoVertex* simplePos = WCSimVertexFinder::Instance()->GetSimplePosition();
00364     fSimplePositionX = simplePos->GetX();
00365     fSimplePositionY = simplePos->GetY();
00366     fSimplePositionZ = simplePos->GetZ();  
00367 
00368     dx = fSimplePositionX - fTrueVtxX;
00369     dy = fSimplePositionY - fTrueVtxY;
00370     dz = fSimplePositionZ - fTrueVtxZ;
00371     ds = sqrt( dx*dx + dy*dy + dz*dz );
00372     px = fTrueDirX;
00373     py = fTrueDirY;
00374     pz = fTrueDirZ;
00375 
00376     if( ds>epsilon ){
00377       px = dx/ds;
00378       py = dy/ds;
00379       pz = dz/ds;
00380     }
00381 
00382     costh = px*fTrueDirX+py*fTrueDirY+pz*fTrueDirZ; 
00383     sinth = 0.0;
00384     if( costh<1.0-epsilon ) sinth = sqrt(1.0-costh*costh);
00385     else costh = 1.0;
00386 
00387     fSimplePositionDelta = ds;
00388     fSimplePositionDeltaTrans = ds*sinth;
00389     fSimplePositionDeltaLong = ds*costh;
00390 
00391 
00392     // point position
00393     WCSimRecoVertex* pointPos = WCSimVertexFinder::Instance()->GetPointPosition();
00394     fPointPositionX = pointPos->GetX();
00395     fPointPositionY = pointPos->GetY();
00396     fPointPositionZ = pointPos->GetZ();
00397     fPointPositionFOM = pointPos->GetFOM();
00398     fPointPositionSteps = pointPos->GetIterations();
00399     fPointPositionPass = pointPos->GetPass();
00400  
00401     dx = fPointPositionX - fTrueVtxX;
00402     dy = fPointPositionY - fTrueVtxY;
00403     dz = fPointPositionZ - fTrueVtxZ;
00404     ds = sqrt( dx*dx + dy*dy + dz*dz );
00405     px = fTrueDirX;
00406     py = fTrueDirY;
00407     pz = fTrueDirZ;
00408 
00409     if( ds>epsilon ){
00410       px = dx/ds;
00411       py = dy/ds;
00412       pz = dz/ds;
00413     }
00414 
00415     costh = px*fTrueDirX+py*fTrueDirY+pz*fTrueDirZ; 
00416     sinth = 0.0;
00417     if( costh<1.0-epsilon ) sinth = sqrt(1.0-costh*costh);
00418     else costh = 1.0;
00419 
00420     fPointPositionDelta = ds;
00421     fPointPositionDeltaTrans = ds*sinth;
00422     fPointPositionDeltaLong = ds*costh;
00423 
00424     // simple direction
00425     WCSimRecoVertex* simpleDir = WCSimVertexFinder::Instance()->GetSimpleDirection();
00426     fSimpleDirectionX = simpleDir->GetDirX();
00427     fSimpleDirectionY = simpleDir->GetDirY();
00428     fSimpleDirectionZ = simpleDir->GetDirZ();
00429  
00430     qx = fSimpleDirectionX;
00431     qy = fSimpleDirectionY;
00432     qz = fSimpleDirectionZ;
00433  
00434     costh2 = qx*fTrueDirX + qy*fTrueDirY + qz*fTrueDirZ;
00435     if( costh2<1.0-epsilon ) angle = (180.0/TMath::Pi())*acos(costh2);
00436     else angle = 0.0;
00437 
00438     fSimpleDirectionDeltaAngle = angle;
00439 
00440     // point direction
00441     WCSimRecoVertex* pointDir = WCSimVertexFinder::Instance()->GetPointDirection();
00442     fPointDirectionX = pointDir->GetDirX();
00443     fPointDirectionY = pointDir->GetDirY();
00444     fPointDirectionZ = pointDir->GetDirZ();
00445     fPointDirectionFOM = pointDir->GetFOM();
00446     fPointDirectionSteps = pointDir->GetIterations();
00447     fPointDirectionPass = pointDir->GetPass();
00448 
00449     qx = fPointDirectionX;
00450     qy = fPointDirectionY;
00451     qz = fPointDirectionZ;
00452 
00453     costh2 = qx*fTrueDirX + qy*fTrueDirY + qz*fTrueDirZ;
00454     if( costh2<1.0-epsilon ) angle = (180.0/TMath::Pi())*acos(costh2);
00455     else angle = 0.0;
00456 
00457     fPointDirectionDeltaAngle = angle;
00458 
00459     // point vertex
00460     WCSimRecoVertex* pointVtx = WCSimVertexFinder::Instance()->GetPointVertex();
00461     fPointVtxX = pointVtx->GetX();
00462     fPointVtxY = pointVtx->GetY();
00463     fPointVtxZ = pointVtx->GetZ(); 
00464     fPointDirX = pointVtx->GetDirX();
00465     fPointDirY = pointVtx->GetDirY();
00466     fPointDirZ = pointVtx->GetDirZ();
00467     fPointConeAngle = pointVtx->GetConeAngle();
00468     fPointTrackLength = pointVtx->GetTrackLength();
00469     fPointVtxFOM = pointVtx->GetFOM();
00470     fPointVtxSteps = pointVtx->GetIterations();
00471     fPointVtxPass = pointVtx->GetPass();
00472 
00473     dx = fPointVtxX - fTrueVtxX;
00474     dy = fPointVtxY - fTrueVtxY;
00475     dz = fPointVtxZ - fTrueVtxZ;
00476     ds = sqrt( dx*dx + dy*dy + dz*dz );
00477     qx = fPointDirX;
00478     qy = fPointDirY;
00479     qz = fPointDirZ;
00480     px = fTrueDirX;
00481     py = fTrueDirY;
00482     pz = fTrueDirZ;
00483 
00484     if( ds>epsilon ){
00485       px = dx/ds;
00486       py = dy/ds;
00487       pz = dz/ds;
00488     }
00489 
00490     costh = px*fTrueDirX+py*fTrueDirY+pz*fTrueDirZ; 
00491     sinth = 0.0;
00492     if( costh<1.0-epsilon ) sinth = sqrt(1.0-costh*costh);
00493     else costh = 1.0;
00494 
00495     costh2 = qx*fTrueDirX + qy*fTrueDirY + qz*fTrueDirZ;
00496     if( costh2<1.0-epsilon ) angle = (180.0/TMath::Pi())*acos(costh2);
00497     else angle = 0.0;
00498 
00499     fPointVtxDelta = ds;
00500     fPointVtxDeltaTrans = ds*sinth;
00501     fPointVtxDeltaLong = ds*costh;
00502     fPointVtxDeltaAngle = angle;
00503 
00504     // extended
00505     WCSimRecoVertex* extendedVtx = WCSimVertexFinder::Instance()->GetExtendedVertex();
00506     fExtendedVtxX = extendedVtx->GetX();
00507     fExtendedVtxY = extendedVtx->GetY();
00508     fExtendedVtxZ = extendedVtx->GetZ();
00509     fExtendedDirX = extendedVtx->GetDirX();
00510     fExtendedDirY = extendedVtx->GetDirY();
00511     fExtendedDirZ = extendedVtx->GetDirZ();   
00512     fExtendedConeAngle = extendedVtx->GetConeAngle();
00513     fExtendedTrackLength = extendedVtx->GetTrackLength();
00514     fExtendedVtxFOM = extendedVtx->GetFOM();
00515     fExtendedVtxSteps = extendedVtx->GetIterations();
00516     fExtendedVtxPass = extendedVtx->GetPass();   
00517 
00518     dx = fExtendedVtxX - fTrueVtxX;
00519     dy = fExtendedVtxY - fTrueVtxY;
00520     dz = fExtendedVtxZ - fTrueVtxZ;
00521     ds = sqrt( dx*dx + dy*dy + dz*dz );
00522     qx = fExtendedDirX;
00523     qy = fExtendedDirY;
00524     qz = fExtendedDirZ;
00525     px = fTrueDirX;
00526     py = fTrueDirY;
00527     pz = fTrueDirZ;
00528 
00529     if( ds>epsilon ){
00530       px = dx/ds;
00531       py = dy/ds;
00532       pz = dz/ds;
00533     }
00534 
00535     costh = px*fTrueDirX+py*fTrueDirY+pz*fTrueDirZ; 
00536     sinth = 0.0;
00537     if( costh<1.0-epsilon ) sinth = sqrt(1.0-costh*costh);
00538     else costh = 1.0;
00539 
00540     costh2 = qx*fTrueDirX + qy*fTrueDirY + qz*fTrueDirZ;
00541     if( costh2<1.0-epsilon ) angle = (180.0/TMath::Pi())*acos(costh2);
00542     else angle = 0.0;
00543 
00544     fExtendedVtxDelta = ds;
00545     fExtendedVtxDeltaTrans = ds*sinth;
00546     fExtendedVtxDeltaLong = ds*costh;
00547     fExtendedVtxDeltaAngle = angle;
00548 
00549     // --- debug ---
00550     fExtendedVtxTimeParam0 = WCSimVertexFinder::Instance()->fTimeParam0;
00551     fExtendedVtxConeParam0 = WCSimVertexFinder::Instance()->fConeParam0;
00552     fExtendedVtxConeParam1 = WCSimVertexFinder::Instance()->fConeParam1;
00553     fExtendedVtxConeParam2 = WCSimVertexFinder::Instance()->fConeParam2;
00554   }
00555 
00556   // write variables
00557   // ===============
00558   this->WriteToFile();
00559 
00560   return;
00561 }
00562 
00563 void WCSimRecoNtuple::WriteToFile()
00564 {
00565   TDirectory* tmpd = 0;
00566 
00567   if( fRecoFile==0 ){
00568     tmpd = gDirectory;
00569     std::cout << " *** WCSimRecoNtuple::OpenFile() *** " << std::endl;
00570     std::cout << "  opening file: " << fNtpFileName.Data() << std::endl;
00571     fRecoFile = new TFile(fNtpFileName.Data(),"recreate");
00572     fRecoTree = new TTree("ntuple","my analysis ntuple");
00573 
00574     fRecoTree->Branch("run",&fRunNum,"run/I");
00575     fRecoTree->Branch("event",&fEventNum,"event/I");
00576     fRecoTree->Branch("trigger",&fTriggerNum,"trigger/I");
00577 
00578     fRecoTree->Branch("truePDG",&fTruePDG,"truePDG/I");
00579     fRecoTree->Branch("trueP",&fTrueP,"trueP/D");
00580     fRecoTree->Branch("trueE",&fTrueE,"trueE/D");  
00581 
00582     fRecoTree->Branch("trueG4VtxX",&fTrueG4VtxX,"trueG4VtxX/D");
00583     fRecoTree->Branch("trueG4VtxY",&fTrueG4VtxY,"trueG4VtxY/D");
00584     fRecoTree->Branch("trueG4VtxZ",&fTrueG4VtxZ,"trueG4VtxZ/D");  
00585     fRecoTree->Branch("trueG4EndX",&fTrueG4EndX,"trueG4EndX/D");
00586     fRecoTree->Branch("trueG4EndY",&fTrueG4EndY,"trueG4EndY/D");
00587     fRecoTree->Branch("trueG4EndZ",&fTrueG4EndZ,"trueG4EndZ/D");  
00588 
00589     fRecoTree->Branch("trueVtxX",&fTrueVtxX,"trueVtxX/D");
00590     fRecoTree->Branch("trueVtxY",&fTrueVtxY,"trueVtxY/D");
00591     fRecoTree->Branch("trueVtxZ",&fTrueVtxZ,"trueVtxZ/D");   
00592     fRecoTree->Branch("trueEndX",&fTrueEndX,"trueEndX/D");
00593     fRecoTree->Branch("trueEndY",&fTrueEndY,"trueEndY/D");
00594     fRecoTree->Branch("trueEndZ",&fTrueEndZ,"trueEndZ/D"); 
00595     fRecoTree->Branch("trueDirX",&fTrueDirX,"trueDirX/D");
00596     fRecoTree->Branch("trueDirY",&fTrueDirY,"trueDirY/D");
00597     fRecoTree->Branch("trueDirZ",&fTrueDirZ,"trueDirZ/D");
00598     
00599     fRecoTree->Branch("trueLength",&fTrueLength,"trueLength/D");
00600 
00601     fRecoTree->Branch("trueVtxTraceFwd",&fTrueVtxTraceFwd,"trueVtxTraceFwd/D");
00602     fRecoTree->Branch("trueVtxTraceBwd",&fTrueVtxTraceBwd,"trueVtxTraceBwd/D");
00603     fRecoTree->Branch("trueVtxDistToEdge",&fTrueVtxDistToEdge,"trueVtxDistToEdge/D");
00604     fRecoTree->Branch("trueEndDistToEdge",&fTrueEndDistToEdge,"trueEndDistToEdge/D");
00605 
00606     fRecoTree->Branch("nDigits",&fNDigits,"nDigits/I");
00607     fRecoTree->Branch("nFilterDigits",&fNFilterDigits,"nFilterDigits/I");
00608     fRecoTree->Branch("nRings",&fNRings,"nRings/I");
00609     fRecoTree->Branch("foundVertex",&fFoundVertex,"foundVertex/I");
00610     fRecoTree->Branch("foundDirection",&fFoundDirection,"foundDirection/I");
00611     fRecoTree->Branch("foundRings",&fFoundRings,"foundRings/I");
00612 
00613     fRecoTree->Branch("recoVtxX",&fRecoVtxX,"recoVtxX/D");
00614     fRecoTree->Branch("recoVtxY",&fRecoVtxY,"recoVtxY/D");
00615     fRecoTree->Branch("recoVtxZ",&fRecoVtxZ,"recoVtxZ/D");
00616     fRecoTree->Branch("recoVtxTime",&fRecoVtxTime,"recoVtxTime/D");
00617     fRecoTree->Branch("recoDirX",&fRecoDirX,"recoDirX/D");
00618     fRecoTree->Branch("recoDirY",&fRecoDirY,"recoDirY/D");
00619     fRecoTree->Branch("recoDirZ",&fRecoDirZ,"recoDirZ/D");
00620     fRecoTree->Branch("recoConeAngle",&fRecoConeAngle,"recoConeAngle/D");
00621     fRecoTree->Branch("recoTrackLength",&fRecoTrackLength,"recoTrackLength/D");
00622     fRecoTree->Branch("recoVtxFOM",&fRecoVtxFOM,"recoVtxFOM/D");
00623     fRecoTree->Branch("recoVtxSteps",&fRecoVtxSteps,"recoVtxSteps/I");  
00624     fRecoTree->Branch("recoVtxPass",&fRecoVtxPass,"recoVtxPass/I");
00625     fRecoTree->Branch("recoVtxStatus",&fRecoVtxStatus,"recoVtxStatus/I");
00626     fRecoTree->Branch("recoVtxTraceFwd",&fRecoVtxTraceFwd,"recoVtxTraceFwd/D");
00627     fRecoTree->Branch("recoVtxTraceBwd",&fRecoVtxTraceBwd,"recoVtxTraceBwd/D");
00628     fRecoTree->Branch("recoVtxDistToEdge",&fRecoVtxDistToEdge,"recoVtxDistToEdge/D");
00629 
00630     fRecoTree->Branch("recoRingVtxX",&fRecoRingVtxX,"recoRingVtxX/D");
00631     fRecoTree->Branch("recoRingVtxY",&fRecoRingVtxY,"recoRingVtxY/D");
00632     fRecoTree->Branch("recoRingVtxZ",&fRecoRingVtxZ,"recoRingVtxZ/D");  
00633     fRecoTree->Branch("recoRingDirX",&fRecoRingDirX,"recoRingDirX/D");
00634     fRecoTree->Branch("recoRingDirY",&fRecoRingDirY,"recoRingDirY/D");
00635     fRecoTree->Branch("recoRingDirZ",&fRecoRingDirZ,"recoRingDirZ/D");
00636     fRecoTree->Branch("recoRingAngle",&fRecoRingAngle,"recoRingAngle/D");
00637     fRecoTree->Branch("recoRingHeight",&fRecoRingHeight,"recoRingHeight/D");
00638     fRecoTree->Branch("recoRingHeightFrac",&fRecoRingHeightFrac,"recoRingHeightFrac/D");
00639 
00640     fRecoTree->Branch("recoDelta",&fRecoDelta,"recoDelta/D");
00641     fRecoTree->Branch("recoDeltaTrans",&fRecoDeltaTrans,"recoDeltaTrans/D");
00642     fRecoTree->Branch("recoDeltaLong",&fRecoDeltaLong,"recoDeltaLong/D");
00643     fRecoTree->Branch("recoDeltaAngle",&fRecoDeltaAngle,"recoDeltaAngle/D");
00644     fRecoTree->Branch("recoRingDeltaAngle",&fRecoRingDeltaAngle,"recoRingDeltaAngle/D");
00645 
00646     fRecoTree->Branch("simplePositionX",&fSimplePositionX,"simplePositionX/D");
00647     fRecoTree->Branch("simplePositionY",&fSimplePositionY,"simplePositionY/D");
00648     fRecoTree->Branch("simplePositionZ",&fSimplePositionZ,"simplePositionZ/D");
00649     fRecoTree->Branch("simplePositionDelta",&fSimplePositionDelta,"simplePositionDelta/D");
00650     fRecoTree->Branch("simplePositionDeltaTrans",&fSimplePositionDeltaTrans,"simplePositionDeltaTrans/D");
00651     fRecoTree->Branch("simplePositionDeltaLong",&fSimplePositionDeltaLong,"simplePositionDeltaLong/D");
00652 
00653     fRecoTree->Branch("pointPositionX",&fPointPositionX,"pointPositionX/D");
00654     fRecoTree->Branch("pointPositionY",&fPointPositionY,"pointPositionY/D");
00655     fRecoTree->Branch("pointPositionZ",&fPointPositionZ,"pointPositionZ/D");
00656     fRecoTree->Branch("pointPositionFOM",&fPointPositionFOM,"pointPositionFOM/D");  
00657     fRecoTree->Branch("pointPositionSteps",&fPointPositionSteps,"pointPositionSteps/I");
00658     fRecoTree->Branch("pointPositionPass",&fPointPositionPass,"pointPositionPass/I");
00659     fRecoTree->Branch("pointPositionDelta",&fPointPositionDelta,"pointPositionDelta/D");
00660     fRecoTree->Branch("pointPositionDeltaTrans",&fPointPositionDeltaTrans,"pointPositionDeltaTrans/D");
00661     fRecoTree->Branch("pointPositionDeltaLong",&fPointPositionDeltaLong,"pointPositionDeltaLong/D");
00662 
00663     fRecoTree->Branch("simpleDirectionX",&fSimpleDirectionX,"simpleDirectionX/D");
00664     fRecoTree->Branch("simpleDirectionY",&fSimpleDirectionY,"simpleDirectionY/D");
00665     fRecoTree->Branch("simpleDirectionZ",&fSimpleDirectionZ,"simpleDirectionZ/D");
00666     fRecoTree->Branch("simpleDirectionDeltaAngle",&fSimpleDirectionDeltaAngle,"simpleDirectionDeltaAngle/D");
00667 
00668     fRecoTree->Branch("pointDirectionX",&fPointDirectionX,"pointDirectionX/D");
00669     fRecoTree->Branch("pointDirectionY",&fPointDirectionY,"pointDirectionY/D");
00670     fRecoTree->Branch("pointDirectionZ",&fPointDirectionZ,"pointDirectionZ/D");
00671     fRecoTree->Branch("pointDirectionFOM",&fPointDirectionFOM,"pointDirectionFOM/D");
00672     fRecoTree->Branch("pointDirectionSteps",&fPointDirectionSteps,"pointDirSteps/I");
00673     fRecoTree->Branch("pointDirectionPass",&fPointDirectionPass,"pointDirPass/I");
00674     fRecoTree->Branch("pointDirectionDeltaAngle",&fPointDirectionDeltaAngle,"pointDirectionDeltaAngle/D");
00675 
00676     fRecoTree->Branch("pointVtxX",&fPointVtxX,"pointVtxX/D");
00677     fRecoTree->Branch("pointVtxY",&fPointVtxY,"pointVtxY/D");
00678     fRecoTree->Branch("pointVtxZ",&fPointVtxZ,"pointVtxZ/D");
00679     fRecoTree->Branch("pointDirX",&fPointDirX,"pointDirX/D");
00680     fRecoTree->Branch("pointDirY",&fPointDirY,"pointDirY/D");
00681     fRecoTree->Branch("pointDirZ",&fPointDirZ,"pointDirZ/D");
00682     fRecoTree->Branch("pointConeAngle",&fPointConeAngle,"pointConeAngle/D");
00683     fRecoTree->Branch("pointTrackLength",&fPointTrackLength,"pointTrackLength/D");
00684     fRecoTree->Branch("pointVtxFOM",&fPointVtxFOM,"pointVtxFOM/D");
00685     fRecoTree->Branch("pointVtxSteps",&fPointVtxSteps,"pointVtxSteps/I");
00686     fRecoTree->Branch("pointVtxPass",&fPointVtxPass,"pointVtxPass/I");
00687     fRecoTree->Branch("pointVtxDelta",&fPointVtxDelta,"pointVtxDelta/D");
00688     fRecoTree->Branch("pointVtxDeltaTrans",&fPointVtxDeltaTrans,"pointVtxDeltaTrans/D");
00689     fRecoTree->Branch("pointVtxDeltaLong",&fPointVtxDeltaLong,"pointVtxDeltaLong/D");
00690     fRecoTree->Branch("pointVtxDeltaAngle",&fPointVtxDeltaAngle,"pointVtxDeltaAngle/D");
00691 
00692     fRecoTree->Branch("extendedVtxX",&fExtendedVtxX,"extendedVtxX/D");
00693     fRecoTree->Branch("extendedVtxY",&fExtendedVtxY,"extendedVtxY/D");
00694     fRecoTree->Branch("extendedVtxZ",&fExtendedVtxZ,"extendedVtxZ/D");
00695     fRecoTree->Branch("extendedDirX",&fExtendedDirX,"extendedDirX/D");
00696     fRecoTree->Branch("extendedDirY",&fExtendedDirY,"extendedDirY/D");
00697     fRecoTree->Branch("extendedDirZ",&fExtendedDirZ,"extendedDirZ/D"); 
00698     fRecoTree->Branch("extendedConeAngle",&fExtendedConeAngle,"extendedConeAngle/D");
00699     fRecoTree->Branch("extendedTrackLength",&fExtendedTrackLength,"extendedTrackLength/D");
00700     fRecoTree->Branch("extendedVtxFOM",&fExtendedVtxFOM,"extendedVtxFOM/D");
00701     fRecoTree->Branch("extendedVtxSteps",&fExtendedVtxSteps,"extendedVtxSteps/I");
00702     fRecoTree->Branch("extendedVtxPass",&fExtendedVtxPass,"extendedVtxPass/I");
00703     fRecoTree->Branch("extendedVtxDelta",&fExtendedVtxDelta,"extendedVtxDelta/D");
00704     fRecoTree->Branch("extendedVtxDeltaTrans",&fExtendedVtxDeltaTrans,"extendedVtxDeltaTrans/D");
00705     fRecoTree->Branch("extendedVtxDeltaLong",&fExtendedVtxDeltaLong,"extendedVtxDeltaLong/D");
00706     fRecoTree->Branch("extendedVtxDeltaAngle",&fExtendedVtxDeltaAngle,"extendedVtxDeltaAngle/D");
00707 
00708     // --- debug ---
00709     fRecoTree->Branch("extendedVtxTimeParam0",&fExtendedVtxTimeParam0,"extendedVtxTimeParam0/D");
00710     fRecoTree->Branch("extendedVtxConeParam0",&fExtendedVtxConeParam0,"extendedVtxConeParam0/D");
00711     fRecoTree->Branch("extendedVtxConeParam1",&fExtendedVtxConeParam1,"extendedVtxConeParam1/D");
00712     fRecoTree->Branch("extendedVtxConeParam2",&fExtendedVtxConeParam2,"extendedVtxConeParam2/D");
00713 
00714     gDirectory = tmpd;
00715   }
00716 
00717   if( fRecoFile ){
00718     tmpd = gDirectory;
00719     fRecoFile->cd();
00720     fRecoTree->Fill();
00721     gDirectory = tmpd;
00722   }
00723 
00724   return;
00725 }
00726 
00727 void WCSimRecoNtuple::CloseFile()
00728 {
00729   TDirectory* tmpd = 0;
00730 
00731   if( fRecoFile ){
00732     tmpd = gDirectory; 
00733     std::cout << " *** WCSimRecoNtuple::CloseFile() *** " << std::endl;
00734     std::cout << "  closing file: " << fNtpFileName.Data() << std::endl;
00735     fRecoFile->cd();
00736     fRecoTree->Write();
00737     fRecoFile->Close();
00738     gDirectory = tmpd;
00739   }
00740   
00741   return;
00742 }