StRoot: macros/analysis/find_vertex.C Source File

00001 // 
00002 //   macro: find_vertex.C
00003 // 
00004 //   author: Marco van Leeuwen
00005 // 
00006 //   notes: placed in CVS by Gene Van Buren
00007 //          Need to load libs before running this macro,
00008 //          because of function with StEvent, StPrimaryVertex
00009 //          When you run this, it will actually add 
00010 //          new vertices to the event, so you have to
00011 //          keep track of how many were there when
00012 //          you read the file at first; I think the
00013 //          new vertices are added after the old one.
00014 // 
00015 
00016 //class StEmcDetector;
00017 class StChain;
00018 //class StEmcRawData;
00019 //class EEfeeDataBlock;
00020 //class Collection;
00021 //class StSPtrVecTrackNodeIterator ;
00022 //class StTriggerIdCollection;
00023 class StMinuitVertexFinder;
00024 class StGenericVertexFinder;
00025 class StGenericVertexMaker;
00026 
00027 StChain *chain=0;
00028 TH1F *dca_z_h=0;
00029 TNtuple *vtx_tuple=0;
00030 
00031 float eval_print_vertex(StEvent *event, StPrimaryVertex *vtx) {
00032   cout << "Vertex at " << vtx->position() << " chisq " << vtx->chiSquared() << endl;
00033   cout << "CTB matches: " << vtx->numMatchesWithCTB() << " BEMC: " 
00034        << vtx->numMatchesWithBEMC() << " cross central membrane " 
00035        << vtx->numTracksCrossingCentralMembrane() 
00036        << " rank " << vtx->ranking() << endl;
00037   int cnt_trk=0;
00038   int cnt_dcaz=0;
00039   int cnt_dca=0;
00040   Int_t n_node = event->trackNodes().size();
00041   Double_t avg_dip = 0;
00042   Double_t rms_dip = 0;
00043   for (int i_node = 0; i_node < n_node; i_node++) {
00044     StTrackNode *node = event->trackNodes()[i_node];
00045     StTrack *track = node->track(global);
00046     if (track &&
00047         track->flag() >= 0 &&
00048         track->fitTraits().numberOfFitPoints() >= 15 &&
00049         !track->topologyMap().trackFtpc() &&
00050         TMath::Finite(track->length()) ) {
00051       //cout << "track " << track << endl;
00052       cnt_trk++;
00053       double pathlength = track->geometry()->helix().pathLength( vtx->position(), false ); // do not scan periods
00054       StThreeVectorF dca = track->geometry()->helix().at(pathlength)-vtx->position(); 
00055       dca_z_h->Fill(dca.z());
00056       if (fabs(dca.z()) < 3) {
00057         cnt_dcaz++;
00058         float dip = track->geometry()->helix().dipAngle();
00059         avg_dip += dip;
00060         rms_dip += dip * dip;
00061       } 
00062       if (dca.mag() < 3) {
00063         cnt_dca++;
00064       } 
00065     }
00066   }
00067   if (cnt_dcaz) {
00068     avg_dip /= cnt_dcaz;
00069     rms_dip = sqrt(rms_dip - cnt_dcaz*avg_dip*avg_dip) / cnt_dcaz;
00070   }
00071   cout << "Tracks used in vtx " << vtx->numTracksUsedInFinder() << ", num with dcaz cut " << cnt_dcaz << ", dca < 3 " << cnt_dca << endl;
00072   cout << "avg dip " << avg_dip << ", rms dip " << rms_dip << endl;
00073   vtx_tuple->Fill(vtx->numTracksUsedInFinder(),vtx->position().z(),vtx->chiSquared(),cnt_dcaz,cnt_dca,vtx->numMatchesWithBEMC(),vtx->numTracksCrossingCentralMembrane(),avg_dip,rms_dip);
00074 
00075   return vtx->chiSquared();
00076 }
00077 
00078 void find_vertex(char * fname="high_053/st_physics_6053108_raw_2020002.event.root", Int_t nevents=1000){
00079 
00080   //fname="high_053/st_physics_6053108_raw_2020002.event.root";// daq1
00081   //fname="low_mb_49/st_physics_adc_6049129_raw_3080002.event.root";
00082 
00083   //      ROOT libs
00084     gSystem->Load("libPhysics");
00085     gSystem->Load("libTable");
00086     gSystem->Load("libGeom");
00087 
00088 //      STAR libs  
00089     gSystem->Load("StarMagField");
00090     gSystem->Load("St_base");
00091     gSystem->Load("StChain");
00092     gSystem->Load("St_Tables");
00093     gSystem->Load("StUtilities");        // new addition 22jul99
00094     gSystem->Load("StTreeMaker");
00095     gSystem->Load("StIOMaker");
00096     gSystem->Load("StarClassLibrary");
00097     gSystem->Load("StTriggerDataMaker"); // new starting from April 2003
00098     gSystem->Load("StBichsel");
00099     gSystem->Load("StEvent");
00100     gSystem->Load("StEventUtilities");
00101     gSystem->Load("StEmcUtil"); 
00102     gSystem->Load("StTofUtil");
00103     gSystem->Load("StPmdUtil");
00104     gSystem->Load("StPreEclMaker");
00105     gSystem->Load("StStrangeMuDstMaker");
00106     gSystem->Load("StMuDSTMaker");
00107     gSystem->Load("StMagF");
00108 
00109     gSystem->Load("StDbLib.so");
00110     gSystem->Load("StDbBroker.so");
00111     gSystem->Load("libStDb_Tables.so");
00112     gSystem->Load("St_db_Maker.so");
00113 
00114 
00115     cout << " loading of shared libraries done" << endl;
00116 
00117 
00118   //gROOT->LoadMacro("$STAR/StRoot/StMuDSTMaker/COMMON/macros/loadSharedLibraries.C");
00119   //loadSharedLibraries();
00120  
00121   // Load my makers
00122  
00123   gSystem->Load("Sti");
00124   gSystem->Load("libStEEmcUtil");
00125   gSystem->Load("StGenericVertexMaker");
00126   
00127   TFile *fout = new TFile("vtx_tree.root","RECREATE");
00128   dca_z_h = new TH1F("dca_z_h","dz_z_h",100,-50,50);
00129   vtx_tuple = new TNtuple("vtx_tuple","vertex info","num_trk:vtx_z:chisq:n_dcaz:n_dca:n_bemc:n_cross:avg_dip:rms_dip");
00130 
00131   // book some histos
00132   //gSystem->Load("libStGenericVertexMaker");
00133   // create chain    
00134   chain = new StChain("bfc"); 
00135   //chain->SetDebug();
00136   
00137   // Now we add Makers to the chain...
00138 
00139   // StIOMaker - to read files ...
00140   StIOMaker* ioMaker = new StIOMaker();  
00141   //ioMaker->SetFile("photon_bemc.event.root"); 
00142   ioMaker->SetFile(fname); 
00143   //ioMaker->SetDebug();
00144   ioMaker->SetIOMode("r");
00145   ioMaker->SetBranch("*",0,"0");             //deactivate all branches
00146   ioMaker->SetBranch("eventBranch",0,"r");   //activate Event Branch
00147   ioMaker->SetIOMode("r"); 
00148  
00149   St_db_Maker* dbMk = new St_db_Maker("db","MySQL:StarDb","$STAR/StarDb","StarDb");
00150 
00151   //StGenericVertexMaker *myfinder=new StGenericVertexMaker("myvertexfinder");
00152   
00153   StGenericVertexMaker *myfinder=new StGenericVertexMaker("myvertexfinder");
00154   //StGenericVertexMaker *myfinder=new StMinuitVertexFinder("myvertexfinder");
00155   //StMaker *myfinder=new StMinuitVertexFinder("myvertexfinder");
00156   myfinder->SetDebug(1);
00157   myfinder->SetMode(1);  // Select Minuit finder
00158   myfinder->SetInternalFind();  // Activate finding
00159   myfinder->Init();
00160   
00161   chain->PrintInfo();
00162   //chain->ls(3);
00163   Int_t initStat = chain->Init(); // This should call the Init() method in ALL makers
00164   if (initStat) chain->Fatal(initStat, "during Init()");
00165     
00166   int istat=0,iev=0;
00167 
00168   // Do the event loop    
00169   while(1) {
00170     if (iev>=nevents) break;
00171     chain->Clear();
00172     cout << "---------------------- Processing Event : " << iev << " ---------------------- " << endl;
00173    istat = chain->Make();
00174     iev++;
00175     if(istat) break;
00176     cout << "istat " << istat<<endl;
00177 
00178     // if(iev<17) continue;
00179     if (istat  == kStEOF || istat == kStFatal) break;
00180 
00181     StEvent* mEvent = (StEvent*)chain->GetInputDS("StEvent");
00182     assert(mEvent);// fix your chain or open the right event file
00183 
00184         //StMinuitVertexFinder* test_finder=(StMinuitVertexFinder*) myfinder->GetGenericFinder();
00185         //myfinder->GetGenericFinder().printInfo();
00186         //StMinuitVertexFinder test_finder;//=new StMinuitVertexFinder();
00187         //test_finder->setPrintLevel(3);
00188         //StThreeVectorD myvertex;
00189         //if (myfinder->fit(mEvent)) {
00190         //myvertex = myfinder->result();
00192           //}
00193           //else
00194           //      cout << "Error: vertex fit failed, no vertex." << endl;
00195           
00196           //delete myfinder;
00197           
00198           //myfinder->Finish();
00199           Int_t nV=mEvent->numberOfPrimaryVertices();
00200           if (nV == 0) continue;
00201           int iv;
00202           Float_t best_rank=1e9;
00203           StPrimaryVertex *best_vtx=0;
00204           for(iv=0;iv<nV;iv++) {
00205             StPrimaryVertex *V=mEvent->primaryVertex(iv);
00206             assert(V);
00207             StThreeVectorF &r=V->position();
00208             StThreeVectorF &er=V->positionError();
00209             Float_t rank=eval_print_vertex(mEvent, V);
00210             if (best_vtx==0 || rank < best_rank) {
00211               best_vtx = V;
00212               best_rank = rank;
00213             }
00214           }
00215           if (best_vtx)
00216             cout << "Best vertex: " << *best_vtx << endl; 
00217 //      printf("  nPrimTr=%d , VFid=%d:: ntrVF=%d nCtb=%d nBemc=%d nEEmc=%d nTpc=%d sumPt=%.1f rank=%g\n"
00218 //             ,V->numberOfDaughters(), V->vertexFinderId() ,V->numTracksUsedInFinder()  ,
00219 //             V->numMatchesWithCTB()  ,V-> numMatchesWithBEMC() ,V->numMatchesWithEEMC()  ,
00220 //             V->numTracksCrossingCentralMembrane()  ,V->sumOfTrackPt()  ,V->ranking());
00221         
00222 //      continue;
00223 //      int nPrTr=0;
00224 //      //.... access prim tracks for given vertex
00225 //      int itr;
00226 //      for(itr=0; itr<V->numberOfDaughters(); itr++) {
00227 //        StTrack *track=V-> daughter(itr);
00228 //        if(track==0)  continue;
00229 //        if (track->flag() <0 ) continue;
00230 //        printf("itr=%d pT=%.1f eta=%.2f nFitP=%d DCA=%.1f\n",itr,
00231 //               track->geometry()->momentum().mag(),
00232 //               track->geometry()->momentum().pseudoRapidity(),
00233 //               track->fitTraits().numberOfFitPoints(),
00234 //               track->geometry()->helix().distance(V->position()));
00235 //        nPrTr++;
00236 //      }
00237         
00238 //      printf("  counted nPrimTr=%d \n",nPrTr);
00239 //       } // end of loop over vertices
00240  
00241   
00242   } // Event Loop
00243   //chain->Finish();
00244   fout->Write();
00245 }
00246 
00247 
00248 //
00249 // $Id: find_vertex.C,v 3.1 2008/03/19 19:39:08 genevb Exp $
00250 // $Log: find_vertex.C,v $
00251 // Revision 3.1  2008/03/19 19:39:08  genevb
00252 // Introduction of macro
00253 //
00254 //