StRoot: StEEmcPool/LCP/StMuLcp2TreeMaker.cxx Source File

00001 
00002 // *-- Author : Jan Balewski
00003 // 
00004 // $Id: StMuLcp2TreeMaker.cxx,v 1.7 2009/12/02 16:35:58 fine Exp $
00005 
00006 #include <TFile.h>
00007 #include <TH2.h>
00008 
00009 #include "StMuLcp2TreeMaker.h"
00010 
00011 #include "TChain.h"
00012 #include "TClonesArray.h"
00013 #include "StL0Trigger.h"
00014 #include "StEventInfo.h"
00015 #include "StEventSummary.h"
00016 #include "StarClassLibrary/StThreeVectorF.hh"
00017 #include "StMuDSTMaker/COMMON/StMuEvent.h"
00018 #include "StMuDSTMaker/COMMON/StMuTrack.h"
00019 #include "StMuDSTMaker/COMMON/StMuDstMaker.h"
00020 //#include "MikesRejector/ExampleUsage.h"
00021 
00022 #include "CtbMatching.h"
00023 
00024 #include "StThreeVector.hh"
00025 
00026 ClassImp(StMuLcp2TreeMaker)
00027 
00028 
00029 //________________________________________________
00030 //________________________________________________
00031 StMuLcp2TreeMaker::StMuLcp2TreeMaker(const char* self ,const char* muDstMakerName) : StMaker(self){
00032   mMuDstMaker = (StMuDstMaker*)GetMaker(muDstMakerName);
00033   assert(mMuDstMaker);
00034 
00035   // rejector=new  ExampleUsage; // Mike's cosmic/laser rejector
00036   primTrA=0;
00037 
00038   runID=101;// default not initialized value
00039   treeName="./";
00040   eve_cosm=-8;
00041   SetMinNFitPoint(15); 
00042   SetMaxDCAxy(3.);// cm
00043   SetMaxEta(1.4);
00044   SetMinPt(0.4); // GeV/c
00045   SetMinFitPfrac(0.55);
00046   SetMaxZvert(100.); // (cm)
00047   ctb=new CtbMatching();
00048   eve_nGlob=-1;
00049 }
00050 
00051 
00052 //___________________ _____________________________
00053 //________________________________________________
00054 StMuLcp2TreeMaker::~StMuLcp2TreeMaker(){
00055   if(tree) tree->Print();
00056  
00057  // Save all objects in this file
00058  hfile->Write();
00059  
00060 }
00061 
00062 
00063 //________________________________________________
00064 //________________________________________________
00065 Int_t StMuLcp2TreeMaker::InitRunFromMake  (int runNumber){
00066   float Pi=3.14159;
00067   if(runID==runNumber)  return kStOK;
00068   if(runID==888999)  return kStOK;// special run number
00069   
00070   if(runID !=101) {
00071     printf("\n\n FATAL %s::InitRunFromMake() runID=%d changed to %d, (not implemented),STOP \n",GetName(),runID,runNumber);
00072     assert(runID==runNumber);
00073   }
00074   
00075   runID=runNumber;
00076   
00077   treeName+="/R";
00078   treeName+=runID;
00079   treeName+=".tree.root";
00080   printf("%s::InitRunFromMake() \n   runID=%d Setup output tree & histo to '%s' , use BXoff48=%d , CUTS: minNFitPoint=%d maxDCAxy/cm=%.2f maxEta=%.2f minPt/GeV/c=%.2f, minFitPfrac=%.2f maxZvert=%f\n",GetName(),runID,treeName.Data(),off48,C_minNFitPoint,C_maxDCAxy,C_maxEta,C_minPt,C_minFitPfrac,C_maxZvertex);
00081   
00082   hfile=new TFile(treeName,"RECREATE"," histograms & trees with LCP ");
00083   assert(hfile->IsOpen());
00084   
00085   // Create some histograms and a profile histogram
00086   h[0]=new TH1F("bX7","rates vs. 7bit bXing",129,-0.5,128.5);
00087   h[1]=new TH1F("bX48","rates vs. 48bit bXing +offset%120 ",129,-0.5,128.5);
00088   h[2]=new TH1F("bX120","rates vs. true bXing [0-119]",129,-0.5,128.5);
00089   h[3]=(TH1F *)new TH2F("sBit120","spinBits vs.  bXing120",129,-0.5,128.5,17,-0.5,16.5);
00090   h[4]   = new TH1F("Vz","Vertex Z/cm",100,-250,250.);
00091   h[5]   = new TH1F("nPrim","No. of prim TPC tracks",50,-0.5,49.5);
00092   h[6]   =new TH1F("eta","eta of LCP",100,-2.,2.);
00093 
00094   h[7]   = new TH1F("phi"," phi/rad  of LCP", 240,-Pi,Pi);
00095   h[8]   = new TH1F("pT","pT/GeV of LCP",100,0,10.);
00096   h[9]   = new TH1F("pTm1","pT/GeV of LCP, CTB match any",100,0,10.);
00097   h[10]   = new TH1F("pTm2","pT/GeV of LCP, CTB match LCP",100,0,10.);
00098   
00099 
00100   if (runNumber>100) {   
00101     // Create tree
00102     tree = new TTree("T-LCP","LCP selected from MinB trigger");
00103     tree->Branch("bx48", &eve_bx48,"bx48/I");
00104     tree->Branch("bx120",&eve_bx120,"bx120/I");
00105     tree->Branch("id",   &eve_id,"id/I");
00106     tree->Branch("sb",   &eve_sb,"sb/I");
00107     tree->Branch("nPrim",&eve_nPrim,"nPrim/I");
00108     tree->Branch("nGlob",&eve_nGlob,"nGlob/I");
00109     tree->Branch("vz",   &eve_vz ,"vz/F");
00110     tree->Branch("cosm", &eve_cosm ,"cosm/F");
00111     
00112     tree->Branch("pt",&lcp_pt ,"pt/F");
00113     tree->Branch("phi",&lcp_phi ,"phi/F");
00114     tree->Branch("eta",&lcp_eta ,"eta/F");
00115     tree->Branch("q",&lcp_q ,"q/I");
00116     tree->Branch("nFit",&lcp_nFit ,"nFit/I");
00117   } else {
00118     printf(" %s::InitRunFromMake(RunNumber=%d) WARN : TTRee NOT created\n",GetName(),runNumber);
00119     tree=0;
00120   }
00121 
00122 #if 0
00123   // cuts study
00124   int nd=0;  float x1=1,x2=2; TString titCore="fixMe2";
00125   
00126   //  titCore="nFitPoint"; nd=15; x1=9.5; x2=39.5;
00127   //  titCore="maxDCAxy/cm"; nd=15; x1=0.; x2=3.0;
00128   // titCore="maxEta"; nd=15; x1=0.35; x2=1.85;
00129   titCore="minFitPfrac"; nd=10; x1=0.325; x2=0.825;
00130   
00131   hc[0]=new TH1F("CnFP-A","Yield vs."+titCore+", for LCP #Delta#phi <1/8 #pi ",nd,x1,x2);
00132   hc[1]=new TH1F("CnFP-B","Yield vs."+titCore+", for  LCP #Delta#phi #in [ 1/8 #pi, 3/8 #pi ] ",nd,x1,x2);
00133   hc[2]=new TH1F("CnFP-C","Yield vs."+titCore+", for LCP #Delta#phi #in [ 3/8 #pi, 5/8 #pi ] ",nd,x1,x2);
00134   hc[3]=new TH1F("CnFP-D","Yield vs. "+titCore+", for LCP #Delta#phi #in [ 5/8 #pi, 7/8 #pi ] ",nd,x1,x2);
00135   hc[4]=new TH1F("CnFP-E","Yield vs."+titCore+", for  LCP #Delta#phi > 7/8 #pi ",nd,x1,x2);
00136   hc[5]=new TH1F("CnFP-L","Yield vs."+titCore+", lost LCP",nd,x1,x2);
00137   hc[6]=new TH1F("CnFP-W","Yield vs."+titCore+", won LCP",nd,x1,x2);
00138 #endif
00139   
00140 
00141   return kStOK;
00142 }
00143  
00144 //________________________________________________
00145 //________________________________________________
00146 Int_t StMuLcp2TreeMaker::Init(){
00147   return StMaker::Init();
00148 }
00149  
00150 
00151 //________________________________________________
00152 //________________________________________________
00153 void StMuLcp2TreeMaker::clearLCP(){
00154   lcp_eta=lcp_phi=lcp_pt=0;
00155   lcp_q=0;
00156   lcp_nFit=0;
00157 }
00158  
00159 
00160 //________________________________________________
00161 //________________________________________________
00162 Int_t StMuLcp2TreeMaker::Make(){
00163   
00164   static int kEve=0;
00165   kEve++;
00166   
00167     printf("%s::Make() is called ..........\n",GetName());
00168   
00169   StMuDst* dst = mMuDstMaker->muDst();
00170   primTrA=dst->primaryTracks();
00171     
00172   StMuEvent* muEve = dst->event();
00173   StL0Trigger &trig=muEve->l0Trigger();
00174   StEventInfo &info=muEve->eventInfo();
00175   StEventSummary &smry=muEve->eventSummary();
00176   ctb->loadHits(muEve);
00177   StThreeVectorF ver=smry.primaryVertexPosition();
00178   if(runID>100) InitRunFromMake(info.runId());
00179 
00180   //  printf("trgWrd=%d nPrim=%d, Vz=%f \n", trig.triggerWord(),primTrA->GetEntries(),ver.z());
00181   // use only minB trigger events 
00182    if(trig.triggerWord()!=0x2000)  return kStOK;
00183   
00184   // reject events without primtracks
00185      if(primTrA->GetEntries()<=0)  return kStOK;
00186  
00187   // fill eve-related tree values
00188   eve_vz=ver.z(); 
00189      if(fabs(eve_vz)>C_maxZvertex) return kStOK;
00190 
00191   int bx7=trig.bunchCrossingId7bit(info.runId());
00192   eve_id=info.id();
00193   eve_bx48 = trig.bunchCrossingId();
00194   eve_bx120 =(trig.bunchCrossingId()+off48 )%120;
00195   eve_sb=trig.spinBits(info.runId());
00196  
00197   // test for false vertex
00198   //eve_cosm=rejector->acceptEvent(dst);
00199   eve_cosm=0;
00200 
00201   // search for LCP  
00202   clearLCP();
00203   const StMuTrack* lcp= findLCP(C_minPt,C_minNFitPoint, C_maxDCAxy, C_maxEta,C_minFitPfrac);
00204 
00205   // fill histo 1:1
00206   h[0]->Fill(bx7);
00207   h[1]->Fill(eve_bx48);
00208   h[2]->Fill(eve_bx120);   
00209   ((TH2F*)h[3])->Fill(eve_bx120,eve_sb);
00210   h[4]->Fill(eve_vz);
00211   h[5]->Fill(eve_nPrim); // <<-- assigned by findLCP()
00212   
00213   if(lcp) { // LCP was found
00214     lcp_eta=lcp->eta();
00215     lcp_phi=lcp->phi();
00216     lcp_pt= lcp->pt();
00217     lcp_q = lcp->charge();
00218     lcp_nFit=lcp->nHitsFit();     
00219 
00220     // fill those histo only for LCP
00221     h[6]->Fill(lcp_eta);
00222     h[7]->Fill(lcp_phi);
00223     h[8]->Fill(lcp_pt);
00224     if(eve_cosm) h[9]->Fill(lcp_pt);
00225     if(eve_cosm>1000) h[10]->Fill(lcp_pt);
00226   }
00227   
00228   // scan for good globals, to match Olg's analysis
00229   {
00230     //    TClonesArray*  globTrA=dst->globalTracks();
00231     TObjArray*  globTrA=dst->primaryTracks();
00232     //printf("tot Glob=%d\n", globTrA->GetEntries());
00233     int nTr=0;
00234     for (int i=0; i<globTrA->GetEntries();i++) {
00235       TObject &oo=globTrA[i];
00236       const StMuTrack* gTr = (StMuTrack*)(&oo);// perhaps it works? JB
00237       if(gTr->flag() <=0) continue;
00238       StThreeVectorF dca=gTr->dcaGlobal();
00239       if(!gTr->topologyMap().trackTpcOnly()) continue;
00240       if(gTr->nHitsFit()<15 )continue;
00241       if(dca.perp()>3.) continue;  
00242       if(fabs(dca.z())>3.) continue;
00243       //printf("nFit=%d,DCA  x=%f y=%f z=%f zV=%f\n",gTr->nHitsFit(),dca.x(),dca.y(),dca.z(),ver.z());
00244       nTr++;
00245     }
00246     eve_nGlob=nTr;
00247   }
00248   
00249   
00250   if(tree)tree->Fill();
00251   //if(kEve%500==0)
00252     printf("#%s::Make(%d) nPrim=%d nGlob=%d lcp: pt=%f nFit=%d eta=%f\n",GetName(),kEve,eve_nPrim, eve_nGlob, lcp_pt,lcp_nFit,lcp_eta);
00253   
00254   // extra tests of cuts
00255   //   examinCut(lcp);
00256   
00257   return kStOK;
00258 }
00259 
00260 //________________________________________________
00261 //________________________________________________
00262 const StMuTrack* StMuLcp2TreeMaker::findLCP( float XminPt,int XminNFitP , float XmaxDCAxy, float XmaxEta, float XminFitPfrac) {
00263   
00264   eve_cosm=0;
00265   const StMuTrack* lp=0;  
00266   int nTr=primTrA->GetEntries();
00267   int nTr0=0,nTr1=0; // counts valid prim TPC tracks
00268   int lpMatch=0;
00269   float maxPt=XminPt;// init pT  
00270   for (int i=0; i<nTr; i++) {
00271     TObject &oo=primTrA[i];
00272     const StMuTrack* pTr = (StMuTrack*)(&oo);// perhaps it works? JB
00273     if(pTr->flag() <=0) continue;
00274     if(!pTr->topologyMap().trackTpcOnly()) continue;
00275     if(pTr->globalTrack()->nHitsFit()>15 ) nTr0++; // was just: nTr0++;
00276     // printf("itr=%d pT=%f nFit=%d\n",i, pTr->pt(),pTr->nHitsFit());
00277     
00278     StThreeVectorF dca=pTr->dcaGlobal();
00279     //printf("xyz=%f %f %f, pt=%f\n",dca.x(),dca.y(),dca.z(),dca.perp());
00280     if(dca.perp()>XmaxDCAxy) continue;
00281 
00282     // switch to oryginal glob-track, since vertex constrain biases pt 
00283     const StMuTrack* gTr=pTr->globalTrack();
00284     assert(gTr);// every prim must have its global
00285         
00286     if ( gTr->nHitsFit() < XminNFitP) continue;
00287     if ( fabs(gTr->eta()) > XmaxEta ) continue;
00288     float frac=1.*gTr->nHitsFit()/gTr->nHitsPoss();
00289     if ( frac < XminFitPfrac) continue;
00290     nTr1++;
00291     // printf("itr=%d pT=%f phi/deg=%f eta=%f \n",i, pTr->pt(),pTr->phi()/3.1416*180,pTr->eta());
00292 
00293     int match=ctb->match(gTr);
00294     if(match>0) eve_cosm++;
00295     // printf("itr=%d pT=%f nFit=%d match=%d\n",i, pTr->pt(),pTr->nHitsFit(),match);
00296     // now search for high pT
00297     if ( gTr->pt() < maxPt ) continue;
00298     maxPt=gTr->pt();
00299     lp=gTr;
00300     if(match>0)lpMatch++;
00301   }
00302   // tmp - apply pT cut on LCP
00303   //if(maxPt<1.5 || maxPt>3) lp=0;
00304 
00305   eve_nPrim=nTr0;
00306   if(lpMatch) eve_cosm+=1000;
00307   //  printf("nTr0=%d nTr1=%d lp=%p match=%d\n",nTr0,nTr1,lp,eve_cosm);
00308   return lp;
00309 }
00310 
00311 
00312 //________________________________________________
00313 //________________________________________________
00314 void StMuLcp2TreeMaker::examinCut(const StMuTrack*lcp0){
00315 
00316   const float Pi=3.1416;
00317   
00318   StThreeVectorF p0;
00319   float phi0=0;
00320 
00321   if(lcp0) {
00322     p0=lcp0->p();
00323     phi0=Pi/2+ atan2(p0.y(),p0.x());
00324   }
00325  
00326   int nb=hc[0]->GetNbinsX();
00327  
00328   int i;
00329   for(i=1;i<=nb;i++) {// vary nFitP ......................
00330     float x=hc[0]->GetBinCenter(i);
00331     const StMuTrack* lcp1=0;
00332 
00333     // lcp1=findLCP(C_minPt,(int)x,C_maxDCAxy,C_maxEta,C_minFitPfrac);   // test nFitPoint
00334     //lcp1=findLCP(C_minPt,C_minNFitPoint,x,C_maxEta,C_minFitPfrac);       // test maxDCAxy
00335     // lcp1=findLCP(C_minPt,C_minNFitPoint,C_maxDCAxy,x,C_minFitPfrac);   // test maxEta
00336     lcp1=findLCP(C_minPt,C_minNFitPoint,C_maxDCAxy,C_maxEta,x);     // test fitPfrac
00337     //printf("i=%d x=%f %p %p\n",i,x,lcp0,lcp1);
00338  
00339     float delPhi=Pi*1.1;
00340     if(lcp0 && lcp1==0) {    // lost LCP
00341       hc[5]->Fill(x);
00342     } else  if(lcp1 && lcp0==0) { // won LCP
00343       hc[6]->Fill(x);
00344     } else if(lcp1 && lcp0) {  // phi could change
00345       StThreeVectorF p1=lcp1->p();
00346       float phi1=Pi/2+ atan2(p1.y(),p1.x());
00347       delPhi=phi1-phi0;
00348       if(delPhi <0) delPhi*=-1;
00349       if(delPhi >Pi) delPhi=2*Pi-delPhi;
00350       // printf("(deg) phi0=%.1f, phi1=%.1f, del=%.1f\n",phi0/Pi*180, phi1/Pi*180, delPhi/Pi*180);
00351       if( delPhi <1./8*Pi) {
00352         hc[0]->Fill(x);
00353       } else if( delPhi <3./8*Pi) {
00354         hc[1]->Fill(x);
00355       } else if( delPhi <5./8*Pi) {
00356         hc[2]->Fill(x);
00357       } else if( delPhi <7./8*Pi) {
00358         hc[3]->Fill(x);
00359       } else {
00360         hc[4]->Fill(x);
00361       }
00362     }
00363   }
00364 }
00365 
00366 
00367 
00368 
00369 // $Log: StMuLcp2TreeMaker.cxx,v $
00370 // Revision 1.7  2009/12/02 16:35:58  fine
00371 // Fix StMuTrack interface
00372 //
00373 // Revision 1.6  2005/10/11 17:28:30  balewski
00374 // fix related L0trig code change
00375 //
00376 // Revision 1.5  2005/08/23 21:09:23  balewski
00377 // fix to follow muDst evolution
00378 //
00379 // Revision 1.4  2004/01/06 17:25:26  balewski
00380 // get LCP from Geant info
00381 //
00382 // Revision 1.3  2003/11/12 18:43:41  balewski
00383 // final for LCP paper
00384 //
00385 // Revision 1.2  2003/10/20 17:04:39  balewski
00386 // LCP analysis code
00387 //
00388 // Revision 1.1  2003/09/16 19:18:35  balewski
00389 // extraction of LCP from muDst
00390 //