Kinfit: h1-k0-2.C Source File

00001 
00002 #define H1Tracks_cxx
00003 
00004 #include <iostream>              // - std::cout
00005 #include <cmath>                  // - math
00006 #include <cstring>                // - add strings
00007 #include "H1Tracks.h"                // - the ntuple
00008 #include <TApplication.h>                // - the ntuple
00009 
00010 #include "jbltools/kinfit/H1K0Event2.h"
00011 #include "jbltools/kinfit/WWFitter.h"
00012 #include "jbltools/kinfit/ftypes.h"
00013 #include "jbltools/kinfit/hbook.h"
00014 #include "jbltools/kinfit/TrackFitObject.h"
00015 #include "jbltools/kinfit/VertexFitObject.h"
00016 
00017 int main(int argc, char** argv) {
00018    
00019   // ntuple variables
00020   FInteger IERR;
00021   FReal FITPROB, CHI2;
00022   FInteger NIT;
00023   FReal V0TRUE[3], V1TRUE[3];
00024   FReal V0FIT[3], V1FIT[3];
00025   FReal P0FIT[4], P1FIT[4], P2FIT[4];
00026   FReal MTRUE, MFIT;
00027   FReal RFIT;
00028   
00029   
00030   FInteger LREC=1024, ISTAT=0;
00031   hropen (21, "K0FIT", "h1-k0-2.hbook", "N", LREC, ISTAT);
00032   hbnt (999,"k0fit","");
00033   hbname (999, "FIT", IERR, "IERR:I");
00034   hbname (999, "FIT", FITPROB, "FITPROB:R");
00035   hbname (999, "FIT", CHI2, "CHI2:R");
00036   hbname (999, "FIT", NIT, "NIT:I");
00037   hbname (999, "FIT", V0TRUE[0], "V0TRUE(3):R");
00038   hbname (999, "FIT", V0FIT[0], "V0FIT(3):R");
00039   hbname (999, "FIT", V1FIT[0], "V1FIT(3):R");
00040   hbname (999, "FIT", P0FIT[0], "P0FIT(4):R");
00041   hbname (999, "FIT", P1FIT[0], "P1FIT(4):R");
00042   hbname (999, "FIT", P2FIT[0], "P2FIT(4):R");
00043   hbname (999, "FIT", MFIT,  "MFIT:R");
00044   hbname (999, "FIT", RFIT,  "RFIT:R");
00045  
00046   // open ntuple (file name is specified in ntuple.h at the moment!)
00047   H1Tracks nt;
00048   if (nt.fChain == 0) {
00049     std::cout << "ERROR opening ntuple" << std::endl;
00050     return 1;  
00051   }
00052   
00053   // needed for graphic
00054   TApplication theApp("main", &argc, argv);
00055   
00056   
00057   int nentries = int(nt.fChain->GetEntriesFast());
00058   for (int ievent = 0; ievent < nentries; ++ievent) {
00059     Int_t ix = nt.LoadTree (ievent);
00060     if (ix < 0) break;
00061     nt.fChain->GetEntry(ievent);
00062     
00063     std::cout << "Run " << nt.RunNumber 
00064               << ", event " << nt.EventNumber << std::endl;
00065     std::cout << "Ntracks = " <<    nt.Ntracks << std::endl;         
00066     for (int t1 = 1; t1 < nt.Ntracks; t1++) {
00067       for (int t2 = t1+1; t2 < nt.Ntracks; t2++) {
00068         if (nt.Par[5*t1]*nt.Par[5*t2] < 0) {
00069           std::cout << "t1 = " << t1 << ", t2 = " << t2 << std::endl;
00070           WWFitter fitter;
00071     
00072           H1K0Event2 evt (nt, t1, t2);
00073           int ierr = evt.fitEvent (fitter);
00074           if (ierr) std::cout << "IERR: " << ierr << std::endl;
00075 
00076             IERR = ierr;
00077             FITPROB = fitter.getProbability();
00078             CHI2    = fitter.getChi2();
00079             NIT     = fitter.getIterations();
00080     
00081           // vertices
00082           for (int i = 0; i < 3; ++i) {
00083             V0TRUE[i] = nt.Vtx[i];
00084             V0FIT[i] = evt.rectrack[0]->getTrajectoryPoint(0).getComponent(i);
00085             V1FIT[i] = evt.k0decvertex->getVertex().getComponent(i);
00086           }
00087     
00088           // 4-momenta
00089           for (int i = 0; i < 4; ++i) {
00090             P0FIT[i] = evt.rectrack[0]->getMomentum(0).getComponent(i);
00091             P1FIT[i] = evt.rectrack[1]->getMomentum(0).getComponent(i);
00092             P2FIT[i] = evt.rectrack[2]->getMomentum(0).getComponent(i);
00093           }
00094           MFIT = (evt.rectrack[1]->getMomentum(0) + evt.rectrack[2]->getMomentum(0)).getM();
00095           ThreeVector primaryvertex (nt.Vtx[0], nt.Vtx[1], nt.Vtx[2]);
00096           RFIT = (evt.k0decvertex->getVertex()-primaryvertex).getR();
00097           hfnt(999);
00098         }
00099       }
00100     }
00101     if (ievent%10 == 9) std::cout << "processed event " << ievent+1 << std::endl;
00102     if (ievent >= 100) break;
00103   }
00104   
00105   hrout(999);
00106   hrend("K0FIT");          
00107 
00108   return 0;
00109    
00110 }