StRoot: StEEmcPool/LCP/StGeant2LcpTreeMaker.cxx Source File

00001 // *-- Author : Jan Balewski
00002 // 
00003 // $Id: StGeant2LcpTreeMaker.cxx,v 1.4 2007/07/12 19:24:55 fisyak Exp $
00004 
00005 #include <math.h>
00006 #include <TFile.h>
00007 #include <TH2.h>
00008 
00009 #include "StGeant2LcpTreeMaker.h"
00010 
00011 #include "TChain.h"
00012 #include "TClonesArray.h"
00013 
00014 #include "StEventInfo.h"
00015 #include "StMuDSTMaker/COMMON/StMuEvent.h"
00016 #include "StMuDSTMaker/COMMON/StMuDstMaker.h"
00017 
00018 #include "St_DataSet.h"
00019 #include "St_DataSetIter.h"
00020 #include "tables/St_g2t_event_Table.h"
00021 #include "tables/St_particle_Table.h"
00022 
00023 #include "StMuDSTMaker/COMMON/StMuTrack.h"
00024 #include "TMath.h"
00025 ClassImp(StGeant2LcpTreeMaker)
00026 
00027 
00028 //________________________________________________
00029 //________________________________________________
00030 StGeant2LcpTreeMaker::StGeant2LcpTreeMaker(const char* self ,const char* muDstMakerName) : StMaker(self){
00031   mMuDstMaker = (StMuDstMaker*)GetMaker(muDstMakerName);
00032   assert(mMuDstMaker);
00033 
00034   runID=101;// default not initialized value 
00035   treeName="./";
00036   SetMaxEta(1.0);
00037   SetMinPt(0.4); // GeV/c
00038 }
00039 
00040 
00041 //___________________ _____________________________
00042 //________________________________________________
00043 StGeant2LcpTreeMaker::~StGeant2LcpTreeMaker(){
00044   if(tree) tree->Print();
00045   
00046   // Save all objects in this file
00047   hfile->Write();
00048   
00049 }
00050 
00051 
00052 //________________________________________________
00053 //________________________________________________
00054 Int_t StGeant2LcpTreeMaker::Init(){
00055 
00056  return StMaker::Init();
00057 }
00058  
00059 
00060 //________________________________________________
00061 //________________________________________________
00062 void StGeant2LcpTreeMaker::clearLCP(){
00063   lcp_eta=lcp_phi=lcp_pt=lcp_e=0;
00064   lcp_idhep=0;
00065 }
00066  
00067 //________________________________________________
00068 //________________________________________________
00069 void StGeant2LcpTreeMaker::InitRunFromMake  (int runNumber){
00070   if(runID==runNumber)  return ;
00071   if(runID==888999)  return ; // special run number
00072   
00073 
00074   if(runID !=101) {
00075     printf("\n\n FATAL %s::InitRunFromMake() runID=%d changed to %d, (not implemented),STOP \n",GetName(),runID,runNumber);
00076     assert(runID==runNumber);
00077   }
00078   
00079   runID=runNumber;
00080   float Pi=3.14159;
00081   treeName+="/G";
00082   treeName+=runID;
00083   treeName+=".tree.root";
00084 
00085   printf("%s::Init() \n   runID=%d Setup output tree & histo to '%s' , CUTS:  maxEta=%.2f minPt/GeV/c=%.2f\n",GetName(),runID,treeName.Data(),C_maxEta,C_minPt);
00086 
00087   hfile=new TFile(treeName,"RECREATE"," histograms & trees with LCP ");
00088   assert(hfile->IsOpen());
00089   
00090   // Create some histograms and a profile histogram
00091   h[0]=0;
00092   h[1]=0;
00093   h[2]=0;
00094   h[3]=0;
00095   h[4]=0;
00096   h[5]   = new TH1F("nPKpi","No. of stable p,pbar,K+/-,pi+/- with |eta|<1 and pT>0.4",50,-0.5,49.5);
00097   h[6]   =new TH1F("eta","eta of LCP",100,-2.,2.);
00098 
00099   h[7]   = new TH1F("phi"," phi/rad  of LCP", 240,-Pi,Pi);
00100   h[8]   = new TH1F("pT","pT/GeV of LCP",100,0,10.);
00101   h[9]   = 0;
00102   h[10]   = 0;
00103   
00104 
00105   if (runID>100) {   
00106     // Create tree
00107     tree = new TTree("G-LCP","LCP selected from Geant data");
00108     tree->Branch("id",   &eve_id,"id/I");
00109     tree->Branch("sub",   &eve_sub,"sub/I");
00110     tree->Branch("nPKPi",&eve_nPKPi,"nPKPi/I");    
00111     tree->Branch("pt",&lcp_pt ,"pt/F");
00112     tree->Branch("phi",&lcp_phi ,"phi/F");
00113     tree->Branch("eta",&lcp_eta ,"eta/F");
00114     tree->Branch("e",   &lcp_e ,"e/F");
00115   } else {
00116     printf(" %s::Init(RunNumber=%d) WARN : TTRee NOT created\n",GetName(),runID);
00117     tree=0;
00118   }
00119 return;
00120 }  
00121 
00122 //________________________________________________
00123 //________________________________________________
00124 Int_t StGeant2LcpTreeMaker::Make(){
00125   
00126   static int kEve=0;
00127   kEve++;
00128   clearLCP();
00129   
00130   printf("%s::Make() is called ..........\n",GetName());
00131   
00132   StMuDst* dst = mMuDstMaker->muDst();
00133     
00134   StMuEvent* muEve = dst->event();
00135   StEventInfo &info=muEve->eventInfo();
00136 
00137   if(runID>100) InitRunFromMake(info.runId());
00138 
00139   printf(" eve ID=%d\n",info.id());
00140     
00141  
00142   // Access to geant-tables .......................
00143   
00144   St_DataSet* Event = GetDataSet("geant");
00145   //Event->ls(3);
00146   St_DataSetIter geantDstI(Event);
00147   
00148   // Event generator info ........................
00149   St_g2t_event *Pg2t_event=(St_g2t_event *) geantDstI("g2t_event");
00150   
00151   g2t_event_st *g2t_event1=Pg2t_event->GetTable();
00152   printf("nr=%d \n",(int)Pg2t_event->GetNRows());
00153   int k1=       g2t_event1->eg_label;
00154   int k2=       g2t_event1->n_event;
00155   int k3=       g2t_event1->subprocess_id;
00156   
00157   printf("eg_label=%d n_event=%d subprocess_id=%d\n", k1,k2,k3);
00158   
00159   assert(info.id()==g2t_event1->n_event);
00160   eve_id=info.id();
00161   eve_sub=g2t_event1->subprocess_id;
00162 
00163   // This is an example to access the particle table directly.
00164   St_particle    *particleTabPtr    =  (St_particle    *) geantDstI("particle");
00165   
00166   findGeantLcp( particleTabPtr);
00167   
00168   if(tree)tree->Fill();
00169   
00170   //if(kEve%500==0)
00171     printf("#%s::Make(%d) nPKPi=%d  lcp: pt=%f idHep=%d \n",GetName(),kEve,eve_nPKPi, lcp_pt,lcp_idhep);
00172   
00173   
00174   return kStOK;
00175 }
00176 
00177 
00178 //================================================
00179 //================================================
00180 void StGeant2LcpTreeMaker::printTrack( particle_st* part) {
00181   cout << "   ist = " << part->isthep;
00182   cout << "   id = " << part->idhep << endl;
00183   cout << "   px = " << part->phep[0];
00184   cout << "   py = " << part->phep[1];
00185   cout << "   pz = " << part->phep[2] << endl;
00186   cout << "   e  = " << part->phep[3];
00187   cout << "   m  = " << part->phep[4] << endl;
00188   //cout << "   moth1  = " << part->jmohep[0] << endl;
00189   //cout << "   moth2  = " << part->jmohep[1] << endl;
00190     
00191 }
00192 
00193 //================================================
00194 //================================================
00195 particle_st* StGeant2LcpTreeMaker::findGeantLcp( St_particle    *tab) {
00196   assert(tab);
00197   particle_st* part = tab->GetTable();
00198   printf("tab dim=%d\n", (int)tab->GetNRows());
00199   int i;
00200   float maxPt=0.4;
00201   float maxEta=1.;
00202   float etaI=99;
00203   int iMax=-1;
00204   int ns=0, np=0,n2=0,n3=0,npkpi=0;
00205   for(i=0;i<tab->GetNRows();i++) {
00206     if(part[i].isthep!=1) continue;
00207     ns++;
00208     int id=part[i].idhep;
00209     int idM=id>0 ? id : -id ;
00210     if(idM!=211 && idM!=321 && idM!=2212) continue;
00211     // accept only p,pbar, pi+/-, K+/-
00212     np++;
00213     
00214     float px=part[i].phep[0];
00215     float py=part[i].phep[1];
00216     float pz=part[i].phep[2];
00217 
00218     float pt=sqrt(px*px+py*py);
00219     if(pt<0.001) continue; // to allow for math
00220     float eta= TMath::ASinH(pz/pt);
00221 
00222 #if 0
00223     {
00224       float phi=atan2(part[i].phep[1],part[i].phep[0]);
00225       printf("i=%d pT=%f phi/deg=%f eta=%f pz=%f idHep=%d\n",i,pt,phi/3.1416*180,eta,pz,id);
00226     }
00227 #endif
00228     //printf("px=%f py=%f pt=%f\n",px,py,pt);
00229     if(pt>0.4 && fabs(eta)<1) npkpi++;
00230     if(pt<maxPt) continue;
00231     n2++;
00232 
00233     if(fabs(eta)>maxEta) continue;
00234     // track is in eta range
00235 
00236     n3++;
00237     maxPt=pt;
00238     etaI=eta;
00239     iMax=i;
00240   }
00241 
00242   printf("#total %3d stable particles, np=%3d n2=%3d n3=%3d maxPt=%f iMax=%d eta=%f npkpi=%d\n",ns,np,n2,n3,maxPt,iMax,etaI,npkpi);
00243 
00244   eve_nPKPi=npkpi;
00245   h[5]->Fill( eve_nPKPi);
00246 
00247   if(iMax>=0) {
00248     printTrack(part+iMax);
00249     lcp_eta=etaI;
00250     lcp_phi=atan2(part[iMax].phep[1],part[iMax].phep[0]);
00251     lcp_pt=maxPt;
00252     lcp_idhep=part[iMax].idhep;
00253     lcp_e=part[iMax].phep[3];
00254     h[6]->Fill(lcp_eta);
00255     h[7]->Fill(lcp_phi);
00256     h[8]->Fill(lcp_pt);
00257     
00258     return part+iMax;
00259   }
00260 
00261   return 0;
00262 
00263 }
00264 
00265 
00266 
00267 // $Log: StGeant2LcpTreeMaker.cxx,v $
00268 // Revision 1.4  2007/07/12 19:24:55  fisyak
00269 // Add includes for TMath for ROOT 5.16
00270 //
00271 // Revision 1.3  2004/02/24 21:33:23  balewski
00272 // *** empty log message ***
00273 //
00274 // Revision 1.2  2004/01/26 22:49:37  perev
00275 // WarnOff
00276 //
00277 // Revision 1.1  2004/01/06 17:25:26  balewski
00278 // get LCP from Geant info
00279 //