StRoot  1
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Groups Pages
Efficiency.cxx
1 
2 #include <assert.h>
3 
4 #include <TColor.h>
5 #include <TMultiGraph.h>
6 #include <TGraphErrors.h>
7 #include <TStyle.h>
8 #include <TMinuit.h>
9 #include <TRandom.h>
10 #include <TTree.h>
11 #include <TF1.h>
12 #include <TFile.h>
13 #include <TCanvas.h>
14 #include <TVector3.h>
15 #include <TLorentzVector.h>
16 #include <TMath.h>
17 #include <Riostream.h>
18 
19 #include <StEmcPool/StPhotonCommon/MyEvent.h>
20 #include <StEmcPool/StPhotonCommon/MyPoint.h>
21 #include <StEmcPool/StPhotonCommon/MyMcTrack.h>
22 
23 #include "AnaCuts.h"
24 #include "Pi0Analysis.h"
25 
26 #include "Efficiency.h"
27 using namespace std;
28 ClassImp(Efficiency)
29 
30 Efficiency::Efficiency(const char *input,const char *dir,const char *flag)
31 {
32  gStyle->SetOptDate(0);
33 
34  dirout=TString(dir);
35 
36  mFile=new TFile(input,"OPEN");
37  myEventTree=(TTree*)mFile->Get("mEventTree");
38  ev=new MyEvent();
39  myEventTree->SetBranchAddress("branch",&ev);
40 
41  pseff=new TPostScript((dirout+"_eff.ps").Data(),-111);
42  //init cuts
43  cuts=new AnaCuts(flag);
44  cuts->printCuts();
45 
46  USEWEIGHT=kTRUE;
47  USEPYTHIAWEIGHT=kFALSE;
48  USEBBCSPREAD=kTRUE;
49 
50  mFlag=flag;
51  isMC=kFALSE;
52  isPythia=kFALSE;
53  isDAU=kFALSE;
54  isPP05=kFALSE;
55  if(strcmp(mFlag,"dAu")==0) isDAU=kTRUE;
56  if(strcmp(mFlag,"pp05")==0) isPP05=kTRUE;
57 
58  NEUTRONS=kFALSE;
59  ANTINEUTRONS=kFALSE;
60  PIONS=kFALSE;
61  ETAS=kFALSE;
62  CHARGEDPIONS=kFALSE;
63  PHOTONS=kFALSE;
64  KZEROLONG=kFALSE;
65 
66  nNeutrons=0;
67 
68  htitle=TString("eff. for ");
69 
70  cout<<endl<<"Efficiency constructed!"<<endl<<endl;
71 }
72 
73 Efficiency::~Efficiency()
74 {
75  cout<<endl<<"Efficiency destructed!"<<endl<<endl;
76 }
77 Int_t Efficiency::init()
78 {
79  //acceptance
80  h_genMB=new TH2F("h_genMB","generated vs y and p_{T}",20,0.,20.,30,-0.5,1.5);
81  h_accMB=new TH2F("h_accMB","acceptance vs y and p_{T}",20,0.,20.,30,-0.5,1.5);
82 
83  //conversions
84  h_convGen=new TH1F("h_convGen","generated vs #eta",45,-0.2,1.3);
85  h_convConv=new TH1F("h_convConv","converted vs #eta",45,-0.2,1.3);
86  h_convConvSSD=new TH1F("h_convConvSSD","converted vs rapidity",45,-0.2,1.3);
87  h_convConvSVT=new TH1F("h_convConvSVT","converted vs rapidity",45,-0.2,1.3);
88  h_convConvIFC=new TH1F("h_convConvIFC","converted vs rapidity",45,-0.2,1.3);
89 
90  h_convConvNotCtb=new TH1F("h_convConvNotCtb","converted vs #eta before CTB",45,-0.2,1.3);
91  h_stopRadius=new TH1F("h_stopRadius","conversion radius",1000,0.,450.);
92 
93 
94  //ht
95  h_HT1adc_id=new TH2F("h_HT1adc_id","ht-1 adc vs id",2400,.5,2400.5,800,200.5,1000.5);
96  h_HT2adc_id=new TH2F("h_HT2adc_id","ht-2 adc vs id",2400,.5,2400.5,800,200.5,1000.5);
97 
98  //split clusters:
99  h_splitClusAll=new TH1F("h_splitClusAll","h_splitClusAll",20,0.0,10.0);
100  h_splitClus=new TH1F("h_splitClus","h_splitClus",20,0.0,10.0);
101 
102  //neutrons:
103  h_nbarDet=new TH1F("h_nbarDet","h_nbarDet",cuts->nPtBinsEffMB,cuts->ptBinsEffMB.GetArray());
104  h_nbarIn=new TH1F("h_nbarIn","h_nbarIn",cuts->nPtBinsEffMB,cuts->ptBinsEffMB.GetArray());
105 
106  //QA
107  h_mcdist=new TH1F("h_mcdist","distance to mctrack",1000,0.,5.);
108  h_mcdist2D=new TH2F("h_mcdist2D","distance to mctrack vs pT",20,0.,20.,100,0.,.25);
109  h_dist=new TH1F("h_dist","distance to charged track",1000,0.,200.);
110  h_dist2D=new TH2F("h_dist2D","distance to charged track vs pT",20,0.,20.,200,0.,50.);
111  //eff:
112  h_inputMB=new TH1F("h_inputMB","mc input;p_{T};dN/dp_{T}",cuts->nPtBinsEffMB,cuts->ptBinsEffMB.GetArray());
113  h_inputHT1=new TH1F("h_inputHT1","mc input;p_{T};dN/dp_{T}",cuts->nPtBinsEffHT1,cuts->ptBinsEffHT1.GetArray());
114  h_inputHT2=new TH1F("h_inputHT2","mc input;p_{T};dN/dp_{T}",cuts->nPtBinsEffHT2,cuts->ptBinsEffHT2.GetArray());
115 
116  h_inputDaughtersMB=new TH1F("h_inputDaughtersMB","mc inputDaughters;p_{T};dN/dp_{T}",cuts->nPtBinsEffMB,cuts->ptBinsEffMB.GetArray());
117  h_inputDaughtersHT1=new TH1F("h_inputDaughtersHT1","mc inputDaughters;p_{T};dN/dp_{T}",cuts->nPtBinsEffHT1,cuts->ptBinsEffHT1.GetArray());
118  h_inputDaughtersHT2=new TH1F("h_inputDaughtersHT2","mc inputDaughters;p_{T};dN/dp_{T}",cuts->nPtBinsEffHT2,cuts->ptBinsEffHT2.GetArray());
119 
120  h_recoMB=new TH1F("h_recoMB","reco.;p_{T};dN/dp_{T}",cuts->nPtBinsEffMB,cuts->ptBinsEffMB.GetArray());
121  h_recoHT1=new TH1F("h_recoHT1","reco.;p_{T};dN/dp_{T}",cuts->nPtBinsEffHT1,cuts->ptBinsEffHT1.GetArray());
122  h_recoHT2=new TH1F("h_recoHT2","reco.;p_{T};dN/dp_{T}",cuts->nPtBinsEffHT2,cuts->ptBinsEffHT2.GetArray());
123 
124  h_recoDaughtersMB=new TH1F("h_recoDaughtersMB","recoDaughters.;p_{T};dN/dp_{T}",cuts->nPtBinsEffMB,cuts->ptBinsEffMB.GetArray());
125  h_recoDaughtersHT1=new TH1F("h_recoDaughtersHT1","recoDaughters.;p_{T};dN/dp_{T}",cuts->nPtBinsEffHT1,cuts->ptBinsEffHT1.GetArray());
126  h_recoDaughtersHT2=new TH1F("h_recoDaughtersHT2","recoDaughters.;p_{T};dN/dp_{T}",cuts->nPtBinsEffHT2,cuts->ptBinsEffHT2.GetArray());
127 
128  h_minvMB=new TH2F("h_minvMB","reco.;p_{T};m_{inv}",
129  cuts->nPtBinsEffMB,cuts->ptBinsEffMB.GetArray(),cuts->nMinvBinsMB,cuts->mInvLowMB,cuts->mInvHighMB);
130  h_minvHT1=new TH2F("h_minvHT1","reco.;p_{T};m_{inv}",
131  cuts->nPtBinsEffHT1,cuts->ptBinsEffHT1.GetArray(),cuts->nMinvBinsHT1,cuts->mInvLowHT1,cuts->mInvHighHT1);
132  h_minvHT2=new TH2F("h_minvHT2","reco.;p_{T};m_{inv}",
133  cuts->nPtBinsEffHT2,cuts->ptBinsEffHT2.GetArray(),cuts->nMinvBinsHT2,cuts->mInvLowHT2,cuts->mInvHighHT2);
134 
135  h_matrixMB=new TH2F("h_matrixMB","pT gen. vs reco. MB",80,0.,20.,80,0.,20.);
136 
137  h_etaphi=new TH2F("h_etaphi","eta/phi neutral points",300,0.,1.,1800,-TMath::Pi(),TMath::Pi());
138  h_pythiaPions=new TH1F("h_pythiaPions","pythia pions vs pT",200,0.,50.);
139  h_pythiaPhotons=new TH1F("h_pythiaPhotons","pythia photons vs pT",200,0.,50.);
140  h_pythiaPartonPt=new TH1F("h_pythiaPartonPt","pT of pythia process",200,0.,50);
141 
142  h_clusterWidth=new TH2F("h_clusterWidth","width BSMD eta+phi",20,0.,20.,100,0.,0.05);
143  h_energyRatio=new TH2F("h_energyRatio","energy ratio BSMD/BTOW",20,0.,20.,160,0.,8.);
144  h_towclusRatio=new TH2F("h_towclusRatio","hightow/cluster",20,0.,20.,50,0.,1.);
145 
146  h_vzMB=new TH1F("h_vzMB","z-vertex MB",480,-120.,120.);
147  h_vzMB->Sumw2();
148 
149  h_asymmMB=new TH2F("h_asymmMB","asymmetry of inv mass pairs MB",cuts->nPtBinsMB,cuts->ptBinsMB.GetArray(),40,0.0,1.0);
150  h_asymmMB->Sumw2();
151  h_asymmHT1=new TH2F("h_asymmHT1","asymmetry of inv mass pairs HT1",cuts->nPtBinsHT1,cuts->ptBinsHT1.GetArray(),40,0.0,1.0);
152  h_asymmHT1->Sumw2();
153  h_asymmHT2=new TH2F("h_asymmHT2","asymmetry of inv mass pairs HT2",cuts->nPtBinsHT2,cuts->ptBinsHT2.GetArray(),40,0.0,1.0);
154  h_asymmHT2->Sumw2();
155 
156  h_smdeta1=new TH2F("h_smdeta1","n eta strips in cluster vs p_{T}",20,0.,10.,10,-0.5,9.5);
157  h_smdphi1=new TH2F("h_smdphi1","n phi strips in cluster vs p_{T}",20,0.,10.,10,-0.5,9.5);
158  h_smdeta2=new TH2F("h_smdeta2","eta energy ratio in cluster vs p_{T}",20,0.,10.,100,0.0,10.);
159  h_smdphi2=new TH2F("h_smdphi2","phi energy ratio in cluster vs p_{T}",20,0.,10.,100,0.0,10.);
160 
161  h_energyeta=new TH1F("h_energyeta","energy BSMDE",40,0.,10.);
162  h_energyphi=new TH1F("h_energyeta","energy BSMDP",40,0.,10.);
163 
164  h_pionsVsEtaMB=new TH1F("h_pionsVsEta","dN_{#pi^{0}}/dy data",40,-0.5,1.5);
165  h_pionsVsEtaMB->Sumw2();
166 
167 
168  h_smdeta1->Sumw2();
169  h_smdphi1->Sumw2();
170  h_smdeta2->Sumw2();
171  h_smdphi2->Sumw2();
172 
173 
174  h_pythiaPartonPt->Sumw2();
175  h_inputMB->Sumw2();
176  h_inputHT1->Sumw2();
177  h_inputHT2->Sumw2();
178  h_inputDaughtersMB->Sumw2();
179  h_inputDaughtersHT1->Sumw2();
180  h_inputDaughtersHT2->Sumw2();
181  h_recoMB->Sumw2();
182  h_recoHT1->Sumw2();
183  h_recoHT2->Sumw2();
184  h_recoDaughtersMB->Sumw2();
185  h_recoDaughtersHT1->Sumw2();
186  h_recoDaughtersHT2->Sumw2();
187  h_minvMB->Sumw2();
188  h_minvHT1->Sumw2();
189  h_minvHT2->Sumw2();
190  h_matrixMB->Sumw2();
191  h_pythiaPions->Sumw2();
192  h_clusterWidth->Sumw2();
193  h_energyRatio->Sumw2();
194  h_towclusRatio->Sumw2();
195 
196 
197 
198  return 1;
199 }
200 Int_t Efficiency::make(Int_t evmax)
201 {
202  Float_t WEIGHT=0.;
203  Int_t pid=0;
204  Int_t i=0;
205  while(myEventTree->GetEntry(i))
206  {
207  if(evmax&&i>evmax)
208  {
209  cout<<"reached evmax,"<<endl;
210  cout<<"abort loop!"<<endl;
211  break;
212  }
213  if(i%10000==0) cout<<"processing event: "<<i<<endl;
214 
215 
216  MyMcTrack *mcTr=0;
217  Float_t ptInput=0;
218  Float_t rapInput=0;
219  if(!isPythia){
220  mcTr=ev->getMcTrack();
221  ptInput=mcTr->momentum().Pt();
222  rapInput=mcTr->momentum().PseudoRapidity();
223  }
224 
225  if(i==0 && !isPythia){
226  cout<<"not pythia, plain MC"<<endl;
227  pid=mcTr->id();
228  if(pid==111){PIONS=kTRUE;htitle=TString("pions");}
229  else if(pid==221) {ETAS=kTRUE;htitle=TString("etas");}
230  else if(pid==2112) {NEUTRONS=kTRUE;htitle=TString("neutrons");}
231  else if(pid==-2112) {ANTINEUTRONS=kTRUE;htitle=TString("antineutrons");}
232  else if(pid==211) {CHARGEDPIONS=kTRUE;htitle=TString("piplus");}
233  else if(pid==22) {PHOTONS=kTRUE;htitle=TString("single gamma");}
234  else if(pid==130) {KZEROLONG=kTRUE;htitle=TString("single kzeroL");}
235  else assert(0);
236 
237  htitle+=TString(" in ");
238  htitle+=TString(mFlag);
239  }
240  else if(isPythia){
241  htitle=TString("pions pythia");
242  }
243 
244  if(NEUTRONS||ANTINEUTRONS){
245  TLorentzVector lv(mcTr->momentum().X(),mcTr->momentum().Y(),mcTr->momentum().Z(),sqrt(mcTr->momentum().Mag2()+1.));
246  rapInput=lv.Rapidity();
247  }
248 
249  //get vertex
250  TVector3 vPos=ev->vertex();
251  //get vertex weight dAu + pp separately:
252  WEIGHT=getWeightVertex(vPos.Z());
253  h_vzMB->Fill(vPos.Z(),WEIGHT);
254 
255 
256  //add spread from bbc:
257  if(USEBBCSPREAD && isPP05 && gRandom->Rndm()<0.36){
258  double bbc_spread=gRandom->Gaus(0.,43.5);
259  vPos.SetZ(vPos.Z()+bbc_spread);
260  }
261 
262  //set trigger from simu for dAu
263  if(isDAU){
264  Int_t tr=1;
265  if(ev->highTowerAdc()>cuts->ht1AdcCUT) tr+=2;
266  if(ev->highTowerAdc()>cuts->ht2AdcCUT) tr+=4;
267  ev->setTrigger(tr);
268  }
269  else if(!isPP05){
270  cout<<"wrong flag"<<endl;
271  assert(0);
272  }
273 
274  if(ev->trigger()&2) h_HT1adc_id->Fill(ev->highTowerId(),ev->highTowerAdc());
275  if(ev->trigger()&4) h_HT2adc_id->Fill(ev->highTowerId(),ev->highTowerAdc());
276 
277 
278  //event cuts
279  ev->setMomentumInTpc(999999.);
280 
281  if(isDAU){
282  if(!cuts->isEventOK(ev,"dAu")) {i++;continue;}
283  }
284  else if(isPP05){
285  if(!cuts->isEventOK(ev,"pp05")) {i++;continue;}
286  }
287  else assert(0);
288 
289 
290  //event ok, count neutrons:
291  nNeutrons=nNeutrons+1.;
292  float pt_rnd=10.*gRandom->Rndm();
293  h_nbarIn->Fill(pt_rnd,exp(-0.5*pt_rnd));
294  h_nbarDet->Fill(mcTr->momentum().Pt());
295 
296  if(ETAS||PIONS){
297  if(ev->numberOfMcPhotons()!=2 && !isPythia){
298  i++;
299  continue;
300  }
301  }
302 
303  //fill input:
304  if(isPythia && 0){
305  Float_t pypT=ev->partonPt();
306  if(pypT<15.){
307  WEIGHT=1.;
308  }
309  else if(pypT<25.){
310  WEIGHT=(216000/208000)*(0.0003895/0.002228);
311  }
312  else if(pypT>25.){
313  WEIGHT=(216000/208000)*(1.016e-005/0.002228);
314  }
315  }
316  if(isPythia){
317  h_pythiaPartonPt->Fill(ev->partonPt(),WEIGHT);
318  for(Int_t it=0;it<ev->getMcPionArray()->GetEntries();it++){
319  MyMcTrack *pyt=(MyMcTrack*)ev->getMcPionArray()->At(it);
320  if(pyt->momentum().PseudoRapidity()<cuts->rapPionMinCUT||
321  pyt->momentum().PseudoRapidity()>cuts->rapPionMaxCUT) continue;
322  h_pythiaPions->Fill(pyt->momentum().Pt(),WEIGHT);
323  if(ev->trigger()&1){
324  h_inputMB->Fill(pyt->momentum().Pt(),WEIGHT);
325  }
326  if(ev->trigger()&1){
327  h_inputHT1->Fill(pyt->momentum().Pt(),WEIGHT);
328  }
329  if(ev->trigger()&1){
330  h_inputHT2->Fill(pyt->momentum().Pt(),WEIGHT);
331  }
332  }
333  for(Int_t it=0;it<ev->getMcPhotonArray()->GetEntries();it++){
334  MyMcTrack *pyt=(MyMcTrack*)ev->getMcPhotonArray()->At(it);
335  if(pyt->momentum().PseudoRapidity()<cuts->rapPionMinCUT||
336  pyt->momentum().PseudoRapidity()>cuts->rapPionMaxCUT) continue;
337  h_pythiaPhotons->Fill(pyt->momentum().Pt(),WEIGHT);
338  if(ev->trigger()&1){
339  h_inputDaughtersMB->Fill(pyt->momentum().Pt(),WEIGHT);
340  }
341  if(ev->trigger()&1){
342  h_inputDaughtersHT1->Fill(pyt->momentum().Pt(),WEIGHT);
343  }
344  if(ev->trigger()&1){
345  h_inputDaughtersHT2->Fill(pyt->momentum().Pt(),WEIGHT);
346  }
347  }
348  }
349 
350 
351 
352  if(NEUTRONS){
353  WEIGHT*=getWeightNeutrons(ptInput);
354  if(rapInput>cuts->rapidityMinCUT&&rapInput<cuts->rapidityMaxCUT){
355  if(ev->trigger()&1) h_inputMB->Fill(mcTr->momentum().Pt(),WEIGHT);
356  if(ev->trigger()&1) h_inputHT1->Fill(mcTr->momentum().Pt(),WEIGHT);
357  if(ev->trigger()&1) h_inputHT2->Fill(mcTr->momentum().Pt(),WEIGHT);
358  }
359  }
360 
361  if(ANTINEUTRONS || KZEROLONG){
362  WEIGHT*=getWeightAntiNeutrons(ptInput);
363  if(rapInput>cuts->rapidityMinCUT&&rapInput<cuts->rapidityMaxCUT){
364  if(ev->trigger()&1) h_inputMB->Fill(mcTr->momentum().Pt(),WEIGHT);
365  if(ev->trigger()&1) h_inputHT1->Fill(mcTr->momentum().Pt(),WEIGHT);
366  if(ev->trigger()&1) h_inputHT2->Fill(mcTr->momentum().Pt(),WEIGHT);
367  }
368  }
369 
370  if(PIONS||ETAS){
371  if(PIONS) WEIGHT*=getWeightPions(ptInput);
372  else if(ETAS) WEIGHT*=getWeightEtas(ptInput);
373 
374  //digamma
375  if(rapInput>cuts->rapPionMinCUT&&rapInput<cuts->rapPionMaxCUT){
376  if(ev->trigger()&1) h_inputMB->Fill(mcTr->momentum().Pt(),WEIGHT);
377  if(ev->trigger()&1) h_inputHT1->Fill(mcTr->momentum().Pt(),WEIGHT);
378  if(ev->trigger()&1) h_inputHT2->Fill(mcTr->momentum().Pt(),WEIGHT);
379  }
380  MyMcTrack *daughterA=(MyMcTrack*)ev->getMcPhotonArray()->At(0);
381  MyMcTrack *daughterB=(MyMcTrack*)ev->getMcPhotonArray()->At(1);
382 
383  //if(TMath::Abs(daughterA->energy()-daughterB->energy())/(daughterA->energy()+daughterB->energy())<0.7){
384  // i++;
385  // continue;
386  //}
387 
388  //for acceptance
389  h_genMB->Fill(mcTr->momentum().Pt(),mcTr->momentum().PseudoRapidity());
390  if(daughterA->position().PseudoRapidity() > cuts->etaMinCUT &&
391  daughterA->position().PseudoRapidity() < cuts->etaMaxCUT){
392  if(daughterB->position().PseudoRapidity() > cuts->etaMinCUT &&
393  daughterB->position().PseudoRapidity() < cuts->etaMaxCUT){
394  h_accMB->Fill(mcTr->momentum().Pt(),mcTr->momentum().PseudoRapidity());
395  }
396  }
397 
398  if(daughterA->momentum().PseudoRapidity() > cuts->rapidityMinCUT &&
399  daughterA->momentum().PseudoRapidity() < cuts->rapidityMaxCUT){
400  if(ev->trigger()&1) h_inputDaughtersMB->Fill(daughterA->momentum().Pt(),WEIGHT);
401  if(ev->trigger()&1) h_inputDaughtersHT1->Fill(daughterA->momentum().Pt(),WEIGHT);
402  if(ev->trigger()&1) h_inputDaughtersHT2->Fill(daughterA->momentum().Pt(),WEIGHT);
403  }
404  if(daughterB->momentum().PseudoRapidity() > cuts->rapidityMinCUT &&
405  daughterB->momentum().PseudoRapidity() < cuts->rapidityMaxCUT){
406  if(ev->trigger()&1) h_inputDaughtersMB->Fill(daughterB->momentum().Pt(),WEIGHT);
407  if(ev->trigger()&1) h_inputDaughtersHT1->Fill(daughterB->momentum().Pt(),WEIGHT);
408  if(ev->trigger()&1) h_inputDaughtersHT2->Fill(daughterB->momentum().Pt(),WEIGHT);
409  }
410  }
411 
412  if(CHARGEDPIONS){
413  WEIGHT*=getWeightPions(ptInput);
414  if(rapInput>cuts->rapidityMinCUT&&rapInput<cuts->rapidityMaxCUT){
415  if(ev->trigger()&1) h_inputMB->Fill(mcTr->momentum().Pt(),WEIGHT);
416  if(ev->trigger()&1) h_inputHT1->Fill(mcTr->momentum().Pt(),WEIGHT);
417  if(ev->trigger()&1) h_inputHT2->Fill(mcTr->momentum().Pt(),WEIGHT);
418  }
419  }
420 
421  if(PHOTONS){
422  WEIGHT*=getWeightPhotons(ptInput);
423  if(rapInput>cuts->rapidityMinCUT&&rapInput<cuts->rapidityMaxCUT){
424  if(ev->trigger()&1) h_inputMB->Fill(mcTr->momentum().Pt(),WEIGHT);
425  if(ev->trigger()&1) h_inputHT1->Fill(mcTr->momentum().Pt(),WEIGHT);
426  if(ev->trigger()&1) h_inputHT2->Fill(mcTr->momentum().Pt(),WEIGHT);
427 
428  h_stopRadius->Fill(mcTr->stopRadius());
429  h_convGen->Fill(mcTr->momentum().PseudoRapidity(),getWeightVertex(vPos.Z()));
430  if(mcTr->stopRadius()>0&&mcTr->stopRadius()<223.5){
431  h_convConv->Fill(mcTr->momentum().PseudoRapidity(),getWeightVertex(vPos.Z()));
432  }
433  if(mcTr->stopRadius()>0&&mcTr->stopRadius()<200.){
434  h_convConvNotCtb->Fill(mcTr->momentum().PseudoRapidity(),getWeightVertex(vPos.Z()));
435  }
436  //SVT
437  if(mcTr->stopRadius()>1.&&mcTr->stopRadius()<20.){
438  h_convConvSVT->Fill(mcTr->momentum().PseudoRapidity(),getWeightVertex(vPos.Z()));
439  }
440  //SSD
441  if(mcTr->stopRadius()>20.&&mcTr->stopRadius()<40.){
442  h_convConvSSD->Fill(mcTr->momentum().PseudoRapidity(),getWeightVertex(vPos.Z()));
443  }
444  //IFC
445  if(mcTr->stopRadius()>40.&&mcTr->stopRadius()<55.){
446  h_convConvIFC->Fill(mcTr->momentum().PseudoRapidity(),getWeightVertex(vPos.Z()));
447  }
448  }
449  }
450 
451  //check split clusters:
452  int nclusters=0;
453 
454  MyPoint *p;
455  MyPoint *pp;
456  TClonesArray *clA=(TClonesArray*)ev->getPointArray();
457  for(Int_t j=0;j<clA->GetEntries();j++)
458  {
459  p=(MyPoint*)clA->At(j);
460  TVector3 pPos=p->position();
461  TVector3 pMom=pPos-vPos;
462  pMom.SetMag(p->energy());
463 
464  h_dist2D->Fill(pMom.Pt(),TMath::Abs(p->distanceToTrack()));
465 
466  if(!cuts->isPointOK(p,vPos)) continue;
467 
468  //chaeck split clusters:
469  nclusters++;
470  if(ev->trigger()&1 && pMom.Pt()>1.){
471  h_etaphi->Fill(pPos.PseudoRapidity(),pPos.Phi(),WEIGHT);
472  }
473  Float_t PdistMC=9999.;
474 
475  h_dist->Fill(TMath::Abs(p->distanceToTrack()));
476 
477  MyMcTrack *closestTrack=0;
478  if(PIONS||ETAS){
479  PdistMC=9999.;//reset
480  MyMcTrack *trA=(MyMcTrack*)ev->getMcPhotonArray()->At(0);
481  MyMcTrack *trB=(MyMcTrack*)ev->getMcPhotonArray()->At(1);
482  TVector3 mcPosA=trA->position();
483  TVector3 mcPosB=trB->position();
484  Float_t An=mcPosA.PseudoRapidity();
485  Float_t Ap=mcPosA.Phi();
486  Float_t Bn=mcPosB.PseudoRapidity();
487  Float_t Bp=mcPosB.Phi();
488  //point p:
489  Float_t Pn=pPos.PseudoRapidity();
490  Float_t Pp=pPos.Phi();
491  Float_t dphiAP=TMath::Abs(Ap-Pp);
492  while(dphiAP>2.*TMath::Pi()) dphiAP-=2.*TMath::Pi();
493  Float_t dphiBP=TMath::Abs(Bp-Pp);
494  while(dphiBP>2.*TMath::Pi()) dphiBP-=2.*TMath::Pi();
495  Float_t dpA=sqrt(pow(An-Pn,2)+pow(dphiAP,2));
496  Float_t dpB=sqrt(pow(Bn-Pn,2)+pow(dphiBP,2));
497  //distance to mc track:
498  PdistMC=dpA<dpB ? dpA : dpB;
499  closestTrack=dpA<dpB ? trA : trB;
500 
501  h_mcdist->Fill(PdistMC);
502  h_mcdist2D->Fill(pMom.Pt(),PdistMC);
503  }
504 
505  //gamma:
506  if(p->distanceToTrack()>cuts->photonCUT){
507 
508  if(pMom.PseudoRapidity()>cuts->rapidityMinCUT && pMom.PseudoRapidity()<cuts->rapidityMaxCUT){
509  if(ETAS||PIONS||isPythia){
510 
511  if(ev->trigger()&1) h_recoDaughtersMB->Fill(pMom.Pt(),WEIGHT);
512  if(ev->trigger()&2) h_recoDaughtersHT1->Fill(pMom.Pt(),WEIGHT);
513  if(ev->trigger()&4) h_recoDaughtersHT2->Fill(pMom.Pt(),WEIGHT);
514 
515  }
516  else{
517  if(ev->trigger()&1) h_recoMB->Fill(pMom.Pt(),WEIGHT);
518  if(ev->trigger()&2) h_recoHT1->Fill(pMom.Pt(),WEIGHT);
519  if(ev->trigger()&4) h_recoHT2->Fill(pMom.Pt(),WEIGHT);
520 
521  if(ev->trigger()&1) h_matrixMB->Fill(pMom.Pt(),ptInput);
522  }
523 
524  if(ev->trigger()&1){
525  h_clusterWidth->Fill(pMom.Pt(),sqrt(p->widthEta()*p->widthEta()+p->widthPhi()*p->widthPhi()),WEIGHT);
526  h_energyRatio->Fill(pMom.Pt(),(p->energyEta()+p->energyPhi())/p->energy(),WEIGHT);
527  h_towclusRatio->Fill(pMom.Pt(),p->towerClusterEnergy(0)/p->energy(),WEIGHT);
528 
529  if(p->distanceToTrack()>30.){
530  h_smdeta1->Fill(pMom.Pt(),p->nHitsEta(),WEIGHT);
531  h_smdphi1->Fill(pMom.Pt(),p->nHitsPhi(),WEIGHT);
532  h_smdeta2->Fill(pMom.Pt(),p->energyEta()/p->energy(),WEIGHT);
533  h_smdphi2->Fill(pMom.Pt(),p->energyPhi()/p->energy(),WEIGHT);
534 
535  h_energyeta->Fill(p->energyEta());
536  h_energyphi->Fill(p->energyPhi());
537  }
538  }
539 
540 
541  }
542 
543  }
544 
545  for(Int_t jj=0;(ETAS||PIONS||PHOTONS||isPythia)&&jj<clA->GetEntries();jj++)
546  {
547  if(jj<=j) continue;
548  pp=(MyPoint*)clA->At(jj);
549  TVector3 ppPos=pp->position();
550  TVector3 ppMom=ppPos-vPos;
551  ppMom.SetMag(pp->energy());
552  if(!cuts->isPointOK(pp,vPos)) continue;
553 
554  /*Float_t PPdistMC=9999.;
555  if(!isPythia){
556  MyMcTrack *trA=(MyMcTrack*)ev->getMcPhotonArray()->At(0);
557  MyMcTrack *trB=(MyMcTrack*)ev->getMcPhotonArray()->At(1);
558  TVector3 mcPosA=trA->position();
559  TVector3 mcPosB=trB->position();
560  Float_t An=mcPosA.PseudoRapidity();
561  Float_t Ap=mcPosA.Phi();
562  Float_t Bn=mcPosB.PseudoRapidity();
563  Float_t Bp=mcPosB.Phi();
564 
565  //point pp:
566  Float_t PPn=ppPos.PseudoRapidity();
567  Float_t PPp=ppPos.Phi();
568  Float_t dphiAPP=TMath::Abs(Ap-PPp);
569  while(dphiAPP>2.*TMath::Pi()) dphiAPP-=2.*TMath::Pi();
570  Float_t dphiBPP=TMath::Abs(Bp-PPp);
571  while(dphiBPP>2.*TMath::Pi()) dphiBPP-=2.*TMath::Pi();
572  Float_t dppA=sqrt(pow(An-PPn,2)+pow(dphiAPP,2));
573  Float_t dppB=sqrt(pow(Bn-PPn,2)+pow(dphiBPP,2));
574  PPdistMC=dppA<dppB ? dppA : dppB;
575  }*/
576 
577  //two neutrals:
578  TVector3 pi0Mom=pMom+ppMom;
579  Float_t angle=pMom.Angle(ppMom);
580  Float_t minv=sqrt(2.*p->energy()*pp->energy()*(1. - cos(angle)));
581  Float_t pTpion=pi0Mom.Pt();
582  Float_t etapion=pi0Mom.PseudoRapidity();
583  Float_t asymm=TMath::Abs(p->energy()-pp->energy())/(p->energy()+pp->energy());
584 
585  if(etapion<cuts->rapPionMinCUT||etapion>cuts->rapPionMaxCUT) continue;
586 
587  if(ev->trigger()&1){
588  if(asymm<=cuts->asymmetryCUT) h_minvMB->Fill(pTpion,minv,WEIGHT);
589  if(pTpion>1.){
590  if(minv>0.08 && minv<0.20){
591  h_asymmMB->Fill(pTpion,asymm,WEIGHT);
592  if(asymm<=cuts->asymmetryCUT){
593  h_pionsVsEtaMB->Fill(pi0Mom.PseudoRapidity(),WEIGHT);
594  h_matrixMB->Fill(pTpion,ptInput,WEIGHT);
595  }
596  }
597  }
598  }
599  if(ev->trigger()&2){
600  if(asymm<=cuts->asymmetryCUT) h_minvHT1->Fill(pTpion,minv,WEIGHT);
601  if(pTpion>4.){
602  if(minv>0.10 && minv<0.20){
603  h_asymmHT1->Fill(pTpion,asymm,WEIGHT);
604  }
605  }
606  }
607  if(ev->trigger()&4){
608  if(asymm<=cuts->asymmetryCUT) h_minvHT2->Fill(pTpion,minv,WEIGHT);
609  if(pTpion>6.){
610  if(minv>0.10 && minv<0.2){
611  h_asymmHT2->Fill(pTpion,asymm,WEIGHT);
612  }
613  }
614  }
615 
616 
617  }
618 
619 
620  }
621 
622  if(PHOTONS && nclusters>0) h_splitClusAll->Fill(mcTr->momentum().Pt());
623  if(PHOTONS && nclusters>1) h_splitClus->Fill(mcTr->momentum().Pt());
624 
625  i++;
626  }
627 
628 
629 
630  if(ETAS||PIONS||isPythia){//extract yield with regular routine:
631  Pi0Analysis *pi0=new Pi0Analysis((const char*)(dirout+"_effbins.ps").Data(),"/dev/null",mFlag);
632  pi0->init("/dev/null");
633  if(!isPythia) pi0->setMC(kTRUE);
634  h_recoMB=new TH1F(*pi0->getYield(h_minvMB,"mb"));
635  h_recoHT1=new TH1F(*pi0->getYield(h_minvHT1,"ht1"));
636  h_recoHT2=new TH1F(*pi0->getYield(h_minvHT2,"ht2"));
637  pi0->storeCanvases((dirout+"_canvases.root").Data());
638  delete pi0;
639  // !!!! calculate from dN/dpt to #/bin:
640  for(Int_t i=1;i<=h_recoMB->GetNbinsX();i++)
641  {
642  h_recoMB->SetBinContent(i,h_recoMB->GetBinContent(i)*h_recoMB->GetBinWidth(i));
643  h_recoMB->SetBinError(i,h_recoMB->GetBinError(i)*h_recoMB->GetBinWidth(i));
644  }
645  for(Int_t i=1;i<=h_recoHT1->GetNbinsX();i++)
646  {
647  h_recoHT1->SetBinContent(i,h_recoHT1->GetBinContent(i)*h_recoHT1->GetBinWidth(i));
648  h_recoHT1->SetBinError(i,h_recoHT1->GetBinError(i)*h_recoHT1->GetBinWidth(i));
649  }
650  for(Int_t i=1;i<=h_recoHT2->GetNbinsX();i++)
651  {
652  h_recoHT2->SetBinContent(i,h_recoHT2->GetBinContent(i)*h_recoHT2->GetBinWidth(i));
653  h_recoHT2->SetBinError(i,h_recoHT2->GetBinError(i)*h_recoHT2->GetBinWidth(i));
654  }
655  }
656 
657  gStyle->SetPalette(1);
658  gStyle->SetStatStyle(0);
659 
660  pseff->On();
661  TCanvas *c=new TCanvas("c","c",600,800);
662  c->Divide(3,4);
663  c->cd(1);
664  gPad->SetLogy();
665  h_inputMB->SetLineColor(1);
666  h_inputMB->SetLineWidth(2);
667  h_inputMB->SetTitle((htitle+TString(" MB")).Data());
668  h_inputMB->Draw("hist");
669  h_recoMB->SetLineColor(1);
670  h_recoMB->SetLineWidth(2);
671  h_recoMB->SetLineStyle(2);
672  if(h_recoMB->GetMinimum()>0.) h_inputMB->SetMinimum(h_recoMB->GetMinimum()/10. );
673  h_recoMB->Draw("histsame");
674  c->cd(2);
675  gPad->SetLogy();
676  h_inputHT1->SetLineColor(TColor::GetColor(24,101,24));
677  h_inputHT1->SetLineWidth(2);
678  h_inputHT1->SetTitle((htitle+TString(" HT1")).Data());
679  h_inputHT1->Draw("hist");
680  h_recoHT1->SetLineColor(TColor::GetColor(24,101,24));
681  h_recoHT1->SetLineWidth(2);
682  h_recoHT1->SetLineStyle(2);
683  if(h_recoHT1->GetMinimum()>0.) h_inputHT1->SetMinimum(h_recoHT1->GetMinimum()/10. );
684  h_recoHT1->Draw("histsame");
685  c->cd(3);
686  gPad->SetLogy();
687  h_inputHT2->SetLineColor(TColor::GetColor(24,28,174));
688  h_inputHT2->SetLineWidth(2);
689  h_inputHT2->SetTitle((htitle+TString(" HT2")).Data());
690  h_inputHT2->Draw("hist");
691  h_recoHT2->SetLineColor(TColor::GetColor(24,28,174));
692  h_recoHT2->SetLineWidth(2);
693  h_recoHT2->SetLineStyle(2);
694  if(h_recoHT2->GetMinimum()>0.) h_inputHT2->SetMinimum(h_recoHT2->GetMinimum()/10. );
695  h_recoHT2->Draw("histsame");
696 
697  c->cd(4);
698  if(NEUTRONS||ANTINEUTRONS||KZEROLONG){
699  gPad->SetLogy();
700  }
701  h_effMB=new TH1F(*h_recoMB);
702  h_effMB->SetNameTitle("h_effMB","efficiency MB;p_{T};#epsilon");
703  h_effMB->SetLineStyle(1);
704  h_effMB->Divide(h_inputMB);
705  h_effMB->Draw("pe");
706  c->cd(5);
707  if(NEUTRONS||ANTINEUTRONS||KZEROLONG){
708  gPad->SetLogy();
709  }
710  h_effHT1=new TH1F(*h_recoHT1);
711  h_effHT1->SetNameTitle("h_effHT1","efficiency HT1;p_{T};#epsilon");
712  h_effHT1->SetLineStyle(1);
713  h_effHT1->Divide(h_inputHT1);
714  h_effHT1->Draw("pe");
715  c->cd(6);
716  if(NEUTRONS||ANTINEUTRONS||KZEROLONG){
717  gPad->SetLogy();
718  }
719  h_effHT2=new TH1F(*h_recoHT2);
720  h_effHT2->SetNameTitle("h_effHT2","efficiency HT2;p_{T};#epsilon");
721  h_effHT2->SetLineStyle(1);
722  h_effHT2->Divide(h_inputHT2);
723  h_effHT2->Draw("pe");
724 
725  c->cd(7);
726  gPad->SetLogy();
727  h_inputDaughtersMB->SetLineColor(1);
728  h_inputDaughtersMB->SetLineWidth(2);
729  h_inputDaughtersMB->Draw("hist");
730  h_recoDaughtersMB->SetLineColor(1);
731  h_recoDaughtersMB->SetLineWidth(2);
732  h_recoDaughtersMB->SetLineStyle(2);
733  h_recoDaughtersMB->Draw("histsame");
734 
735  c->cd(8);
736  gPad->SetLogy();
737  h_inputDaughtersHT1->SetLineColor(TColor::GetColor(24,101,24));
738  h_inputDaughtersHT1->SetLineWidth(2);
739  h_inputDaughtersHT1->Draw("hist");
740  h_recoDaughtersHT1->SetLineColor(TColor::GetColor(24,101,24));
741  h_recoDaughtersHT1->SetLineWidth(2);
742  h_recoDaughtersHT1->SetLineStyle(2);
743  h_recoDaughtersHT1->Draw("histsame");
744  c->cd(9);
745  gPad->SetLogy();
746  h_inputDaughtersHT2->SetLineColor(TColor::GetColor(24,28,174));
747  h_inputDaughtersHT2->SetLineWidth(2);
748  h_inputDaughtersHT2->Draw("hist");
749  h_recoDaughtersHT2->SetLineColor(TColor::GetColor(24,28,174));
750  h_recoDaughtersHT2->SetLineWidth(2);
751  h_recoDaughtersHT2->SetLineStyle(2);
752  h_recoDaughtersHT2->Draw("histsame");
753 
754  c->cd(10);
755  h_effDaughtersMB=new TH1F(*h_recoDaughtersMB);
756  h_effDaughtersMB->SetNameTitle("h_effDaughtersMB","efficiency daughters MB;p_{T};#epsilon");
757  h_effDaughtersMB->SetLineStyle(1);
758  h_effDaughtersMB->Divide(h_inputDaughtersMB);
759  h_effDaughtersMB->Draw("pe");
760  c->cd(11);
761  h_effDaughtersHT1=new TH1F(*h_recoDaughtersHT1);
762  h_effDaughtersHT1->SetNameTitle("h_effDaughtersHT1","efficiency daughters HT1;p_{T};#epsilon");
763  h_effDaughtersHT1->SetLineStyle(1);
764  h_effDaughtersHT1->Divide(h_inputDaughtersHT1);
765  h_effDaughtersHT1->Draw("pe");
766  c->cd(12);
767  h_effDaughtersHT2=new TH1F(*h_recoDaughtersHT2);
768  h_effDaughtersHT2->SetNameTitle("h_effDaughtersHT2","efficiency daughters HT2;p_{T};#epsilon");
769  h_effDaughtersHT2->SetLineStyle(1);
770  h_effDaughtersHT2->Divide(h_inputDaughtersHT2);
771  h_effDaughtersHT2->Draw("pe");
772 
773 
774  c->cd(0);
775  c->Update();
776  c->SaveAs((dirout+"_effplots.root").Data());
777  pseff->Close();
778 
779 
780  return 1;
781 }
782 Int_t Efficiency::finish()
783 {
784  //save histograms
785  mFileOut=new TFile((dirout+"_eff.root").Data(),"RECREATE");
786  mFileOut->cd();
787 
788  h_genMB->Write();
789  h_accMB->Write();
790 
791  h_HT1adc_id->Write();
792  h_HT2adc_id->Write();
793 
794  h_splitClus->Write();
795  h_splitClusAll->Write();
796 
797  h_convGen->Write();
798  h_convConv->Write();
799  h_convConvSSD->Write();
800  h_convConvSVT->Write();
801  h_convConvIFC->Write();
802 
803  h_convConvNotCtb->Write();
804  h_stopRadius->Write();
805 
806  h_nbarDet->Write();
807  h_nbarIn->Write();
808 
809  h_clusterWidth->Write();
810  h_energyRatio->Write();
811  h_towclusRatio->Write();
812 
813  h_smdeta1->Write();
814  h_smdphi1->Write();
815  h_smdeta2->Write();
816  h_smdphi2->Write();
817 
818 
819  h_pythiaPions->Write();
820  h_pythiaPhotons->Write();
821  h_pythiaPartonPt->Write();
822  //
823  h_mcdist->Write();
824  h_mcdist2D->Write();
825  h_dist->Write();
826  h_dist2D->Write();
827 
828  h_effMB->Write();
829  h_effHT1->Write();
830  h_effHT2->Write();
831 
832  h_effDaughtersMB->Write();
833  h_effDaughtersHT1->Write();
834  h_effDaughtersHT2->Write();
835 
836  h_matrixMB->Write();
837 
838  h_etaphi->Write();
839 
840  h_minvMB->Write();
841  h_minvHT1->Write();
842  h_minvHT2->Write();
843 
844  h_asymmMB->Write();
845  h_asymmHT1->Write();
846  h_asymmHT2->Write();
847 
848  h_pionsVsEtaMB->Write();
849 
850  h_vzMB->Write();
851 
852  h_energyeta->Write();
853  h_energyphi->Write();
854 
855  //
856  mFileOut->Close();
857  return 1;
858 }
859 Float_t Efficiency::getWeightPions(Float_t pT)
860 {
861  if(pT<0.2) pT=0.2;
862  float weight=1.;
863  if(isPP05){
864  float p[]={592.,-9.31,5.03,-7.55,-1.7,6.06};
865  float WS=1. - 1./(1.+TMath::Exp(TMath::Abs(p[4])*(pT-p[5])));
866  weight=WS*p[2]*pow(pT,p[3]);
867  weight+=(1.- WS)*p[0]*pow(1.+pT,p[1]);
868  weight*=pT;
869  weight=7.0e+05*pT*pow(1.+pT,-9.3);
870  }
871  else if(isDAU){
872  float p[]={39.3,-8.58,0.828,-7.27,0.835,5.18};
873  float WS=1. - 1./(1.+exp(p[4]*(pT-p[5])));
874  weight=WS*p[2]*pow(1.+pT,p[3]);
875  weight+=(1.- WS)*p[0]*pow(1.+pT,p[1]);
876  weight*=pT;
877  weight=7.0e+05*pT*pow(1.+pT,-9.3);
878  }
879 
880  if(!USEWEIGHT) weight=1.;
881  if(USEPYTHIAWEIGHT) weight=exp(-0.401*pT);
882  return weight;
883 }
884 Float_t Efficiency::getWeightEtas(Float_t pT)
885 {
886  Float_t ME=0.548;//eta meson mass
887  Float_t MP=0.135;//pi0 mass
888  Float_t weight=0.45*pT/sqrt(pT*pT+ME*ME-MP*MP)*getWeightPions(sqrt(pT*pT+ME*ME-MP*MP));
889  if(!USEWEIGHT) weight=1.;
890  return weight;
891 }
892 Float_t Efficiency::getWeightAntiNeutrons(Float_t pT)
893 {
894  //hadron data for Levy function, nucl-ex/0601033:
895  //--> d^2N/(2*pi*pT*N*dpT*dy) = B/((1+((mT - m0)/nT))^n)
896  // {p-dAu; pbar-dAu; p-pp; pbar-pp} and m0 = m_neutron = 1.0 GeV.
897  //Double_t B[]={0.3,0.23,0.072,0.061};//->[0]
898  //Double_t eB[]={0.01,0.01,0.005,0.005};
899  Double_t T[]={0.205,0.215,0.179,0.173};//->[1]
900  //Double_t eT[]={0.004,0.005,0.006,0.006};
901  Double_t n[]={11.00,12.55,10.87,10.49};//->[2]
902  //Double_t en[]={0.29,0.41,0.43,0.40};
903 
904  if(pT<0.5) pT=0.5;
905  Float_t weight=0.;
906  Int_t bin=-1;
907  if(isDAU) bin=1;
908  if(isPP05) bin=3;
909  weight=pT/pow((1.+(sqrt(pT*pT+1.) - 1.)/(T[bin]*n[bin])),n[bin]);
910  if(isPP05 && ANTINEUTRONS) weight*=exp(0.5*pT);
911  if(KZEROLONG) weight=getWeightPions(pT);
912  return weight;
913 }
914 Float_t Efficiency::getWeightNeutrons(Float_t pT)
915 {
916  //hadron data for Levy function, nucl-ex/0601033:
917  //--> d^2N/(2*pi*pT*N*dpT*dy) = B/((1+((mT - m0)/nT))^n)
918  // {p-dAu; pbar-dAu; p-pp; pbar-pp} and m0 = m_neutron = 1.0 GeV.
919  //Double_t B[]={0.3,0.23,0.072,0.061};//->[0]
920  //Double_t eB[]={0.01,0.01,0.005,0.005};
921  Double_t T[]={0.205,0.215,0.179,0.173};//->[1]
922  //Double_t eT[]={0.004,0.005,0.006,0.006};
923  Double_t n[]={11.00,12.55,10.87,10.49};//->[2]
924  //Double_t en[]={0.29,0.41,0.43,0.40};
925 
926  if(pT<0.5) pT=0.5;
927  Float_t weight=0.;
928  Int_t bin=-1;
929  if(isDAU) bin=0;
930  if(isPP05) bin=2;
931  weight=pT/pow((1.+(sqrt(pT*pT+1.) - 1.)/(T[bin]*n[bin])),n[bin]);
932  return weight;
933 }
934 Float_t Efficiency::getWeightPhotons(Float_t pT)
935 {
936  if(pT<0.2) pT=0.2;
937  return getWeightPions(pT);
938 }
939 Float_t Efficiency::getWeightVertex(Float_t z)
940 {
941  Float_t weight=1.;
942  if(isDAU){
943  weight=1.04 + 0.00547*z - 0.00016*z*z - 2.8e-06*z*z*z;
944  weight=weight + 4.7e-08*z*z*z*z + 4.5e-10*z*z*z*z*z;
945  }
946  if(isPP05){
947  float sigma_data2=50.*50.;
948  float sigma_mc2=60.*60.;
949  float factor=(1./sigma_data2)-(1./sigma_mc2);
950  weight=exp(-0.5*z*z*factor);
951  }
952  return weight;
953 }
MyPoint
Definition: MyPoint.h:7
Generated by   doxygen 1.8.5