StRoot  1
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Groups Pages
StFmsDiPi0.cxx
1 // \class StFmsDiPi0
2 // \author Akio Ogawa
3 //
4 // $Id: StFmsDiPi0.cxx,v 1.5 2017/09/05 17:42:19 akio Exp $
5 // $Log: StFmsDiPi0.cxx,v $
6 // Revision 1.5 2017/09/05 17:42:19 akio
7 // update
8 //
9 // Revision 1.4 2016/11/17 18:57:23 akio
10 // Many updates
11 //
12 // Revision 1.3 2016/10/10 19:17:40 akio
13 // *** empty log message ***
14 //
15 // Revision 1.2 2016/06/08 16:28:09 akio
16 // *** empty log message ***
17 //
18 // Revision 1.1 2016/01/20 19:58:55 akio
19 // *** empty log message ***
20 //
21 // Revision 1.1 2016/01/20 19:50:04 akio
22 // *** empty log message ***
23 //
24 // Revision 1.1 2015/10/20 19:55:51 akio
25 // Initial version of FMS di-pi0 analysis
26 //
27 
28 #include "StFmsDiPi0.h"
29 
30 #include "StMessMgr.h"
31 #include "Stypes.h"
32 #include "StMuDSTMaker/COMMON/StMuTypes.hh"
33 
34 #include "StThreeVectorF.hh"
35 #include "StFmsDbMaker/StFmsDbMaker.h"
36 #include "StEnumerations.h"
37 //#include "StEventTypes.h"
38 #include "StEvent/StEvent.h"
39 #include "StEvent/StFmsCollection.h"
40 #include "StEvent/StFmsHit.h"
41 #include "StEvent/StFmsPoint.h"
42 #include "StEvent/StFmsPointPair.h"
43 #include "StEvent/StTriggerData.h"
44 
45 #include "StSpinPool/StFmsFastSimMaker/StFmsFastSimMaker.h"
46 //#include "StarGenerator/BASE/StarPrimaryMaker.h"
47 //#include "StarGenerator/EVENT/StarGenEvent.h"
48 //#include "StarGenerator/EVENT/StarGenParticle.h"
49 
50 #include "TFile.h"
51 #include "TH1F.h"
52 #include "TH2F.h"
53 
54 static const double PI=TMath::Pi();
55 static const double twoPI=PI*2.0;
56 
57 const float mPtBin[7]={0.5,1.0,1.5,2.0,2.5,3.0,4.0};
58 const float energyCut=1.0;
59 const float ptCut=0.1;
60 const float ZggCut=0.7;
61 const float MassCut0=0.07;
62 const float MassCut1=0.2;
63 const float MassCut2=0.35;
64 
65 static const int NTRG=12; //123=FMS-sm-bs123,456=FMS-lg-bs123,7=FMS-DiBS,8910=FMS-JP012,11=FMS-DiJP,13=LED
66 static const char *CTRG[NTRG] = {"SmBS1","SmBS2","SmBS3","LgBS1",
67  "LgBS2","LgBS3","DiBS" ,"JP2",
68  "JP1" ,"JP0" ,"DiJP" ,"LED"};
69 
70 static const short MAXP=256;
71 unsigned int tNEVT;
72 unsigned short tTRG,tBBCE,tZDCE;
73 unsigned char tBC, tBBCM, tTOFM, tNP;
74 unsigned char tDet[MAXP];
75 unsigned short tId[MAXP];
76 float tX[MAXP], tY[MAXP];
77 float tPx[MAXP], tPy[MAXP], tPz[MAXP], tE[MAXP];
78 
79 int getFmsTrigId(const StTriggerId& trgid, int print=0){
80  const int TIDBASE=480800;
81  const int NBEAM=4;
82  const int MAXVERSION=3;
83  int trig=0;
84  if(print) LOG_INFO << "TRGID = ";
85  for(int k=0; k<NBEAM; k++){
86  for(int j=0; j<MAXVERSION; j++){
87  for(int i=1; i<=NTRG; i++){
88  int l=i;
89  if(i==12) l=13;
90  int id=TIDBASE + 10000*k + 20*j + l;
91  if(trgid.isTrigger(id)){
92  trig |= (1<<(i-1));
93  if(print)LOG_INFO << CTRG[i-1] << " ";
94  }
95  }
96  }
97  }
98  if(print) LOG_INFO << endm;
99  return trig;
100 }
101 
102 double wrapAround(double phi){
103  double res=phi;
104  while(res>=1.5*PI) {res-=twoPI;}
105  while(res<-0.5*PI) {res+=twoPI;}
106  return res;
107 }
108 
109 ClassImp(StFmsDiPi0);
110 
111 StFmsDiPi0::StFmsDiPi0(const Char_t* name):
112  StMaker(name),mFilename((char *)"fmsDiPi0.root") {}
113 
114 StFmsDiPi0::~StFmsDiPi0(){}
115 
116 Int_t StFmsDiPi0::Init(){
117  if(!mReadTree && !mPythia){
118  mFmsDbMaker=static_cast<StFmsDbMaker*>(GetMaker("fmsDb"));
119  if(!mFmsDbMaker){
120  LOG_ERROR << "StFmsDiPi0::InitRun Failed to get StFmsDbMaker" << endm;
121  return kStFatal;
122  }
123  }
124 
125  if(!mWriteTree){
126  mFile=new TFile(mFilename,"RECREATE");
127  char c[100];
128  mBC=new TH1F("BC","BC",120,0.0,120.0);
129  mBBC=new TH1F("BBCE","BBCE",100,0.0,70000.0);
130  mBBCAG=new TH1F("BBCEAG","BBCEAG",100,0.0,70000.0);
131  mBBCM=new TH1F("BBCMult","BBCMult",17,0.0,17.0);
132  mBBCMAG=new TH1F("BBCMultAG","BBCMultAG",17,0.0,17.0);
133  mTOF=new TH1F("TOF","TOF",150,0.0,300.0);
134  mTOFAG=new TH1F("TOFAG","TOFAG",150,0.0,300.0);
135  mBBCTOF=new TH2F("BBC_TOF","BBC_TOF",50,0.0,70000.0,50,0.0,300.0);
136  mBBCMTOF=new TH2F("BBCMult_TOF","BBCMult_TOF",17,0.0,17.0,50,0.0,300.0);
137  mBBCBBCM=new TH2F("BBC_BBCMult","BBC_BBCMult",50,0.0,70000.0,17,0.0,17.0);
138  mTOFTOF=new TH2F("TOF_TOF","TOF_TOF",50,0.0,300.0,50,0.0,300.0);
139 
140  mMass0=new TH2F("Mass0","Mass0",50,0.0, 0.4,16,0.0,16.0);
141  mMass1=new TH2F("Mass1","Mass1",50,0.0, 0.4,16,0.0,16.0);
142  mMass2=new TH2F("Mass2","Mass2",50,0.0, 0.4,16,0.0,16.0);
143  mEne =new TH2F("Ene", "Ene", 50,0.0,100.0,16,0.0,16.0);
144  mPt =new TH2F("pT", "pT", 50,0.0, 10.0,16,0.0,16.0);
145 
146  for(int i=0; i<kNPtBin; i++){
147  for(int j=0; j<=i; j++){
148  for(int k=0; k<=kNCut; k++){
149  if(j==0){
150  sprintf(c,"m0_%1d_c%d",i,k); mM0[i][k]=new TH1F(c,c,50,0.0,0.4);
151  sprintf(c,"phi0_%1d_c%d",i,k); mPhi0[i][k]=new TH1F(c,c,128,-0.5*PI,1.5*PI);
152  sprintf(c,"etaphi0_%1d_c%d",i,k); mEtaPhi0[i][k]=new TH2F(c,c,50,2.5,4.5,64,-0.5*PI,1.5*PI);
153  }
154  sprintf(c,"m1_%1d%1d_c%d",i,j,k); mM1[i][j][k]=new TH1F(c,c,50,0.0,0.4);
155  sprintf(c,"m2_%1d%1d_c%d",i,j,k); mM2[i][j][k]=new TH1F(c,c,50,0.0,0.4);
156  sprintf(c,"z1_%1d%1d_c%d",i,j,k); mZ1[i][j][k]=new TH1F(c,c,50,0.0,1.0);
157  sprintf(c,"z2_%1d%1d_c%d",i,j,k); mZ2[i][j][k]=new TH1F(c,c,50,0.0,1.0);
158  sprintf(c,"e1_%1d%1d_c%d",i,j,k); mE1[i][j][k]=new TH1F(c,c,50,0.0,100.0);
159  sprintf(c,"e2_%1d%1d_c%d",i,j,k); mE2[i][j][k]=new TH1F(c,c,50,0.0,100.0);
160  sprintf(c,"pt1_%1d%1d_c%d",i,j,k); mPt1[i][j][k]=new TH1F(c,c,50,0.0,10.0);
161  sprintf(c,"pt2_%1d%1d_c%d",i,j,k); mPt2[i][j][k]=new TH1F(c,c,50,0.0,10.0);
162  sprintf(c,"eta1_%1d%1d_c%d",i,j,k); mEta1[i][j][k]=new TH1F(c,c,50,2.0,5.0);
163  sprintf(c,"eta2_%1d%1d_c%d",i,j,k); mEta2[i][j][k]=new TH1F(c,c,50,2.0,5.0);
164  sprintf(c,"phi1_%1d%1d_c%d",i,j,k); mPhi1[i][j][k]=new TH1F(c,c,128,-0.5*PI,1.5*PI);
165  sprintf(c,"phi2_%1d%1d_c%d",i,j,k); mPhi2[i][j][k]=new TH1F(c,c,128,-0.5*PI,1.5*PI);
166  sprintf(c,"dphi_%1d%1d_c%d",i,j,k); mDphi[i][j][k]=new TH1F(c,c,128,-0.5*PI,1.5*PI);
167  sprintf(c,"bbce_%1d%1d_c%d",i,j,k); mBbce[i][j][k]=new TH1F(c,c,50,0.0,70000.0);
168  sprintf(c,"tofm_%1d%1d_c%d",i,j,k); mTofm[i][j][k]=new TH1F(c,c,50,0.0,250.0);
169  sprintf(c,"phi1dphi_%1d%1d_c%d",i,j,k); mPhi1Dphi[i][j][k]=new TH2F(c,c,50,-0.5*PI,1.5*PI,50,-0.5*PI,1.5*PI);
170  }
171  }
172  }
173  }
174 
175  if(mWriteTree){
176  mTreeFile=new TFile(mTreeFilename,"recreate");
177  mTree=new TTree("dipi0","dipi0 tree");
178  mTree->Branch("trg",&tTRG,"trg/s");
179  mTree->Branch("bbce",&tBBCE,"bbce/s");
180  mTree->Branch("zdce",&tZDCE,"zdce/s");
181  mTree->Branch("bc",&tBC,"bc/b");
182  mTree->Branch("bbcm",&tBBCM,"bbcm/b");
183  mTree->Branch("tofm",&tTOFM,"tofm/b");
184  mTree->Branch("np",&tNP,"np/b");
185  mTree->Branch("det",tDet,"det[np]/b");
186  mTree->Branch("id",tId,"id[np]/s");
187  mTree->Branch("e",tE,"e[np]/F");
188  mTree->Branch("x",tX,"x[np]/F");
189  mTree->Branch("y",tY,"y[np]/F");
190  mTree->Branch("px",tPx,"px[np]/F");
191  mTree->Branch("py",tPy,"py[np]/F");
192  mTree->Branch("pz",tPz,"pz[np]/F");
193  }
194 
195  if(mReadTree){
196  //mTreeFile=new TFile(mTreeFilename);
197  //mTree = (TTree*)mTreeFile->Get("dipi0");
198  //tNEVT=mChain->GetEntries();
199  //LOG_INFO << Form("Found %d entries in the input tree",tNEVT) << endm;
200  mChain->SetBranchAddress("trg",&tTRG);
201  mChain->SetBranchAddress("bbce",&tBBCE);
202  mChain->SetBranchAddress("zdce",&tZDCE);
203  mChain->SetBranchAddress("bc",&tBC);
204  mChain->SetBranchAddress("bbcm",&tBBCM);
205  mChain->SetBranchAddress("tofm",&tTOFM);
206  mChain->SetBranchAddress("np",&tNP);
207  mChain->SetBranchAddress("det",tDet);
208  mChain->SetBranchAddress("id",tId);
209  mChain->SetBranchAddress("e",tE);
210  mChain->SetBranchAddress("x",tX);
211  mChain->SetBranchAddress("y",tY);
212  mChain->SetBranchAddress("px",tPx);
213  mChain->SetBranchAddress("py",tPy);
214  mChain->SetBranchAddress("pz",tPz);
215  mFile->cd();
216  }
217  return kStOK;
218 }
219 
220 Int_t StFmsDiPi0::Finish(){
221  if(!mWriteTree){
222  LOG_INFO << Form("Writing and Closing %s",mFilename) << endm;
223  mFile->Write();
224  mFile->Close();
225  }else{
226  LOG_INFO << Form("Writing and Closing %s",mTreeFilename) << endm;
227  mTreeFile->cd();
228  mTreeFile->Write();
229  mTreeFile->Close();
230  }
231  return kStOK;
232 }
233 
234 Int_t StFmsDiPi0::ptbin(float pt){
235  if(pt < mPtBin[0]) return -1;
236  for(int i=0; i<kNPtBin; i++){
237  if(pt >= mPtBin[i] && pt < mPtBin[i+1]) return i;
238  }
239  return kNPtBin;
240 }
241 
242 Int_t StFmsDiPi0::Make(){
243  int print=0;
244  int ag=0, bbcEsum=0, zdce=0, bbcmult=0, tofmul1=0, tofmul2=0, bsTrg=0, bsSTrg=0, bsLTrg=0, jetTrg=0, doubleTrg=0;
245  int bunch=-1, trigger=-1;
246 
247  //LOG_INFO << "BBCCuts = " << mBBCCut1 <<" " << mBBCCut2 <<" "<< mBBCCut3 <<" "<< mBBCCut4 <<" "<<endm;
248  if(mPythia==1){
249  bsTrg=1;
250  }else{
251  if(mReadTree==1){
252  static unsigned int ievt=0;
253  //if(ievt>tNEVT) return kStEOF;
254  //mTreeFile->cd();
255  int ret=mChain->GetEntry(ievt++);
256  if(ret==0) return kStEOF;
257  bunch=tBC;
258  bbcEsum=tBBCE;
259  zdce=tZDCE;
260  bbcmult=tBBCM;
261  tofmul2=tTOFM;
262  trigger=tTRG;
263  mFile->cd();
264  }else{ //Reading from MuDst
265  StMuDst* mudst = (StMuDst*)GetInputDS("MuDst");
266  if(!mudst) {LOG_ERROR << "StFmsDiPi0::Make did not find MuDst"<<endm; return kStErr;}
267  StMuEvent* mueve = mudst->event();
268  if(!mueve) {LOG_ERROR << "StFmsDiPi0::Make did not find MuEvent"<<endm; return kStErr;}
269  const StTriggerData* trgd=mueve->triggerData();
270  if(!trgd) {LOG_ERROR << "StFmsDiPi0::Make did not find StTriggerData from MuDst"<<endm; return kStErr;}
271  //bunch and abort gap
272  bunch = trgd->bunchId7Bit();
273  //bbc east sum
274  for(int i=1; i<=16; i++){
275  int tac=trgd->bbcTDC(east,i);
276  if(tac>100 && tac<2400){
277  int adc=trgd->bbcADC(east,i);
278  if(adc>50) bbcmult++;
279  bbcEsum += trgd->bbcADC(east,i);
280  }
281  }
282  zdce=trgd->zdcAttenuated(east);
283  tofmul1=trgd->tofMultiplicity();
284  tofmul2=mudst->event()->btofTrayMultiplicity();
285  //triggers
286  trigger = getFmsTrigId(mudst->event()->triggerIdCollection().nominal());
287  if(print)
288  LOG_INFO << Form("Run=%8d Event=%10d",mueve->runNumber(),mueve->eventNumber())<<endm;
289  }
290  if (bunch>=30 && bunch< 40) {ag=1;}
291  else if(bunch>=110 && bunch<120) {ag=2;}
292  else if(bbcmult==0 || tofmul2<3) {ag=3;}
293  if(trigger & 0x03f) bsTrg=1;
294  if(trigger & 0x007) bsSTrg=1;
295  if(trigger & 0x038) bsLTrg=1;
296  if(trigger & 0x380) jetTrg=1;
297  if(trigger & 0x440) doubleTrg=1;
298 
299  if(!mWriteTree){
300  mBC->Fill(float(bunch));
301  if(ag==0 || ag==3){
302  mBBC->Fill(float(bbcEsum));
303  mBBCM->Fill(float(bbcmult));
304  mTOF->Fill(float(tofmul2));
305  }else if(ag==1){
306  mBBCAG->Fill(float(bbcEsum));
307  mBBCMAG->Fill(float(bbcmult));
308  mTOFAG->Fill(float(tofmul2));
309  }
310  mBBCTOF->Fill(bbcEsum,tofmul2);
311  mBBCMTOF->Fill(bbcmult,tofmul2);
312  mBBCBBCM->Fill(bbcEsum,bbcmult);
313  mTOFTOF->Fill(tofmul1,tofmul2);
314  }
315  if(print)
316  LOG_INFO << Form("Bunch=%3d AGap=%1d BbcEsum=%6d Tof=%6d | %6d",
317  bunch,ag,bbcEsum,tofmul1,tofmul2) << endm;
318  }
319 
320  int np=0;
321  StPtrVecFmsPoint pcut;
322  if(mReadTree>0){
323  np=tNP;
324  for(int i=0; i<np; i++){
325  StFmsPoint* point = new StFmsPoint();
326  point->setDetectorId(tDet[i]);
327  point->setId(tId[i]);
328  point->setEnergy(tE[i]);
329  point->setX(tX[i]);
330  point->setX(tY[i]);
331  point->setFourMomentum(StLorentzVectorF(tPx[i],tPy[i],tPz[i],tE[i]));
332  pcut.push_back(point);
333  }
334  }else{
335  StEvent* event = (StEvent*)GetInputDS("StEvent");
336  if(!event) {LOG_ERROR << "StFmsDiPi0::Make did not find StEvent"<<endm; return kStErr;}
337  mFmsColl = event->fmsCollection();
338  if(!mFmsColl) {LOG_ERROR << "StFmsDiPi0::Make did not find StEvent->FmsCollection"<<endm; return kStErr;}
339  StSPtrVecFmsPoint &points = mFmsColl->points();
340  np=mFmsColl->numberOfPoints();
341  if(print) cout << Form("Npoint = %d",np) << endl;
342  //FV cut and minimal energy/pT cut for points
343  for(int i=0; i<np; i++){
344  int edgeType;
345  float distance=mFmsDbMaker->distanceFromEdge(points[i],edgeType);
346  if(mPythia || distance<-0.51 || (mFmsColl->isMergeSmallToLarge() && edgeType==4)){ //FV cut
347  if(mPythia || (points[i]->energy()>energyCut && points[i]->fourMomentum().perp()>ptCut)){ //minimal E and pT cuts
348  pcut.push_back(points[i]);
349  }
350  }
351  }
352  }
353 
354  if(print) cout << Form("Ncut = %d",pcut.size()) << endl;
355  if(pcut.size()>1){
356  //sort by pT
357  std::sort(pcut.begin(), pcut.end(), [](StFmsPoint* a, StFmsPoint* b) {
358  return b->fourMomentum().perp() < a->fourMomentum().perp();
359  });
360  //pair
361  vector<StFmsPointPair*> pair;
362  std::vector<int> pcutuse(pcut.size(), 0);
363  for(unsigned int i=0; i<pcut.size()-1; i++){
364  for(unsigned int j=i+1; j<pcut.size(); j++){
365  StFmsPointPair* pp = new StFmsPointPair(pcut[i],pcut[j]);
366  if(mPythia || (pp->zgg() < ZggCut && pp->pT() > mPtBin[0])){
367  pair.push_back(pp);
368  pcutuse[i]=1;
369  pcutuse[j]=1;
370  }else{
371  delete pp;
372  }
373  }
374  }
375  if(print) cout << Form("Npair = %d",pair.size()) << endl;
376  if(pair.size()>0) {
377 
378  // write Tree file
379  if(mWriteTree){
380  mTreeFile->cd();
381  tBC=bunch;
382  tTRG=trigger;
383  tBBCM=bbcmult;
384  tTOFM=tofmul2;
385  tBBCE=bbcEsum;
386  tZDCE=zdce;
387  tNP=pcut.size();
388  for(unsigned int i=0; i<pcut.size(); i++){
389  if(pcutuse[i]==1){
390  tDet[i]=pcut[i]->detectorId();
391  tId[i]=pcut[i]->id();
392  tE[i]=pcut[i]->energy();;
393  tX[i]=pcut[i]->x();
394  tY[i]=pcut[i]->y();
395  tPx[i]=pcut[i]->fourMomentum().px();
396  tPy[i]=pcut[i]->fourMomentum().py();
397  tPz[i]=pcut[i]->fourMomentum().pz();
398  }
399  }
400  mTree->Fill();
401  return kStOK; //exit if writing TTree
402  }
403 
404  //sort by pT of pair
405  std::sort(pair.begin(), pair.end(), [](StFmsPointPair* a, StFmsPointPair* b) {
406  return b->pT() < a->pT();
407  });
408  if(print){
409  for(unsigned int i=0; i<pair.size(); i++){
410  LOG_INFO << Form("%3d m=%6.3f e=%6.2f pt=%6.2f eta=%6.2f phi=%6.3f zgg=%6.3f id=%4d %4d",
411  i,
412  pair[i]->mass(), pair[i]->energy(),pair[i]->pT(),
413  pair[i]->eta(), pair[i]->phi(), pair[i]->zgg(),
414  pair[i]->point(0)->id(), pair[i]->point(1)->id() ) << endm;
415  }
416  }
417  //pair of pair (4 photons or 2 pairs)
418  int ndipi0=0;
419  int firstpair=0;
420  int used1[2000],used2[2000],used3[2000],used4[2000];
421  memset(used1,0,sizeof(used1)); //pi0-pi0
422  memset(used2,0,sizeof(used2)); //pi0-highmass
423  memset(used3,0,sizeof(used3)); //highmass-pi0
424  memset(used4,0,sizeof(used4)); //highmass-highmass
425 
426  for(unsigned int i=0; i<pair.size(); i++){ //loop over 1st pair
427  int ptbin0=ptbin(pair[i]->pT());
428  if(ptbin0<0) continue; //below lowest pt bin
429 
430  int id0=pair[i]->point(0)->id();
431  int id1=pair[i]->point(1)->id();
432  float m0=pair[i]->mass();
433  float z0=pair[i]->zgg();
434  float p0=wrapAround(pair[i]->phi());
435  int topbot=1; //bottom
436  if(p0>=0.0 && p0<PI) topbot=0;//top
437 
438  //QA to look for mass/pt/energy for each octet
439  float pp=p0;
440  while(pp<0.0) pp+=twoPI;
441  while(pp>=twoPI) pp-=twoPI;
442  int ii=int(8.0*pp/twoPI);
443  int dd=pair[i]->point(0)->detectorId();
444  if(dd==kFmsNorthSmallDetId || dd==kFmsSouthSmallDetId) ii+=8;
445  mMass0->Fill(m0,float(ii));
446  if(pair[i]->pT()>1.5) mMass1->Fill(m0,float(ii));
447  if(pair[i]->energy()>30.0) mMass2->Fill(m0,float(ii));
448  if(m0>=MassCut0 && m0<MassCut1){
449  mEne->Fill(pair[i]->energy(),float(ii));
450  mPt->Fill(pair[i]->pT(),float(ii));
451  }
452 
453  if(ptbin0>=kNPtBin) continue; //above highest pt bin
454  int cut1[kNCut];
455  memset(cut1,0,sizeof(cut1));
456 
457  //pi0 mass and highmass
458  int masscut1=0;
459  if(m0>=MassCut0 && m0<MassCut1) {masscut1=1;}
460  else if(m0>=MassCut1 && m0<MassCut2) {masscut1=2;}
461 
462  //4 exclusiveness for pi0-pi0, pi0-hightmass, highmass-pi0, highmass-highmass
463  int exclusive1=0, exclusive2=0, exclusive3=0, exclusive4=0;
464  if(masscut1==1 && ((used1[id0]==0 && used1[id1]==0) || (used1[id0]==id1 && used1[id1]==id0))) {exclusive1=1; used1[id0]=id1; used1[id1]=id0;}
465  if(masscut1==1 && ((used2[id0]==0 && used2[id1]==0) || (used2[id0]==id1 && used2[id1]==id0))) {exclusive2=1; used2[id0]=id1; used2[id1]=id0;}
466  if(masscut1==2 && ((used3[id0]==0 && used3[id1]==0) || (used3[id0]==id1 && used3[id1]==id0))) {exclusive3=1; used3[id0]=id1; used3[id1]=id0;}
467  if(masscut1==2 && ((used4[id0]==0 && used4[id1]==0) || (used4[id0]==id1 && used4[id1]==id0))) {exclusive4=1; used4[id0]=id1; used4[id1]=id0;}
468 
469  //setup cuts
470  if(ag==0 && bsTrg==1){
471  cut1[0]=1;
472  if(masscut1==1){
473  cut1[1]=1;
474  if(exclusive1==1){
475  cut1[2]=1;
476  if (bbcEsum<mBBCCut1) {cut1[3]=1;}
477  else if(bbcEsum<mBBCCut2) {cut1[4]=1;}
478  else if(bbcEsum<mBBCCut3) {cut1[5]=1;}
479  else if(bbcEsum<mBBCCut4) {cut1[15]=1;}
480  else {cut1[16]=1;}
481  }
482  }
483  }
484  if(ag==1 && bsTrg==1 && exclusive1) cut1[6]=1;
485  if(ag==3 && bsTrg==1 && exclusive1) cut1[17]=1;
486  if(ag==0 && jetTrg==1 && exclusive1) cut1[7]=1;
487  if(ag==0 && doubleTrg==1 && exclusive1) cut1[8]=1;
488  if(ag==0 && bsTrg==1 && exclusive2) cut1[9]=1; //this is for pi0-highmass pair
489  if(ag==0 && bsTrg==1 && exclusive3) cut1[10]=1; //this is for highmass-pi0 pair
490  if(ag==0 && bsTrg==1 && exclusive4) cut1[11]=1; //this is for highmass-highmass pair
491  if(ag==0 && bsSTrg==1 && exclusive1) cut1[12]=1;
492  if(ag==0 && bsLTrg==1 && exclusive1) cut1[13]=1;
493  if(ag==0 && bsTrg==1 && exclusive1 && firstpair==0) cut1[14]=1;
494  if(ag==0 && bsTrg==1 && exclusive1 && topbot==0) cut1[18]=1;
495  if(ag==0 && bsTrg==1 && exclusive1 && topbot==1) cut1[19]=1;
496 
497  for(unsigned int k=0; k<kNCut; k++){
498  if(cut1[k]==1){
499  mM0[ptbin0][k]->Fill(m0); //this is for normalization
500  mPhi0[ptbin0][k]->Fill(p0); //this is for mix event
501  mEtaPhi0[ptbin0][k]->Fill(pair[i]->eta(),p0);
502  }
503  }
504  //loop over 2nd pair
505  for(unsigned int j=i+1; j<pair.size(); j++){
506  if(i==j) continue; //same pair
507  int ptbin1=ptbin(pair[j]->pT());
508  if(ptbin1<0) continue; //below lowest pt bin
509  int id2=pair[j]->point(0)->id();
510  int id3=pair[j]->point(1)->id();
511  if(id0==id2 || id0==id3 || id1==id2 || id1==id3) continue; //using same photon
512  float dphi = wrapAround(pair[i]->phi() - pair[j]->phi());
513  float m1=pair[j]->mass();
514  float z1=pair[j]->zgg();
515  float p1=wrapAround(pair[j]->phi());
516  int cut2[kNCut];
517  memset(cut2,0,sizeof(cut2));
518  int masscut2=0;
519  if (m1>=MassCut0 && m1<MassCut1) {masscut2=1;}
520  else if(m1>=MassCut1 && m1<MassCut2) {masscut2=2;}
521 
522  int exc1=0,exc2=0, exc3=0, exc4=0;
523  if(masscut2==1 && exclusive1 &&
524  ((used1[id2]==0 && used1[id3]==0) || (used1[id2]==id3 && used1[id3]==id2)) ){
525  used1[id2]=id3; used1[id3]=id2;
526  exc1=1;
527  }
528  if(masscut2==2 && exclusive2 &&
529  ((used2[id2]==0 && used2[id3]==0) || (used2[id2]==id3 && used2[id3]==id2)) ){
530  used1[id2]=id3; used2[id3]=id2;
531  exc2=1;
532  }
533  if(masscut2==1 && exclusive3 &&
534  ((used3[id2]==0 && used3[id3]==0) || (used3[id2]==id3 && used3[id3]==id2)) ){
535  used3[id2]=id3; used3[id3]=id2;
536  exc3=1;
537  }
538  if(masscut2==2 && exclusive4 &&
539  ((used4[id2]==0 && used4[id3]==0) || (used4[id2]==id3 && used4[id3]==id2)) ){
540  used4[id2]=id3; used4[id3]=id2;
541  exc4=1;
542  }
543 
544  if(ag==0 && bsTrg==1){
545  cut2[0]=1;
546  if(masscut2==1){
547  cut2[1]=1;
548  if(exc1==1){
549  cut2[2]=1;
550  if (bbcEsum<mBBCCut1) {cut2[3]=1;}
551  else if(bbcEsum<mBBCCut2) {cut2[4]=1;}
552  else if(bbcEsum<mBBCCut3) {cut2[5]=1;}
553  else if(bbcEsum<mBBCCut4) {cut2[15]=1;}
554  else {cut2[16]=1;}
555  }
556  }
557  }
558  if(ag==1 && bsTrg==1 && exc1) cut2[6]=1;
559  if(ag==3 && bsTrg==1 && exc1) cut2[17]=1;
560  if(ag==0 && jetTrg==1 && exc1) cut2[7]=1;
561  if(ag==0 && doubleTrg==1 && exc1) cut2[8]=1;
562  if(ag==0 && bsTrg==1 && exc2) cut2[9]=1; //this is for pi0-highmass pair
563  if(ag==0 && bsTrg==1 && exc3) cut2[10]=1; //this is for highmass-pi0 pair
564  if(ag==0 && bsTrg==1 && exc4) cut2[11]=1; //this is for highmass-highmass pair
565  if(ag==0 && bsSTrg==1 && exc1) cut2[12]=1;
566  if(ag==0 && bsLTrg==1 && exc1) cut2[13]=1;
567  if(ag==0 && bsTrg==1 && exc1 && firstpair==0) {firstpair=1;cut2[14]=1;}
568  if(ag==0 && bsTrg==1 && exc1 && topbot==0) cut2[18]=1;
569  if(ag==0 && bsTrg==1 && exc1 && topbot==1) cut2[19]=1;
570 
571  //filling histos
572  for(unsigned int k=0; k<kNCut; k++){
573  if(cut2[k]==1){
574  mM1[ptbin0][ptbin1][k]->Fill(m0);
575  mM2[ptbin0][ptbin1][k]->Fill(m1);
576  mZ1[ptbin0][ptbin1][k]->Fill(z0);
577  mZ2[ptbin0][ptbin1][k]->Fill(z1);
578  mE1[ptbin0][ptbin1][k]->Fill(pair[i]->energy());
579  mE2[ptbin0][ptbin1][k]->Fill(pair[j]->energy());
580  mPt1[ptbin0][ptbin1][k]->Fill(pair[i]->pT());
581  mPt2[ptbin0][ptbin1][k]->Fill(pair[j]->pT());
582  mEta1[ptbin0][ptbin1][k]->Fill(pair[i]->eta());
583  mEta2[ptbin0][ptbin1][k]->Fill(pair[j]->eta());
584  mPhi1[ptbin0][ptbin1][k]->Fill(p0);
585  mPhi2[ptbin0][ptbin1][k]->Fill(p1);
586  mDphi[ptbin0][ptbin1][k]->Fill(dphi);
587  mBbce[ptbin0][ptbin1][k]->Fill(bbcEsum);
588  mPhi1Dphi[ptbin0][ptbin1][k]->Fill(p0,dphi);
589  }
590  }
591  if(print) LOG_INFO << Form("diPair=%3d pair=%2d %2d dphi=%6.3f ptbin=%1d %1d masscut=%1d ag=%1d exc=%1d",
592  ndipi0,i,j,dphi,ptbin0,ptbin1,masscut2,ag,exclusive2)<<endm;
593  ndipi0++;
594  }
595  }
596  }
597  for(unsigned int i=0; i<pair.size(); i++) delete pair[i];
598  }
599  if(mReadTree){
600  for(unsigned int i=0; i<pcut.size(); i++) delete pcut[i];
601  }
602 
603  //Pythia events, so we know real pi0
604  if(mPythia){
605  StFmsFastSimMaker* fmssim = (StFmsFastSimMaker*)GetMaker("FmsFastSimMaker");
606  if(!fmssim){
607  LOG_ERROR << "Did not find FmsFastSimMaker!!!" << endm;
608  return kStErr;
609  }
610  int n=fmssim->nPi0();
611  if(n>1){
612  for(int i=0; i<n; i++){
613  int ptbin0=ptbin(fmssim->pi0(i)->perp());
614  if(ptbin0<0) continue; //below lowest pt bin
615  mM0[ptbin0][kNCut]->Fill(0.135); //this is for normalization
616  if(print)
617  cout << Form("RealPi0 n=%2d pt=%6.2f eta=%6.2f phi=%6.3f",
618  n,fmssim->pi0(i)->perp(),
619  fmssim->pi0(i)->pseudoRapidity(),fmssim->pi0(i)->phi())<<endl;
620  for(int j=i+1; j<n; j++){
621  if(i==j) continue;
622  int ptbin1=ptbin(fmssim->pi0(j)->perp());
623  if(ptbin1<0) continue; //below lowest pt bin
624  float phi1=wrapAround(fmssim->pi0(i)->phi());
625  float phi2=wrapAround(fmssim->pi0(j)->phi());
626  float dphi=wrapAround(phi1-phi2);
627  mM1[ptbin0][ptbin1][kNCut]->Fill(0.135);
628  mM2[ptbin0][ptbin1][kNCut]->Fill(0.135);
629  mPt1[ptbin0][ptbin1][kNCut]->Fill(fmssim->pi0(i)->perp());
630  mPt2[ptbin0][ptbin1][kNCut]->Fill(fmssim->pi0(j)->perp());
631  mE1[ptbin0][ptbin1][kNCut]->Fill(fmssim->pi0(i)->e());
632  mE2[ptbin0][ptbin1][kNCut]->Fill(fmssim->pi0(j)->e());
633  mEta1[ptbin0][ptbin1][kNCut]->Fill(fmssim->pi0(i)->pseudoRapidity());
634  mEta2[ptbin0][ptbin1][kNCut]->Fill(fmssim->pi0(j)->pseudoRapidity());
635  mPhi1[ptbin0][ptbin1][kNCut]->Fill(phi1);
636  mPhi2[ptbin0][ptbin1][kNCut]->Fill(phi2);
637  mDphi[ptbin0][ptbin1][kNCut]->Fill(dphi);
638  if(print) LOG_INFO << Form("RealPi0Pair %d %d dphi=%6.3f",i,j,dphi) << endm;
639  }
640  }
641  }
642  }
643  return kStOK;
644 }
StFmsDiPi0::Finish
Int_t Finish()
Definition: StFmsDiPi0.cxx:220
StFmsDiPi0::Make
Int_t Make()
Definition: StFmsDiPi0.cxx:242
StFmsDbMaker::distanceFromEdge
Float_t distanceFromEdge(Int_t det, Float_t x, Float_t y, int &edge)
get 4 vector assuing m=0 and taking beamline from DB
Definition: StFmsDbMaker.cxx:2142
kStEOF
Definition: Stypes.h:43
StMuDst::event
static StMuEvent * event()
returns pointer to current StMuEvent (class holding the event wise information, e.g. event number, run number)
Definition: StMuDst.h:320
kStOK
Definition: Stypes.h:40
kStErr
Definition: Stypes.h:44
Generated by   doxygen 1.8.5