StRoot: StEEmcPool/StEEmcPi0Mixer/StEEmcPi0Maker.cxx Source File

00001 #include "StEEmcPi0Maker.h"
00002 #include "StMessMgr.h"
00003 
00004 #include "StEEmcPair.h"
00005 
00006 #include "StMuDSTMaker/COMMON/StMuDstMaker.h"
00007 #include "StMuDSTMaker/COMMON/StMuDst.h"
00008 #include "StMuDSTMaker/COMMON/StMuEvent.h"
00009 #include "StMuDSTMaker/COMMON/StMuPrimaryVertex.h"
00010 #include "StMuDSTMaker/COMMON/StMuTriggerIdCollection.h"
00011 
00012 #include "StEvent/StEvent.h"
00013 #include "StEvent/StTriggerIdCollection.h"
00014 #include "StEvent/StTriggerId.h"
00015 
00016 #include "StarClassLibrary/StThreeVectorF.hh"
00017 
00018 #include "TH1F.h"
00019 #include "TH2F.h"
00020 
00021 #include "TFile.h"
00022 #include "TTree.h"
00023 #include "StEEmcMixEvent.h"
00024 
00025 ClassImp(StEEmcPi0Maker);
00026 
00027 //34567890123456789012345678901234567890123456789012345678901234567890123456789012345678901234567890
00028 //       1         2         3         4         5         6         7       *
00029 
00030 // ---------------------------------------------------------------------------
00031 StEEmcPi0Maker::StEEmcPi0Maker( const Char_t *name,
00032                                 StEEmcA2EMaker *a2e,
00033                                 StEEmcGenericClusterMaker *cl,
00034                                 StEEmcGenericPointMaker   *pt ) : StMaker(name)
00035 {
00036   mEEanalysis=a2e; assert(a2e);
00037   mEEclusters=cl;  assert(cl);
00038   mEEpoints=pt;    assert(pt);
00039   mCheckTrigger=false;
00040   setFile(0);
00041   mTree=0;
00042 }
00043 
00044 void   StEEmcPi0Maker::setFile( TFile *f ){ mFile=f; }
00045 TTree *StEEmcPi0Maker::tree(){ return mTree; }
00046 
00047 // ---------------------------------------------------------------------------
00048 Int_t StEEmcPi0Maker::Init()
00049 {
00050   hMass   = new TH2F("hMass","Diphoton invariant mass; M [GeV]",120,0.,1.2,20,0.,20.);
00051   hMass_cut = new TH2F("hMass_cut","Diphoton invariant mass; M [GeV]",120,0.,1.2,20,0.,20.);
00052   hMass_split = new TH2F("hMass_split","Mass for events w/ a split smd cluster",120,0.,1.2,20,0.,20.);
00053 
00054   hPT     = new TH1F("hPT","Diphoton p_{T}; p_{T} [GeV]",20,0.,20.);
00055   hPT_cut = new TH1F("hPT_cut","Diphoton p_{T}; p_{T} [GeV]",20,0.,20.);
00056   hXF     = new TH1F("hXF","x_{F}=#frac{p_{L}}{#sqrt{s}/2}",50,0.,1.);
00057   hEnergy = new TH1F("hEnergy","Diphoton energy; E [GeV]",50,0.,50.);
00058   hEta    = new TH1F("hEta","Diphoton #eta",24,0.8,2.2);
00059   hPhi    = new TH1F("hPhi","Diphoton #phi",30,-3.142,3.142);
00060   hZvertex= new TH1F("hZvertex","z vertex",150,-150.,150.);
00061   hZgg    = new TH1F("hZgg","energy sharing",50,0.,1.);
00062   hZgg_cut    = new TH1F("hZgg_cut","energy sharing",50,0.,1.);
00063 
00064   hEChi2  = new TH1F("hEChi2","#sum_{points} #frac{(Eu-Ev)^{2}}{Nmip}",50,0.,10.);
00065   hE1Chi2  = new TH1F("hE1Chi2","point1 #frac{(Eu-Ev)^{2}}{Nmip}",50,0.,10.);
00066   hE2Chi2  = new TH1F("hE2Chi2","point2 #frac{(Eu-Ev)^{2}}{Nmip}",50,0.,10.);
00067 
00068   hEChi2_low  = new TH1F("hEChi2_low","#sum_{points} #frac{(Eu-Ev)^{2}}{Nmip}",50,0.,10.);
00069   hE1Chi2_low  = new TH1F("hE1Chi2_low","point1 #frac{(Eu-Ev)^{2}}{Nmip}",50,0.,10.);
00070   hE2Chi2_low  = new TH1F("hE2Chi2_low","point2 #frac{(Eu-Ev)^{2}}{Nmip}",50,0.,10.);
00071 
00072   hEChi2_hi  = new TH1F("hEChi2_hi","#sum_{points} #frac{(Eu-Ev)^{2}}{Nmip}",50,0.,10.);
00073   hE1Chi2_hi  = new TH1F("hE1Chi2_hi","point1 #frac{(Eu-Ev)^{2}}{Nmip}",50,0.,10.);
00074   hE2Chi2_hi  = new TH1F("hE2Chi2_hi","point2 #frac{(Eu-Ev)^{2}}{Nmip}",50,0.,10.);
00075 
00076   hRatio = new TH1F("hRatio","E_{smd} / E_{towers}",20,0.,0.05);
00077   hRatio_low = new TH1F("hRatio_low","E_{smd} / E_{towers}",20,0.,0.05);
00078   hRatio_hi = new TH1F("hRatio_hi","E_{smd} / E_{towers}",20,0.,0.05);
00079   
00080   hEvents = new TH1F("hEvents","event counter",1,0.,1.); hEvents->SetBit(TH1::kCanRebin);
00081 
00082   //hdEds   = new TH2F("hdEds","smd energy derivative; #frac{1}{E}#frac{dE}{ds}; rel. strip",30,-2.0,1.0,10,-5.,5.);
00083 
00084   if ( mFile ) 
00085     {
00086       mTree=new TTree("mPi0Tree","EEmc pi0 events");
00087       mPi0Event=new StEEmcMixEvent();
00088       mTree->Branch("pi0event",&mPi0Event,32000,99);
00089       mTree->SetDirectory( mFile );
00090     }
00091 
00092   return StMaker::Init();
00093 }
00094 
00095 // ---------------------------------------------------------------------------
00096 Int_t StEEmcPi0Maker::Make()
00097 {
00098 
00099   StMuDstMaker *mudst=(StMuDstMaker*)GetMaker("MuDst");
00100   assert(mudst);
00101 
00102   hEvents->Fill("event",1.0);
00103 
00104   if ( !checkTrigger() ) return kStOK;
00105   hEvents->Fill("trigger",1.0);
00106   
00107   StThreeVectorF pv=mudst->muDst()->event()->primaryVertexPosition();
00108   TVector3 vertex( pv.x(), pv.y(), pv.z() );
00109 
00110   if ( pv.z() < -150. || pv.z() > 150.0 || pv.z()==0.0 ) return kStOK;
00111   hEvents->Fill("vertex",1.0);
00112   
00113   // get the points from the point maker
00114   StEEmcPointVec_t points = mEEpoints->points();
00115 
00116   if ( points.size() < 2 ) return kStOK;
00117   hEvents->Fill("two+ points",1.0);
00118 
00119   mPi0Event->setEvent( mudst->muDst()->event() );
00120 
00121   // loop over all pairs of points
00122   for ( UInt_t ipoint=0;ipoint<points.size();ipoint++ )
00123     {
00124       StEEmcPoint p1 = points[ipoint];
00125 
00126 
00127       for ( UInt_t jpoint=0;jpoint<points.size();jpoint++ )
00128         {
00129           if (ipoint>=jpoint) continue;
00130           StEEmcPoint p2=points[jpoint];
00131 
00132 
00133           Int_t phibin1 = p1.tower(0).phibin();
00134           Int_t phibin2 = p2.tower(0).phibin();
00135 
00136           //Int_t etabin1 = p1.tower(0).etabin();
00137           //Int_t etabin2 = p2.tower(0).etabin();
00138 
00139           //Int_t deta = TMath::Abs(etabin1-etabin2);
00140           Int_t dphi = phibin1-phibin2;
00141           if ( dphi < 0 ) dphi+=60;
00142           if ( dphi > 30 ) dphi=60-dphi;
00143 
00144           if ( dphi > 10 ) continue; 
00147           StEEmcPair pair(p1,p2,vertex,vertex);
00148 
00149           
00150           /***********************************************
00151            **
00152            ** Fill the TTre and the list of pi0 pairs
00153            **
00154            ***********************************************/
00155 
00156           mPairs.push_back( pair );
00157           mPi0Event->addPair( pair );
00158 
00159           //Int_t index1 = p1.tower(0).index();
00160           //Int_t index2 = p2.tower(0).index();
00161           
00162           StEEmcClusterVec_t tower_clusters; // tower clusters associated w/ pair (should be 1 or 2)
00163           StEEmcCluster c1 = p1.clusters(0)[0];
00164           StEEmcCluster c2 = p2.clusters(0)[0];
00165           tower_clusters.push_back( c1 );
00166           if ( c1.key() != c2.key() ) {
00167             tower_clusters.push_back(c2);
00168           }
00169 
00170           mPi0Event->mNumberOfTowerClusters.push_back( (Int_t)tower_clusters.size() );
00171 
00172           Int_t ntracks = 0;
00173           Int_t npoints = 0;
00174           for ( UInt_t i=0;i<tower_clusters.size();i++ )
00175             {
00176               ntracks += mEEclusters->numberOfTracks( tower_clusters[i] );
00177               npoints += mEEpoints->numberOfPoints( tower_clusters[i] );
00178             }
00179           mPi0Event->mNumberOfTracks.push_back( ntracks );
00180           mPi0Event->mNumberOfPoints.push_back( npoints );
00181 
00182 
00183 
00184           /***********************************************/
00185           
00186 
00187 
00188 
00190 
00191 
00192           hEvents->Fill("pair",1.0);
00193 
00194           if ( pair.momentum().Perp() < 3.0 ) continue;
00195 
00196           hEvents->Fill("pair pt>3",1.0);
00197 
00198 
00199           //      printf("------------------- retrieve tower clusters from points -------------------\n");
00200           StEEmcCluster cluster1=p1.clusters(0).at(0); //mEEpoints->cluster( p1 );
00201           StEEmcCluster cluster2=p2.clusters(0).at(0); //mEEpoints->cluster( p2 );
00202 
00203           // require both points to match same tower cluster
00204           //555  if ( cluster1.key() != cluster2.key() ) continue;
00205 
00206           hEvents->Fill("pair in cl",1.0);
00207           // require 2 and only 2 smd clusters matched to tower cluster
00208           // in each view
00209 
00210           if ( cluster1.numberOfMatchingClusters(4) > 2 ||
00211                cluster1.numberOfMatchingClusters(5) > 2 ) continue;
00212 
00213           StEEmcSmdCluster u1=p1.cluster(0);
00214           StEEmcSmdCluster u2=p2.cluster(0);
00215           StEEmcSmdCluster v1=p1.cluster(1);
00216           StEEmcSmdCluster v2=p2.cluster(1);
00217 
00218           Bool_t didSplit = false;
00219           if ( u1.key() == u2.key() || v1.key() == v2.key() ) didSplit=true;
00220 
00221 
00222           //      printf("cluster1 key=%i E=%5.2f cluster2 key=%i E=%5.2f\n",
00223           //             cluster1.key(),cluster1.energy(),cluster2.key(),cluster2.energy());
00224 
00225           Float_t esum_smd = u1.energy()+u2.energy()+v1.energy()+v2.energy(); // [GeV]
00226           Float_t esum_tow = cluster1.energy();
00227           if ( cluster1.key() != cluster2.key() ) esum_tow += cluster2.energy();
00228             
00229           
00230 
00231           hEvents->Fill("pair ncl",1.0);
00232 
00233           hMass->Fill( pair.mass(), pair.pt() );
00234           if( didSplit ) hMass_split->Fill( pair.mass(), pair.pt() );
00235 
00236           Float_t ediff1 = 1000.0*(u1.energy()-v1.energy());
00237           Float_t ediff2 = 1000.0*(u2.energy()-v2.energy());
00238           Float_t nmip1  = 1000.0*(u1.energy()+v1.energy())/1.3;
00239           Float_t nmip2  = 1000.0*(u2.energy()+v2.energy())/1.3;
00240           if ( nmip1>0. && nmip2>0. ) 
00241             {
00242               Float_t e1chi2 = ediff1*ediff1/nmip1;
00243               Float_t e2chi2 = ediff2*ediff2/nmip2;
00244               Float_t echi2  = e1chi2+e2chi2;
00245 
00246               if ( e1chi2 < 4.0 && e2chi2 < 4.0 ) {
00247                 hMass_cut -> Fill( pair.mass(), pair.pt() );
00248                 if ( pair.mass()>0.1 && pair.mass()<0.18 ) hPT_cut->Fill( pair.pt() );
00249                 if ( pair.mass()>0.1 && pair.mass()<0.18 ) hZgg_cut->Fill( pair.zgg() );
00250               }
00251 
00252               if ( pair.mass() < 0.1 ) {            
00253                 hEChi2_low->Fill(echi2);
00254                 hE1Chi2_low->Fill(e1chi2);
00255                 hE2Chi2_low->Fill(e2chi2);
00256                 hRatio_low -> Fill( esum_smd/esum_tow );                
00257               }
00258               if ( pair.mass() > 0.2 ) {            
00259                 hEChi2_hi->Fill(echi2);
00260                 hE1Chi2_hi->Fill(e1chi2);
00261                 hE2Chi2_hi->Fill(e2chi2);
00262                 hRatio_hi -> Fill( esum_smd/esum_tow );
00263               }
00264             }
00265 
00266 
00267           if ( pair.mass() > 0.1 && pair.mass() < 0.18 ) {
00268 
00269             hPT->Fill( pair.momentum().Perp() );
00270             hXF->Fill( pair.momentum().Z()/100.0 );
00271             hEnergy->Fill( pair.energy() );
00272             hEta->Fill( pair.momentum().Eta() );
00273             hPhi->Fill( pair.momentum().Phi() );
00274             hZvertex->Fill( pair.vertex().Z() );
00275             hZgg->Fill( pair.zgg() );
00276 
00277             if ( nmip1==0. || nmip2==0. ) continue;
00278             Float_t e1chi2 = ediff1*ediff1/nmip1;
00279             Float_t e2chi2 = ediff2*ediff2/nmip2;
00280             Float_t echi2  = e1chi2+e2chi2;
00281  
00282             hEChi2->Fill(echi2);
00283             hE1Chi2->Fill(e1chi2);
00284             hE2Chi2->Fill(e2chi2);
00285             hRatio -> Fill( esum_smd/esum_tow );
00286             printf("esum_smd/esum_tow = %5.3f\n", esum_smd/esum_tow );
00287             
00288           }
00289 
00290         }
00291     }
00292 
00293   if ( mTree )
00294     mTree->Fill();
00295   
00296   return kStOK;
00297 }
00298 
00299 
00300 
00301 // ---------------------------------------------------------------------------
00302 void StEEmcPi0Maker::Clear(Option_t *opts)
00303 {
00304   mPairs.clear();
00305   mPi0Event->Clear();
00306 }
00307 
00308 
00309 
00310 // ----------------------------------------------------------------------------                                                     
00311 void  StEEmcPi0Maker::setCheckTrigger(Bool_t t){ mCheckTrigger=t; }
00312 void  StEEmcPi0Maker::addTrigger(Int_t t){ mTriggerList.push_back(t); }
00313 
00314 Bool_t StEEmcPi0Maker::checkTrigger()
00315 {
00316 
00318   if ( mTriggerList.size() == 0 ) return 1;
00319 
00320   StTriggerId nominal;
00321 
00323   StMuDstMaker *mumk = (StMuDstMaker*)GetMaker("MuDst");
00324   StEvent *event = (StEvent*)GetInputDS("StEvent");
00325 
00326   if ( mumk )
00327     {
00328       nominal = mumk->muDst()->event()->triggerIdCollection().nominal();
00329       goto CHECK_TRIGGER;
00330     }
00331 
00333 
00334   if ( event )
00335     {
00336       nominal=*event->triggerIdCollection()->nominal();
00337       goto CHECK_TRIGGER;
00338     }
00339 
00341   goto NO_DATA;
00342 
00343 
00344  CHECK_TRIGGER:
00345 
00346   for ( UInt_t ii=0;ii<mTriggerList.size();ii++ )
00347     {
00348       if ( nominal.isTrigger( mTriggerList[ii] ) ) return mTriggerList[ii];
00349     }
00350   return 0;
00351 
00352  NO_DATA:
00353   assert(2+2==5); // noooo data                                                                                                    $
00354   return 0;
00355 
00356 }