StRoot  1
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Groups Pages
Pi0Analysis.cxx
1 #include <stdio.h>
2 #include <TMultiGraph.h>
3 #include <TStyle.h>
4 #include <TMinuit.h>
5 #include <TTree.h>
6 #include <TF1.h>
7 #include <TFile.h>
8 #include <TCanvas.h>
9 #include <TVector3.h>
10 #include <Riostream.h>
11 #include <TMath.h>
12 
13 #include <StEmcPool/StPhotonCommon/MyEvent.h>
14 #include <StEmcPool/StPhotonCommon/MyPoint.h>
15 #include <StEmcPool/StPhotonCommon/MyMcTrack.h>
16 
17 #include "AnaCuts.h"
18 #include "EventMixer.h"
19 
20 #include "Pi0Analysis.h"
21 using namespace std;
22 
23 //temp solution
24 Bool_t MINBIAS=kTRUE;
25 Bool_t EXCLUDERANGE=kTRUE;
26 Bool_t EXCLUDEETA=kFALSE;
27 Bool_t doNEUTRONS=kFALSE;
28 Float_t rightPI=0.22;
29 Float_t leftPI=0.08;//was 0.12
30 Double_t CONSTANT=0.;
31 Double_t combFit(Double_t *x,Double_t *par)
32 {
33  //old settings:
34  //leftPI=0.08;
35  //rightPI=0.24;
36 
37  if( (x[0]>leftPI&&x[0]<rightPI) && EXCLUDERANGE )
38  {
39  TF1::RejectPoint();
40  return 0;
41  }
42  if( (x[0]>0.5&&x[0]<0.65) && EXCLUDEETA )//thesis result: 0.48
43  {
44  TF1::RejectPoint();
45  return 0;
46  }
47 
48 
49  Double_t ret=par[0]*x[0]+par[1]*x[0]*x[0]+par[2]*x[0]*x[0]*x[0]+CONSTANT;
50 
51  return ret;
52 }
53 Double_t sumFit(Double_t *x,Double_t *par)
54 {
55  Double_t ret=par[0]*x[0]+par[1]*x[0]*x[0]+par[2]*x[0]*x[0]*x[0];
56  Double_t w=x[0]-par[4];
57  Double_t s=par[5];
58  Double_t ret2=par[3]*TMath::Exp(-0.5*(w/s)*(w/s));
59  return ret+ret2;
60 }
61 
62 ClassImp(Pi0Analysis)
63 
64 Pi0Analysis::Pi0Analysis(const Char_t *psout, const Char_t *psout2, const Char_t *flag)
65 {
66  gStyle->SetOptDate(0);
67  ps=new TPostScript(psout,-111);
68  ps2=new TPostScript(psout2,-111);
69 
70  c_array=new TObjArray(100);
71 
72 
73  cout<<"DEFAULT CONSTRUCTOR FOR PI0ANALYSIS!!!!"<<endl;
74  cout<<"storing ps in: "<<psout<<endl;
75 
76  //init cuts
77  cuts=new AnaCuts(flag);
78  //cuts->printCuts();
79 
80  //eventmixer
81  mixer=new EventMixer(flag);
82 
83  runPrev=0;
84  startdatePrev=0;
85  starttimePrev=0;
86  numberOfMB=0;
87  nMBinRun=0;
88  numberOfHT1=0;
89  nHT1inRun=0;
90  numberOfHT2=0;
91  nHT2inRun=0;
92  psMB=0;
93  psHT1=0;
94  psHT2=0;
95 
96  psHT1_eff=0;
97  psHT1_eff2=0;
98  nMB_eff=0;
99 
100  isMC=kFALSE;
101  isPythia=kFALSE;
102  isDAU=kFALSE;
103  isPP05=kFALSE;
104  isAUAU200=kFALSE;
105  isHIJING=kFALSE;
106  if(strcmp(flag,"dAu")==0){
107  isDAU=kTRUE;
108  }
109  if(strcmp(flag,"pp05")==0){
110  isPP05=kTRUE;
111  }
112  if(strcmp(flag,"auau200")==0){
113  isAUAU200=kTRUE;
114  }
115 
116  noMINBIAS=kFALSE;
117 
118  WEIGHT=1.;
119  iev_0=0;
120  iev_1=0;
121  iev_2=0;
122 
123 }
124 Pi0Analysis::~Pi0Analysis()
125 {
126  cout<<endl<<"Pi0Analysis destructed!"<<endl<<endl;
127 }
128 Int_t Pi0Analysis::init(const Char_t *output)
129 {
130 
131  mFileOut=new TFile(output,"RECREATE");
132 
133  h_bbcVsTpc=new TH2F("h_bbcVsTpc","BBC vs TPC vertex z",481,-240.5,240.5,480,-240.,240.);
134  h_bbcVsTpcCorr=new TH2F("h_bbcVsTpcCorr","corrected BBC vs TPC vertex z",481,-240.5,240.5,480,-240.,240.);
135  h_bbcRes=new TH1F("h_bbcres","bbc resolution",100,-2.,2.);
136 
137  h_smdeta1=new TH2F("h_smdeta1","n eta strips in cluster vs p_{T}",20,0.,10.,10,0.5,10.5);
138  h_smdphi1=new TH2F("h_smdphi1","n phi strips in cluster vs p_{T}",20,0.,10.,10,0.5,10.5);
139  h_smdeta2=new TH2F("h_smdeta2","eta energy ratio in cluster vs p_{T}",20,0.,10.,100,0.0,10.);
140  h_smdphi2=new TH2F("h_smdphi2","phi energy ratio in cluster vs p_{T}",20,0.,10.,100,0.0,10.);
141 
142  h_ratioMB=new TH1F("ratioMB","E neutral over E total",60,-.2,1.);
143  h_ratioMB->Sumw2();
144  h_ratioHT1=new TH1F("ratioHT1","E neutral over E total",60,-.2,1.);
145  h_ratioHT1->Sumw2();
146  h_ratioHT2=new TH1F("ratioHT2","E neutral over E total",60,-.2,1.);
147  h_ratioHT2->Sumw2();
148  h_vzMB=new TH1F("h_vzMB","z-vertex MB",480,-120.,120.);
149  h_vzMB->Sumw2();
150  h_vzHT1=new TH1F("h_vzHT1","z-vertex HT1",480,-120.,120.);
151  h_vzHT1->Sumw2();
152  h_vzHT2=new TH1F("h_vzHT2","z-vertex HT2",480,-120.,120.);
153  h_vzHT2->Sumw2();
154  h_trigidHT1=new TH1F("h_trigidHT1","id of trigger HT1",4800,0.5,4800.5);
155  h_trigidHT1->Sumw2();
156  h_trigidHT2=new TH1F("h_trigidHT2","id of trigger HT2",4800,0.5,4800.5);
157  h_trigidHT2->Sumw2();
158  h_trigidHT1aftercut=new TH1F("h_trigidHT1aftercut","id of trigger HT1 after cut",4800,0.5,4800.5);
159  h_trigidHT1aftercut->Sumw2();
160  h_trigidHT2aftercut=new TH1F("h_trigidHT2aftercut","id of trigger HT2 after cut",4800,0.5,4800.5);
161  h_trigidHT2aftercut->Sumw2();
162  h_events=new TH1F("h_events","trigger counts: mb=+1,ht1=+2,ht2=+4,sn=+8",21,-0.5,20.5);
163  h_events->Sumw2();
164 
165  h_HT1adc_id=new TH2F("h_HT1adc_id","ht-1 adc vs id",2400,.5,2400.5,800,200.5,1000.5);
166  h_HT2adc_id=new TH2F("h_HT2adc_id","ht-2 adc vs id",2400,.5,2400.5,800,200.5,1000.5);
167 
168  h_etaphi=new TH2F("h_etaphi","eta/phi of neutral points",300,0.0,1.0,1800,-3.14,3.14);
169  h_etaphi->Sumw2();
170  h_rapphi=new TH2F("h_rapphi","rap/phi of neutral points",300,-0.2,1.2,1800,-3.14,3.14);
171  h_rapphi->Sumw2();
172  h_dist=new TH1F("h_dist","distance of point to track",1000,0.0,200.0);
173  h_dist->Sumw2();
174 
175  h_dist2DMB=new TH2F("h_dist2DMB","dist. vs pt MB",cuts->nPtBinsMB,cuts->ptBinsMB.GetArray(),1000,0.0,1000.0);
176  h_dist2DHT1=new TH2F("h_dist2DHT1","dist. vs pt HT1",cuts->nPtBinsHT1,cuts->ptBinsHT1.GetArray(),1000,0.0,1000.0);
177  h_dist2DHT2=new TH2F("h_dist2DHT2","dist. vs pt HT2",cuts->nPtBinsHT2,cuts->ptBinsHT2.GetArray(),1000,0.0,1000.0);
178 
179  h_dist2DMBpions=new TH2F("h_dist2DMBpions","dist. vs pt MB",cuts->nPtBinsMB,cuts->ptBinsMB.GetArray(),1000,0.0,1000.0);
180  h_dist2DHT1pions=new TH2F("h_dist2DHT1pions","dist. vs pt HT1",cuts->nPtBinsHT1,cuts->ptBinsHT1.GetArray(),1000,0.0,1000.0);
181  h_dist2DHT2pions=new TH2F("h_dist2DHT2pions","dist. vs pt HT2",cuts->nPtBinsHT2,cuts->ptBinsHT2.GetArray(),1000,0.0,1000.0);
182 
183  h_asymmMB=new TH2F("h_asymmMB","asymmetry of inv mass pairs MB",cuts->nPtBinsMB,cuts->ptBinsMB.GetArray(),40,0.0,1.0);
184  h_asymmMB->Sumw2();
185  h_asymmHT1=new TH2F("h_asymmHT1","asymmetry of inv mass pairs HT1",cuts->nPtBinsHT1,cuts->ptBinsHT1.GetArray(),40,0.0,1.0);
186  h_asymmHT1->Sumw2();
187  h_asymmHT2=new TH2F("h_asymmHT2","asymmetry of inv mass pairs HT2",cuts->nPtBinsHT2,cuts->ptBinsHT2.GetArray(),40,0.0,1.0);
188  h_asymmHT2->Sumw2();
189  h_asymmMBbg=new TH2F("h_asymmMBbg","asymmetry of inv mass pairs MB (bg)",cuts->nPtBinsMB,cuts->ptBinsMB.GetArray(),40,0.0,1.0);
190  h_asymmMBbg->Sumw2();
191  h_asymmHT1bg=new TH2F("h_asymmHT1bg","asymmetry of inv mass pairs HT1 (bg)",cuts->nPtBinsHT1,cuts->ptBinsHT1.GetArray(),40,0.0,1.0);
192  h_asymmHT1bg->Sumw2();
193  h_asymmHT2bg=new TH2F("h_asymmHT2bg","asymmetry of inv mass pairs HT2 (bg)",cuts->nPtBinsHT2,cuts->ptBinsHT2.GetArray(),40,0.0,1.0);
194  h_asymmHT2bg->Sumw2();
195 
196  h_minvMB=new TH2F("h_minvMB","invariant mass vs pT MB",cuts->nPtBinsMB,cuts->ptBinsMB.GetArray(),cuts->nMinvBinsMB,cuts->mInvLowMB,cuts->mInvHighMB);
197  h_minvMB->Sumw2();
198  h_minvHT1=new TH2F("h_minvHT1","invariant mass vs pT HT1",cuts->nPtBinsHT1,cuts->ptBinsHT1.GetArray(),cuts->nMinvBinsHT1,cuts->mInvLowHT1,cuts->mInvHighHT1);
199  h_minvHT1->Sumw2();
200  h_minvHT2=new TH2F("h_minvHT2","invariant mass vs pT HT2",cuts->nPtBinsHT2,cuts->ptBinsHT2.GetArray(),cuts->nMinvBinsHT2,cuts->mInvLowHT2,cuts->mInvHighHT2);
201  h_minvHT2->Sumw2();
202 
203  h_pionsVsEtaMB=new TH1F("h_pionsVsEta","dN_{#pi^{0}}/dy data",40,-0.5,1.5);
204 
205  h_gammaMB=new TH1F("h_gammaMB","photon yield vs pT MB",cuts->nPtBinsMB,cuts->ptBinsMB.GetArray());
206  h_gammaMB->Sumw2();
207  h_gammaHT1=new TH1F("h_gammaHT1","photon yield vs pT HT1",cuts->nPtBinsHT1,cuts->ptBinsHT1.GetArray());
208  h_gammaHT1->Sumw2();
209  h_gammaHT2=new TH1F("h_gammaHT2","photon yield vs pT HT2",cuts->nPtBinsHT2,cuts->ptBinsHT2.GetArray());
210  h_gammaHT2->Sumw2();
211 
212  h_pythiaPartonPt=new TH1F("h_pythiaPartonPt","pT of pythia process",200,0.,50);
213  h_pythiaPartonPt->Sumw2();
214  h_pythiaPions=new TH1F("h_pythiaPions","pythia pions vs pT",200,0.,50.);
215  h_pythiaPions->Sumw2();
216  h_pythiaPionsMB=new TH1F("h_pythiaPionsMB","pythia pions vs pT MB",cuts->nPtBinsMB,cuts->ptBinsMB.GetArray());
217  h_pythiaPionsMB->Sumw2();
218  h_pythiaPionsHT1=new TH1F("h_pythiaPionsHT1","pythia pions vs pT HT1",cuts->nPtBinsHT1,cuts->ptBinsHT1.GetArray());
219  h_pythiaPionsHT1->Sumw2();
220  h_pythiaPionsHT2=new TH1F("h_pythiaPionsHT2","pythia pions vs pT HT2",cuts->nPtBinsHT2,cuts->ptBinsHT2.GetArray());
221  h_pythiaPionsHT2->Sumw2();
222 
223  isHot=(TArrayI)cuts->isHot;
224 
225  h_adcHT1=new TH1F("h_adcHT1","highest adc en>0. for HT1 after cuts",1000,0.,1000.);
226  h_adcHT2=new TH1F("h_adcHT2","highest adc en>0. for HT2 after cuts",1000,0.,1000.);
227 
228  h_mcneutronsMB=new TH1F("h_mcneutronsMB","gen. neutrons vs pT MB",cuts->nPtBinsEffMB,cuts->ptBinsEffMB.GetArray());
229  h_mcneutronsHT1=new TH1F("h_mcneutronsHT1","gen. neutrons vs pT HT1",cuts->nPtBinsEff,cuts->ptBinsEff.GetArray());
230  h_mcneutronsHT2=new TH1F("h_mcneutronsHT2","gen. neutrons vs pT HT2",cuts->nPtBinsEff,cuts->ptBinsEff.GetArray());
231  h_mcneutronsWeightMB=new TH1F("h_mcneutronsWeightMB","gen. neutrons vs pT (weighted) MB",cuts->nPtBinsEff,cuts->ptBinsEff.GetArray());
232  h_mcneutronsWeightHT1=new TH1F("h_mcneutronsWeightHT1","gen. neutrons vs pT (weighted) HT1",cuts->nPtBinsEff,cuts->ptBinsEff.GetArray());
233  h_mcneutronsWeightHT2=new TH1F("h_mcneutronsWeightHT2","gen. neutrons vs pT (weighted) HT2",cuts->nPtBinsEff,cuts->ptBinsEff.GetArray());
234  h_neutronsMB=new TH1F("h_neutronsMB","reco. neutrons vs pT MB",cuts->nPtBinsEff,cuts->ptBinsEff.GetArray());
235  h_neutronsHT1=new TH1F("h_neutronsHT1","reco. neutrons vs pT HT1",cuts->nPtBinsEff,cuts->ptBinsEff.GetArray());
236  h_neutronsHT2=new TH1F("h_neutronsHT2","reco. neutrons vs pT HT2",cuts->nPtBinsEff,cuts->ptBinsEff.GetArray());
237 
238  h_EvsE=new TH2F("h_EvsE","E neutral vs E in TPC",80,0.,40,80,0.,40.);
239 
240  h_nstripsETA=new TH2F("h_nstripsETA","n strips vs pt bsmde",15,0.,15.,6,0.,6.);
241  h_nstripsPHI=new TH2F("h_nstripsPHI","n strips vs pt bsmdp",15,0.,15.,6,0.,6.);
242 
243  h_clusterWidth=new TH2F("h_clusterWidth","width BSMD eta+phi",20,0.,20.,100,0.,0.05);
244  h_clusterWidth->Sumw2();
245  h_energyRatio=new TH2F("h_energyRatio","energy ratio BSMD/BTOW",20,0.,20.,160,0.,8.);
246  h_energyRatio->Sumw2();
247 
248  if(isPP05){
249  fout_mb.open("timestamps/pp_timestamps_mb.list");
250  fout_ht1.open("timestamps/pp_timestamps_ht1.list");
251  fout_ht2.open("timestamps/pp_timestamps_ht2.list");
252  }
253  if(isDAU){
254  fout_mb.open("timestamps/dau_timestamps_mb.list");
255  fout_ht1.open("timestamps/dau_timestamps_ht1.list");
256  fout_ht2.open("timestamps/dau_timestamps_ht2.list");
257  }
258 
259  return 1;
260 }
261 Int_t Pi0Analysis::make(Int_t evmax, const Char_t *input)
262 {
263  mFile=new TFile(input,"OPEN");
264  myEventTree=(TTree*)mFile->Get("mEventTree");
265  ev=new MyEvent();
266  myEventTree->SetBranchAddress("branch",&ev);
267 
268  Int_t i=0;
269  while(myEventTree->GetEntry(i)){
270  if(evmax&&i>evmax){
271  cout<<"reached evmax,"<<endl;
272  cout<<"abort loop!"<<endl;
273  break;
274  }
275  if(i%500000==0) cout<<"processing "<<input<<" evt: "<<i<<endl;
276 
277  if(i==0){
278  startdatePrev=ev->date();
279  starttimePrev=ev->time();
280  }
281 
282  //assign trigger for MC event:
283  Int_t mcTr=1;//mb
284  if(isMC&&!isPythia){
285  if(ev->highTowerAdc()>cuts->ht1AdcCUT) mcTr+=2;
286  if(ev->highTowerAdc()>cuts->ht2AdcCUT) mcTr+=4;
287  ev->setTrigger(mcTr);
288  }
289 
290  TVector3 vPos=ev->vertex();
291 
292  float bbc_vertz=(ev->bbcVertexZ()-10.77)/2.82;
293  //now bbc_vert is west-east, then from correlation:
294  bbc_vertz+=4.04;
295  bbc_vertz/=0.354;
296  ev->setBbcVertexZ(bbc_vertz);
297 
298  if(ev->trigger()&1 && TMath::Abs(vPos.Z())>0.00000001){
299  h_bbcVsTpc->Fill(vPos.Z(),ev->bbcVertexZ());
300  }
301 
302  if(ev->trigger()&1) {h_vzMB->Fill(vPos.Z());}
303  if(ev->trigger()&2) {h_vzHT1->Fill(vPos.Z());h_trigidHT1->Fill(ev->highTowerId());}
304  if(ev->trigger()&4) {h_vzHT2->Fill(vPos.Z());h_trigidHT2->Fill(ev->highTowerId());}
305 
306  if(ev->trigger()&4 && TMath::Abs(vPos.Z())<60.){h_EvsE->Fill(ev->momentumInTpc(),ev->energyInBarrel());}
307 
308  //set tpc momentum for simu:
309  if(isMC) ev->setMomentumInTpc(99999.);
310 
311  //exclude depleted minbias
312  if(noMINBIAS && ev->trigger()&1){
313  int trigg=ev->trigger();
314  ev->setTrigger(trigg-1);
315  }
316 
317  //event cuts + hot tower:
318  if(isDAU && !cuts->isEventOK(ev,"dAu")) {i++;continue;}
319  if(isPP05 && !cuts->isEventOK(ev,"pp05")) {i++;continue;}
320 
321 
322 
323  //fill mixer
324  //if(isDAU || isPP05) mixer->addEvent(ev);
325 
326  //check resolution bbc:
327  if(isPP05 && TMath::Abs(vPos.Z())>0.){
328  float bbcres=(ev->bbcVertexZ()-vPos.Z())/vPos.Z();
329  h_bbcRes->Fill(bbcres);
330  }
331 
332  Float_t ratioE=TMath::Abs(ev->energyInBarrel()+ev->momentumInTpc())>0. ? ev->energyInBarrel()/(ev->energyInBarrel()+ev->momentumInTpc()) : -.1;
333 
334  if(i>0 && ev->runId()!=runPrev){
335  psHT1_eff+=psHT1*nHT1inRun;
336  nMB_eff+=psMB*nMBinRun;
337 
338  fout_mb<<runPrev<<" "<<startdatePrev<<" "<<starttimePrev<<" "<<nMBinRun<<endl;
339  fout_ht1<<runPrev<<" "<<startdatePrev<<" "<<starttimePrev<<" "<<nHT1inRun<<endl;
340  fout_ht2<<runPrev<<" "<<startdatePrev<<" "<<starttimePrev<<" "<<nHT2inRun<<endl;
341 
342  nMBinRun=0;
343  nHT1inRun=0;
344  nHT2inRun=0;
345  startdatePrev=ev->date();
346  starttimePrev=ev->time();
347  }
348  if(ev->trigger()&1) {nMBinRun++;}
349  if(ev->trigger()&2) {nHT1inRun++;}
350  if(ev->trigger()&4) {nHT2inRun++;}
351  psMB=(ev->prescale(0) > ev->prescale(1)) ? (Float_t)ev->prescale(0) : (Float_t)ev->prescale(1);
352  psHT1=(Float_t)ev->prescale(2);
353  psHT2=(Float_t)ev->prescale(3);
354 
355  //get neutrons:
356  for(Int_t j=0;j<ev->numberOfMcPions();j++){
357  MyMcTrack *pyt=(MyMcTrack*)ev->getMcPionArray()->At(j);
358  //check if neutron in rapidity bite:
359  //change this to rapidity...
360  if(pyt->momentum().PseudoRapidity()<cuts->rapidityMinCUT||pyt->momentum().PseudoRapidity()>cuts->rapidityMaxCUT) continue;
361  h_mcneutronsMB->Fill(pyt->momentum().Pt());
362  h_mcneutronsHT1->Fill(pyt->momentum().Pt());
363  h_mcneutronsHT2->Fill(pyt->momentum().Pt());
364  Float_t weight=getWeight(pyt->momentum().Pt());
365  h_mcneutronsWeightMB->Fill(pyt->momentum().Pt(),weight);
366  h_mcneutronsWeightHT1->Fill(pyt->momentum().Pt(),weight);
367  h_mcneutronsWeightHT2->Fill(pyt->momentum().Pt(),weight);
368  }
369 
370  if(ev->trigger()&1){
371  h_ratioMB->Fill(ratioE);
372  numberOfMB++;
373  }
374  if(ev->trigger()&2){
375  h_ratioHT1->Fill(ratioE);
376  numberOfHT1++;
377  h_trigidHT1aftercut->Fill(ev->highTowerId());
378 
379  h_adcHT1->Fill(ev->highTowerAdc());
380  h_HT1adc_id->Fill(ev->highTowerId(),ev->highTowerAdc());
381  }
382  if(ev->trigger()&4){
383  h_ratioHT2->Fill(ratioE);
384  numberOfHT2++;
385  h_trigidHT2aftercut->Fill(ev->highTowerId());
386 
387  h_adcHT2->Fill(ev->highTowerAdc());
388  h_HT2adc_id->Fill(ev->highTowerId(),ev->highTowerAdc());
389  }
390 
391  h_events->Fill(ev->trigger());
392 
393  if(isPythia){
394  Float_t pypT=ev->partonPt();
395  if(pypT<15.){
396  WEIGHT=1.;
397  iev_0++;
398  }
399  else if(pypT<25.){
400  WEIGHT=(216000/208000)*(0.0003895/0.002228);
401  iev_1++;
402  }
403  else if(pypT>25.){
404  iev_2++;
405  WEIGHT=(216000/208000)*(1.016e-005/0.002228);
406  }
407  h_pythiaPartonPt->Fill(ev->partonPt(),WEIGHT);
408  for(Int_t it=0;it<ev->getMcPionArray()->GetEntries();it++){
409  MyMcTrack *pyt=(MyMcTrack*)ev->getMcPionArray()->At(it);
410  if(pyt->momentum().PseudoRapidity()<cuts->rapPionMinCUT||pyt->momentum().PseudoRapidity()>cuts->rapPionMaxCUT) continue;
411  h_pythiaPions->Fill(pyt->momentum().Pt(),WEIGHT);
412  if(ev->trigger()&1){
413  h_pythiaPionsMB->Fill(pyt->momentum().Pt(),WEIGHT);
414  }
415  if(ev->trigger()&2){
416  h_pythiaPionsHT1->Fill(pyt->momentum().Pt(),WEIGHT);
417  }
418  if(ev->trigger()&4){
419  h_pythiaPionsHT2->Fill(pyt->momentum().Pt(),WEIGHT);
420  }
421  }
422  }
423 
424  runPrev=ev->runId();
425  //loop on points
426  MyPoint *p;
427  MyPoint *pp;
428  TClonesArray *clA=(TClonesArray*)ev->getPointArray();
429  for(Int_t j=0;j<clA->GetEntries();j++){
430  p=(MyPoint*)clA->At(j);
431  TVector3 pPos=p->position();
432  TVector3 pMom=pPos - vPos;
433  pMom.SetMag(p->energy());
434 
435  if(!cuts->isPointOK(p,vPos)) continue;
436 
437  h_dist->Fill(p->distanceToTrack());
438  if(ev->trigger()&1) h_dist2DMB->Fill(pMom.Pt(),p->distanceToTrack());
439  if(ev->trigger()&2) h_dist2DHT1->Fill(pMom.Pt(),p->distanceToTrack());
440  if(ev->trigger()&4) h_dist2DHT2->Fill(pMom.Pt(),p->distanceToTrack());
441 
442  if(ev->trigger()&1 && pMom.Pt()>1.){
443  h_etaphi->Fill(pPos.PseudoRapidity(),pPos.Phi());
444  h_rapphi->Fill(pMom.PseudoRapidity(),pMom.Phi());
445  }
446 
447  if(pMom.PseudoRapidity()>cuts->rapidityMinCUT && pMom.PseudoRapidity()<cuts->rapidityMaxCUT){
448  //first neutrons:
449  if(doNEUTRONS){
450  Float_t D=1000000.;
451  Float_t pt4weight=10.;
452  for(Int_t j=0;j<ev->numberOfMcPions();j++){
453  MyMcTrack *pyt=(MyMcTrack*)ev->getMcPionArray()->At(j);
454  TVector3 mcPos=pyt->position();
455  Float_t n=mcPos.PseudoRapidity();
456  Float_t p=mcPos.Phi();
457  Float_t nn=pPos.PseudoRapidity();
458  Float_t pp=pPos.Phi();
459  Float_t d=sqrt(pow(p-pp,2)+pow(n-nn,2));
460  if(d<D){
461  D=d;
462  pt4weight=pyt->momentum().Pt();
463  }
464  }
465  if(D<0.025){
466  if(ev->trigger()&1) h_neutronsMB->Fill(pMom.Pt(),getWeight(pt4weight));
467  if(ev->trigger()&2) h_neutronsHT1->Fill(pMom.Pt(),getWeight(pt4weight));
468  if(ev->trigger()&4) h_neutronsHT2->Fill(pMom.Pt(),getWeight(pt4weight));
469  }
470  }
471 
472  //fill photons:
473  if(p->distanceToTrack()>cuts->photonCUT){
474  if(ev->trigger()&1) h_gammaMB->Fill(pMom.Pt(),WEIGHT);
475  if(ev->trigger()&2) h_gammaHT1->Fill(pMom.Pt(),WEIGHT);
476  if(ev->trigger()&4) h_gammaHT2->Fill(pMom.Pt(),WEIGHT);
477 
478  if(ev->trigger()&1){
479  h_clusterWidth->Fill(pMom.Pt(),sqrt(p->widthEta()*p->widthEta()+p->widthPhi()*p->widthPhi()));
480  h_energyRatio->Fill(pMom.Pt(),(p->energyEta()+p->energyPhi())/p->energy());
481  h_nstripsETA->Fill(pMom.Pt(),p->nHitsEta(),WEIGHT);
482  h_nstripsPHI->Fill(pMom.Pt(),p->nHitsPhi(),WEIGHT);
483  //super isolated:
484  if(p->distanceToTrack()>30.){
485  h_smdeta1->Fill(pMom.Pt(),p->nHitsEta());
486  h_smdphi1->Fill(pMom.Pt(),p->nHitsPhi());
487  h_smdeta2->Fill(pMom.Pt(),p->energyEta()/p->energy());
488  h_smdphi2->Fill(pMom.Pt(),p->energyPhi()/p->energy());
489  }
490  }
491  }
492 
493  }
494 
495  for(Int_t jj=0;jj<clA->GetEntries();jj++){
496  if(jj<=j) continue;
497  pp=(MyPoint*)clA->At(jj);
498  TVector3 ppPos=pp->position();
499  TVector3 ppMom=ppPos - vPos;
500  ppMom.SetMag(pp->energy());
501  if(!cuts->isPointOK(pp,vPos)) continue;
502 
503  if(p->distanceToTrack()<pp->distanceToTrack()){
504  if(ev->trigger()&1) h_dist2DMBpions->Fill(pMom.Pt(),p->distanceToTrack());
505  if(ev->trigger()&2) h_dist2DHT1pions->Fill(pMom.Pt(),p->distanceToTrack());
506  if(ev->trigger()&4) h_dist2DHT2pions->Fill(pMom.Pt(),p->distanceToTrack());
507  }
508  else{
509  if(ev->trigger()&1) h_dist2DMBpions->Fill(pMom.Pt(),pp->distanceToTrack());
510  if(ev->trigger()&2) h_dist2DHT1pions->Fill(pMom.Pt(),pp->distanceToTrack());
511  if(ev->trigger()&4) h_dist2DHT2pions->Fill(pMom.Pt(),pp->distanceToTrack());
512  }
513 
514  //two neutrals
515  TVector3 pi0Mom=pMom+ppMom;
516  Float_t angle=pMom.Angle(ppMom);
517  Float_t Epion=p->energy()+pp->energy();
518  Float_t asymm=TMath::Abs(p->energy()-pp->energy())/Epion;
519  Float_t minv=sqrt(2.*p->energy()*pp->energy()*(1. - cos(angle)));
520  Float_t pTpion=pi0Mom.Pt();
521 
522  //select rapidity bite / asymm cut for pions:
523  if(pi0Mom.PseudoRapidity()<cuts->rapPionMinCUT||pi0Mom.PseudoRapidity()>cuts->rapPionMaxCUT) continue;
524 
525  if(ev->trigger()&1){
526  if(asymm<=cuts->asymmetryCUT) h_minvMB->Fill(pTpion,minv,WEIGHT);
527  if(pTpion>1.){
528  if(minv>0.08 && minv<0.20){
529  h_asymmMB->Fill(pTpion,asymm);
530  if(asymm<=cuts->asymmetryCUT) h_pionsVsEtaMB->Fill(pi0Mom.PseudoRapidity());
531  }
532  else h_asymmMBbg->Fill(pTpion,asymm);
533  }
534  }
535  if(ev->trigger()&2){
536  if(asymm<=cuts->asymmetryCUT) h_minvHT1->Fill(pTpion,minv,WEIGHT);
537  if(pTpion>4.){
538  if(minv>0.10 && minv<0.20){
539  h_asymmHT1->Fill(pTpion,asymm);
540  }
541  else h_asymmHT1bg->Fill(pTpion,asymm);
542  }
543  }
544  if(ev->trigger()&4){
545  if(asymm<=cuts->asymmetryCUT) h_minvHT2->Fill(pTpion,minv,WEIGHT);
546  if(pTpion>6.){
547  if(minv>0.10 && minv<0.2){
548  h_asymmHT2->Fill(pTpion,asymm);
549  }
550  else h_asymmHT2bg->Fill(pTpion,asymm);
551  }
552  }
553 
554  }
555 
556 
557  }
558 
559 
560  i++;
561  }
562 
563  return 1;
564 }
565 
566 void Pi0Analysis::printPrescales()
567 {
568  cout<<"--------------------------------"<<endl;
569  cout<<" stats processed: "<<endl;
570  cout<<" "<<endl;
571  cout<<" number of mb: "<<numberOfMB<<endl;
572  cout<<" number of ht1: "<<numberOfHT1<<endl;
573  cout<<" number of ht2: "<<numberOfHT2<<endl;
574  cout<<"--------------------------------"<<endl;
575  cout<<" prescale corrections: "<<endl;
576  cout<<" "<<endl;
577  cout<<" psHT1_eff "<<psHT1_eff/(1.*numberOfHT1)<<endl;
578  cout<<endl;
579  cout<<" N_mb*ps_mb: "<<nMB_eff<<endl;
580  cout<<"--------------------------------"<<endl;
581 }
582 Int_t Pi0Analysis::finish()
583 {
584  h_minvMB_mixed=new TH2F(*(TH2F*)mixer->getMinvMB());
585  mFileOut->Add(h_minvMB_mixed);
586  mFileOut->Write();
587 
588  return 1;
589 }
590 
591 void Pi0Analysis::getYield()
592 {
593  h_yieldMB=new TH1F(*getYield(NULL,"mb"));
594  h_yieldHT1=new TH1F(*getYield(NULL,"ht1"));
595  h_yieldHT2=new TH1F(*getYield(NULL,"ht2"));
596 }
597 
598 //new routine:
599 TH1F *Pi0Analysis::getYield(TH2F *h, const Char_t *flag)
600 {
601  ps->On();
602  if(!h){
603  cout<<"zero pointer"<<endl;
604  if(strcmp(flag,"mb")==0) h=new TH2F(*h_minvMB);
605  if(strcmp(flag,"ht1")==0) h=new TH2F(*h_minvHT1);
606  if(strcmp(flag,"ht2")==0) h=new TH2F(*h_minvHT2);
607  }
608  TString retname(h->GetName());
609  retname.Append("_yield");
610  TString rettitle(h->GetTitle());
611  rettitle.Append("_yield");
612  TH1F *retYield=new TH1F(*((TH1F*)h->ProjectionX()));
613  retYield->Reset();
614  retYield->Sumw2();
615  retYield->SetNameTitle(retname.Data(),rettitle.Data());
616 
617  Char_t canvasname[100];
618  int page=0;
619  TCanvas *c_copy;
620 
621  Float_t NSIGMAHI=3.;
622  Float_t NSIGMALO=3.;
623 
624  const int nx=h->GetNbinsX();
625  Float_t *x=new Float_t[nx];
626  Float_t *ex=new Float_t[nx];
627  Float_t *y=new Float_t[nx];
628  Float_t *ey=new Float_t[nx];
629  Float_t *s=new Float_t[nx];
630  Float_t *es=new Float_t[nx];
631  Float_t *m=new Float_t[nx];
632  Float_t *em=new Float_t[nx];
633 
634  cout<<"------------- "<<flag<<" --------------"<<endl;
635 
636  Int_t color=1;
637  Int_t cbgcolor=1;
638  Int_t peakcolor=1;
639  Int_t fillcolor=16;
640  MINBIAS=kTRUE;
641  if(strcmp(flag,"ht1")==0) {color=4; cbgcolor=2; peakcolor=1; fillcolor=38; MINBIAS=kFALSE;}
642  if(strcmp(flag,"ht2")==0) {color=2; cbgcolor=4; peakcolor=1; fillcolor=46; MINBIAS=kFALSE;}
643 
644  TCanvas *c=new TCanvas(flag,flag,600,600);
645  c->Divide(3,3);
646 
647  TF1 *combbg[50];
648  Char_t combname[100];
649 
650  TF1 *pi0peak[50];
651  Char_t peakname[100];
652  TF1 *sumfit[50];
653  Char_t sumname[100];
654 
655  TH1D *ptslice[50];
656  Char_t name[100];
657 
658  TH1D *ptslice_bg[50];
659 
660  TH1D *ptslice_bgsub[50];
661  Char_t subname[100];
662  Char_t title[100];
663 
664  for(Int_t i=1;i<=nx;i++){
665  if((i-1)%9==0){
666  if(i>1){
667  c->Update();
668  sprintf(canvasname,"c_%s_page_%d",flag,page++);
669  c_copy=(TCanvas*)c->Clone();
670  c_copy->SetName(canvasname);
671  c_array->Add(c_copy);
672  }
673  ps->NewPage();
674  }
675  c->cd((i-1)%9 + 1);
676 
677  sprintf(combname,"comb_%d_%s",i,flag);
678 
679  sprintf(peakname,"gaus_%d_%s",i,flag);
680  sprintf(sumname,"sum_%d_%s",i,flag);
681  sprintf(name,"ptslice_%d_%s",i,flag);
682  sprintf(subname,"subslice_%d_%s",i,flag);
683  sprintf(title,"%.2f < p_{T} < %.2f",h->GetXaxis()->GetBinLowEdge(i),h->GetXaxis()->GetBinUpEdge(i));
684 
685  //get bin
686  ptslice[i]=h->ProjectionY(name,i,i,"e");
687 
688  //for next iter.
689  ptslice_bg[i]=new TH1D(*ptslice[i]);
690 
691 
692  //set fitfunction
693  combbg[i]=new TF1(combname,combFit,0.0,1.0,3);
694 
695 
696  sumfit[i]=new TF1(sumname,sumFit,0.0,1.0,6);
697 
698 
699  float firstbin=0;
700  for(int ib=1;ib<=ptslice[i]->GetNbinsX();ib++){
701  if(ptslice[i]->GetBinContent(ib)>0.){
702  firstbin=(float)ptslice[i]->GetXaxis()->GetBinCenter(ib);
703  break;
704  }
705  }
706 
707  combbg[i]->SetParameters(1.0,-1.0,1.0);
708 
709  EXCLUDERANGE=kTRUE;//leave out pi0 range, interval:
710  EXCLUDEETA=kTRUE;
711  if(strcmp(flag,"mb")==0) rightPI=isDAU ? 0.2 : 0.2;//dau was 0.26 : 0.22
712  if(strcmp(flag,"ht1")==0) rightPI=isDAU ? 0.281 : 0.281;
713  if(strcmp(flag,"ht2")==0) rightPI=isDAU ? 0.26 : 0.26;
714 
715  if(isDAU && strcmp(flag,"ht1")==0){
716  if(h->GetXaxis()->GetBinCenter(i)>9.){
717  combbg[i]->FixParameter(2,0.);//was 9.
718  //rightPI=0.3;
719  }
720  }
721  if(isDAU && strcmp(flag,"ht2")==0){
722  if(h->GetXaxis()->GetBinCenter(i)>11.){//was 11.
723  combbg[i]->FixParameter(2,0.);//was 9.
724  //rightPI=0.3;
725  }
726  }
727 
728  if(isPP05 && h->GetXaxis()->GetBinCenter(i)>9.) combbg[i]->FixParameter(2,0.);
729 
730  combbg[i]->SetRange(0.0,0.8);
731  if(isPP05){
732  combbg[i]->SetRange(0.28,0.8);//was 0.0-0.8 (24/01/08)
733  if(h->GetXaxis()->GetBinCenter(i)<7. && strcmp(flag,"ht2")==0){
734  combbg[i]->SetRange(0.3,0.8);
735  }
736  }
737  if(isDAU) combbg[i]->SetRange(0.28,0.8);
738 
739  sumfit[i]->SetRange(0.0,0.8);//was 0.0-0.8, this to show background
740 
741 
742  //MARKER
743  if(isPP05 && strcmp(flag,"mb")==0){
744  combbg[i]->SetRange(0.0,0.4);//0.8
745  if(h->GetXaxis()->GetBinCenter(i)<2.){
746  combbg[i]->SetRange(0.05,0.4);
747  leftPI=0.081;
748  rightPI=0.19;
749  }
750  if(h->GetXaxis()->GetBinCenter(i)>2.){
751  combbg[i]->FixParameter(2,0.);
752  combbg[i]->SetRange(0.05,0.4);
753  leftPI=0.10;
754  rightPI=0.22;
755  }
756  if(h->GetXaxis()->GetBinCenter(i)>3.){
757  leftPI=0.06;
758  rightPI=0.22;
759  }
760  sumfit[i]->SetRange(0.0,0.4);
761  }
762  if(isDAU && strcmp(flag,"mb")==0){
763  combbg[i]->SetRange(0.0,0.4);//0.8
764  if(h->GetXaxis()->GetBinCenter(i)<2.){
765  combbg[i]->SetRange(0.0,0.4);
766  leftPI=0.0;
767  rightPI=0.2;
768  }
769  if(h->GetXaxis()->GetBinCenter(i)>2.){
770  leftPI=0.0;
771  rightPI=0.22;
772  }
773  if(h->GetXaxis()->GetBinCenter(i)>4.){
774  combbg[i]->FixParameter(2,0.);
775  leftPI=0.0;
776  rightPI=0.22;
777  }
778  sumfit[i]->SetRange(0.0,0.4);
779  }
780 
781 
782  if(!isMC||isPythia) ptslice[i]->Fit(combname,"QR0");
783  EXCLUDERANGE=kFALSE;
784  EXCLUDEETA=kFALSE;
785  combbg[i]->SetRange(0.0,0.8);
786  if(strcmp(flag,"mb")==0) combbg[i]->SetRange(0.0,0.4);
787  //cout<<gMinuit->fCstatu.Data()<<endl;
788 
789  ptslice[i]->SetMarkerStyle(4);
790  ptslice[i]->SetMarkerSize(.7);
791  ptslice[i]->SetMarkerColor(color);
792  ptslice[i]->SetLineColor(color);
793  ptslice[i]->Draw();
794  combbg[i]->SetLineStyle(2);
795  combbg[i]->SetLineWidth(1);
796  combbg[i]->SetLineColor(cbgcolor);
797  sumfit[i]->SetLineWidth(1);
798  sumfit[i]->SetLineColor(cbgcolor);
799 
800 
801  ptslice_bgsub[i]=new TH1D(*ptslice[i]);
802  ptslice_bgsub[i]->SetName(subname);
803  if(!isMC||isPythia){
804  ptslice_bgsub[i]->Add(combbg[i],-1.0);
805  combbg[i]->Draw("same");
806  }
807 
808  //fit peak
809  pi0peak[i]=new TF1(peakname,"gaus",0.05,0.2);
810  pi0peak[i]->SetParLimits(1,.1,.2);
811 
812  if(!isMC||isPythia) ptslice_bgsub[i]->Fit(peakname,"QR0");
813  else ptslice[i]->Fit(peakname,"QR0");
814  pi0peak[i]->SetLineColor(peakcolor);
815  pi0peak[i]->SetLineStyle(3);
816  pi0peak[i]->SetLineWidth(2);
817  pi0peak[i]->SetRange(0.0,0.8);
818 
819  //draw sum of comb. + sign.
820  Double_t params[6];
821  combbg[i]->GetParameters(&params[0]);
822  pi0peak[i]->GetParameters(&params[3]);
823  sumfit[i]->SetParameters(params);
824  if(!isMC||isPythia){
825  sumfit[i]->Draw("same");
826  }
827 
828  ptslice[i]->SetTitle(title);
829  TAxis *ax=ptslice[i]->GetXaxis();
830  ax->SetRangeUser(0.0,.8);
831  if(strcmp(flag,"mb")==0) ax->SetRangeUser(0.0,0.4);
832 
833  if(!isMC||isPythia){
834  Double_t MINIMUM=ptslice_bgsub[i]->GetMinimum();
835  ptslice[i]->SetMinimum(MINIMUM);
836  }
837 
838  //get 1st yield
839  Float_t mean=pi0peak[i]->GetParameter(1);
840  Float_t meanerr=pi0peak[i]->GetParError(1);
841  Float_t sigma=pi0peak[i]->GetParameter(2);
842  Float_t sigmaerr=pi0peak[i]->GetParError(2);
843 
844  Float_t yield=pi0peak[i]->Integral(mean-2.*sigma,mean+2.*sigma);
845  yield/=ptslice_bgsub[i]->GetBinWidth(1);
846  Float_t erryield=0.;
847 
848  if(yield>0.) erryield=sqrt(yield);
849  else {yield=0.;erryield=0.;}
850 
851  x[i-1]=h->GetXaxis()->GetBinCenter(i);
852  ex[i-1]=0.;
853  y[i-1]=yield/h->GetXaxis()->GetBinWidth(i);// divide dpt
854  ey[i-1]=erryield/h->GetXaxis()->GetBinWidth(i);
855  m[i-1]=mean;
856  em[i-1]=meanerr;
857  s[i-1]=sigma;
858  es[i-1]=sigmaerr;
859 
860  }
861  c->Update();
862  sprintf(canvasname,"c_%s_page_%d",flag,page++);
863  c_copy=(TCanvas*)c->Clone();
864  c_copy->SetName(canvasname);
865  c_array->Add(c_copy);
866  ps->NewPage();
867 
868  Float_t pTmin=1.;
869  Float_t pTmax=7.;
870  Float_t thres=0;
871  if(strcmp(flag,"ht1")==0) {pTmin=3.;pTmax=11.;thres=2.5;}
872  if(strcmp(flag,"ht2")==0) {pTmin=5.;pTmax=15.;thres=4.5;}
873  //fit mean
874  TF1 *meanfit=new TF1("meanfit","[0]*(1.0 - exp(-[1]*(x-[2])))");
875  meanfit->SetParameter(0,.145);
876  meanfit->SetParameter(1,2.);
877  meanfit->SetParameter(2,thres);
878  meanfit->SetParLimits(2,thres-2.0,10.);
879  meanfit->SetRange(pTmin,pTmax);
880  meanfit->SetLineColor(color);
881  meanfit->SetLineStyle(2);
882  //fit sigma
883  TF1 *sigmafit=new TF1("sigmafit","[0]+[1]*exp(-[2]*(x-[3]))");
884  sigmafit->SetParameter(0,0.02);
885  sigmafit->SetParameter(1,0.05);
886  sigmafit->SetParameter(2,1.0);
887  sigmafit->SetParameter(3,thres);
888  sigmafit->SetParLimits(3,thres-2.0,10.);
889  sigmafit->SetRange(pTmin,pTmax);
890  sigmafit->SetLineColor(color);
891  sigmafit->SetLineStyle(2);
892 
893  TCanvas *c2=new TCanvas(TString("c2_")+TString(flag),TString("c2_")+TString(flag),600,600);
894  c2->Divide(2,2);
895  c2->cd(1);
896 
897  //draw mean,sigma,yield vs pT
898  TGraphErrors *gry=new TGraphErrors(nx,x,y,ex,ey);
899  gry->SetName((TString(flag)+"_yieldfit").Data());
900  gPad->SetLogy();
901  gry->SetMarkerStyle(25);
902  gry->SetMarkerSize(.9);
903  gry->SetMarkerColor(color);
904  gry->GetXaxis()->SetRangeUser(pTmin,pTmax);
905  gry->Draw("ap");
906  gry->SetMaximum(10.*gry->GetHistogram()->GetMaximum());
907  Float_t minimum=.1*gry->GetHistogram()->GetMinimum();
908  gry->SetMinimum(minimum > 1. ? minimum : 1.);
909 
910  c2->cd(2);
911  TGraphErrors *grm=new TGraphErrors(nx,x,m,ex,em);
912  grm->SetName((TString(flag)+"_meanfit").Data());
913  grm->SetMarkerStyle(25);
914  grm->SetMarkerSize(.9);
915  grm->SetMarkerColor(color);
916  grm->GetXaxis()->SetRangeUser(pTmin,pTmax);
917  grm->Draw("ap");
918  grm->Fit("meanfit","QR");
919  meanfit->Draw("same");
920  c2->cd(3);
921  TGraphErrors *grs=new TGraphErrors(nx,x,s,ex,es);
922  grs->SetName((TString(flag)+"_sigmafit").Data());
923  grs->SetMarkerStyle(25);
924  grs->SetMarkerSize(.9);
925  grs->SetMarkerColor(color);
926  grs->SetMaximum(0.15);
927  grs->SetMinimum(0.);
928  grs->GetXaxis()->SetRangeUser(pTmin,pTmax);
929  grs->Draw("ap");
930  grs->Fit("sigmafit","QR");
931  c2->cd(0);
932 
933  //c2->Update();
934  //ps->NewPage();
935  delete c2;
936 
937  TCanvas *cr=new TCanvas(TString("cr_")+TString(flag),TString("cr_")+TString(flag),1000,1000);
938  cr->Divide(4,5);
939 
940  //re-integrate peak region (with correct error..):
941  Float_t *x2=new Float_t[nx];
942  Float_t *ex2=new Float_t[nx];
943  Float_t *y2=new Float_t[nx];
944  Float_t *ey2=new Float_t[nx];
945  Float_t *m2=new Float_t[nx];
946  Float_t *em2=new Float_t[nx];
947  Float_t *s2=new Float_t[nx];
948  Float_t *es2=new Float_t[nx];
949  Float_t *bgy=new Float_t[nx];
950 
951 
952  TF1 *finalpeak[100];
953  TH1D *intregion[100];
954  Char_t newname[100];
955  Char_t newname2[100];
956  Char_t newtitle[100];
957  Char_t newfit[100];
958  for(Int_t g=1;g<=nx;g++)
959  {
960  sprintf(newname,"ptbin_%d_%s",g,flag);
961  sprintf(newname,"ptbin2_%d_%s",g,flag);
962  sprintf(newtitle,"%.2f < p_{T} < %.2f",h->GetXaxis()->GetBinLowEdge(g),h->GetXaxis()->GetBinUpEdge(g));
963  sprintf(newfit,"fit_peak_%d_%s",g,flag);
964  finalpeak[g]=new TF1(newfit,"gaus");
965  //get pT
966  Float_t xval=h->GetXaxis()->GetBinCenter(g);
967  //limit peak to +/- 2 sigma from extrapolation
968  Float_t left=meanfit->Eval(xval)-2.*sigmafit->Eval(xval);
969  Float_t right=meanfit->Eval(xval)+2.*sigmafit->Eval(xval);
970 
971  finalpeak[g]->SetParLimits(1,left,right);
972  //define fit range +/- 3*sigma
973  left=meanfit->Eval(xval)-3*sigmafit->Eval(xval);
974  right=meanfit->Eval(xval)+3*sigmafit->Eval(xval);
975 
976  finalpeak[g]->SetRange(left,right);
977  finalpeak[g]->SetLineColor(color);
978  finalpeak[g]->SetLineStyle(2);
979  finalpeak[g]->SetLineWidth(1);
980  cr->cd(g);
981  ptslice_bgsub[g]->SetName(newname);
982  ptslice_bgsub[g]->SetTitle(newtitle);
983  ptslice_bgsub[g]->GetXaxis()->SetTitle("M_{inv} (GeV/c^{2})");
984  ptslice_bgsub[g]->GetYaxis()->SetTitle("counts (10 MeV/c^{2})^{-1}");
985  if(strcmp(flag,"mb")!=0){
986  ptslice_bgsub[g]->GetYaxis()->SetTitle("counts (20 MeV/c^{2})^{-1}");
987  }
988  ptslice_bgsub[g]->SetFillColor(0);
989  ptslice_bgsub[g]->SetLineColor(color);
990  ptslice_bgsub[g]->Draw("hist");
991  ptslice_bgsub[g]->GetXaxis()->SetRangeUser(0.0,.8);
992  if(strcmp(flag,"mb")==0){
993  ptslice_bgsub[g]->GetXaxis()->SetRangeUser(0.0,0.4);
994  //fit residual:
995  }
996 
997  ptslice_bgsub[g]->Fit(newfit,"QR0");
998  finalpeak[g]->SetRange(0.0,1.0);
999  finalpeak[g]->Draw("same");
1000 
1001 
1002  //get integration region
1003  intregion[g]=new TH1D(*ptslice_bgsub[g]);
1004  intregion[g]->SetName(newname2);
1005  intregion[g]->SetFillColor(fillcolor);
1006  Float_t LEFT=finalpeak[g]->GetParameter(1)-NSIGMALO*finalpeak[g]->GetParameter(2);
1007  Float_t RIGHT=finalpeak[g]->GetParameter(1)+(NSIGMAHI)*finalpeak[g]->GetParameter(2);
1008  if(isMC) RIGHT=finalpeak[g]->GetParameter(1)+(NSIGMAHI+3)*finalpeak[g]->GetParameter(2);
1009  if(!isMC){
1010  if(isDAU){
1011  if(strcmp(flag,"mb")==0){
1012  if(RIGHT>0.2) RIGHT=0.2;
1013  if(LEFT<0.08) LEFT=0.081;
1014  }
1015  if(strcmp(flag,"ht1")==0&&xval>9.){
1016  LEFT=0.081;
1017  RIGHT=0.26;
1018  }
1019  if(strcmp(flag,"ht2")==0&&xval>9.){
1020  LEFT=0.081;
1021  if(xval>12.) RIGHT=0.3;
1022  }
1023  }
1024 
1025  if(isPP05){
1026  if(strcmp(flag,"mb")==0){
1027  //if(RIGHT>0.2) RIGHT=0.19;
1028  if(xval>2.) RIGHT=0.199;
1029  if(LEFT<0.08) LEFT=0.081;
1030  }
1031  if(strcmp(flag,"ht1")==0&&xval>9.){
1032  LEFT=0.081;
1033  RIGHT=0.241;
1034  }
1035  if(strcmp(flag,"ht2")==0&&xval>10.){
1036  LEFT=0.081;
1037  }
1038 
1039 
1040  }
1041 
1042  }
1043 
1044 
1045  x2[g-1]=xval;
1046  ex2[g-1]=0.;
1047 
1048  Float_t dpt=h->GetXaxis()->GetBinWidth(g);
1049  intregion[g]->GetXaxis()->SetRangeUser(LEFT,RIGHT);
1050 
1051 
1052  ptslice_bg[g]->GetXaxis()->SetRangeUser(LEFT,RIGHT);
1053 
1054  y2[g-1]=intregion[g]->Integral();
1055  bgy[g-1]=ptslice_bg[g]->Integral() - intregion[g]->Integral();
1056  if(!isMC) ey2[g-1]=sqrt(y2[g-1]+2.*bgy[g-1]);
1057  else{
1058  //calculate montecarlo error:
1059  ey2[g-1]=0;
1060  for(Int_t i=1;i<=intregion[g]->GetNbinsX();i++){
1061  ey2[g-1]+=pow(intregion[g]->GetBinError(i),2);
1062  }
1063  ey2[g-1]=sqrt(ey2[g-1]);
1064  }
1065  y2[g-1]/=dpt;
1066  ey2[g-1]/=dpt;
1067  //cout<<xval<<"\t"<<y2[g-1]<<" +/- "<<ey2[g-1]<<endl;
1068  if(strcmp(flag,"mb")==0) cout<<xval<<"\t"<<dpt*y2[g-1]/bgy[g-1]<<"\t"<<ey2[g-1]/y2[g-1]<<endl;
1069 
1070  intregion[g]->Draw("histsame");
1071 
1072  //get mean from final gaus fit
1073  m2[g-1]=finalpeak[g]->GetParameter(1);
1074  em2[g-1]=finalpeak[g]->GetParError(1);
1075  s2[g-1]=finalpeak[g]->GetParameter(2);
1076  es2[g-1]=finalpeak[g]->GetParError(2);
1077 
1078  retYield->SetBinContent(g,y2[g-1]);
1079  retYield->SetBinError(g,ey2[g-1]);
1080 
1081  }
1082  cr->cd(0);
1083  cr->Update();
1084  sprintf(canvasname,"c_%s_page_%d",flag,page++);
1085  c_copy=(TCanvas*)cr->Clone();
1086  c_copy->SetName(canvasname);
1087  c_array->Add(c_copy);
1088  ps->NewPage();
1089 
1090 
1091  TCanvas *c4=new TCanvas(TString("c4_")+TString(flag),TString("c4_")+TString(flag),600,900);
1092  c4->Divide(2,3);
1093  c4->cd(1);
1094  TGraphErrors *regry=new TGraphErrors(nx,x2,y2,ex2,ey2);
1095  regry->SetName((TString(flag)+"_yield").Data());
1096  gPad->SetLogy();
1097  regry->SetMarkerStyle(21);
1098  regry->SetMarkerSize(.9);
1099  regry->SetMarkerColor(color);
1100  regry->GetXaxis()->SetRangeUser(pTmin,pTmax);
1101  gry->GetXaxis()->SetRangeUser(pTmin,pTmax);
1102  Float_t maxy=gry->GetHistogram()->GetMaximum();
1103  Float_t miny=gry->GetHistogram()->GetMinimum() + 1.;
1104  if(isMC) miny=gry->GetHistogram()->GetMinimum() + 0.001;
1105  if(isPythia) miny=gry->GetHistogram()->GetMinimum();
1106  regry->SetMaximum(maxy);
1107  regry->SetMinimum(miny);
1108  gry->SetMaximum(maxy);
1109  gry->SetMinimum(miny);
1110  regry->Draw("ap");
1111  c4->cd(2);
1112  gPad->SetLogy();
1113  gry->Draw("ap");
1114  c4->cd(3);
1115  TGraphErrors *regrm=new TGraphErrors(nx,x2,m2,ex2,em2);
1116  regrm->SetName((TString(flag)+"_mean").Data());
1117  regrm->SetMarkerStyle(21);
1118  regrm->SetMarkerSize(.9);
1119  regrm->SetMarkerColor(color);
1120  regrm->GetXaxis()->SetRangeUser(pTmin,pTmax);
1121  grm->GetXaxis()->SetRangeUser(pTmin,pTmax);
1122  regrm->SetMaximum(.2);
1123  regrm->SetMinimum(.1);
1124  grm->SetMaximum(.2);
1125  grm->SetMinimum(.1);
1126  regrm->Draw("ap");
1127  c4->cd(4);
1128  grm->Draw("ap");
1129  c4->cd(5);
1130  TGraphErrors *regrs=new TGraphErrors(nx,x2,s2,ex2,es2);
1131  regrs->SetName((TString(flag)+"_sigma").Data());
1132  regrs->SetMarkerStyle(21);
1133  regrs->SetMarkerSize(.9);
1134  regrs->SetMarkerColor(color);
1135  regrs->GetXaxis()->SetRangeUser(pTmin,pTmax);
1136  grs->GetXaxis()->SetRangeUser(pTmin,pTmax);
1137  regrs->SetMaximum(.05);
1138  regrs->SetMinimum(0.);
1139  grs->SetMaximum(.05);
1140  grs->SetMinimum(0.);
1141  regrs->Draw("ap");
1142  c4->cd(6);
1143  grs->Draw("ap");
1144  c4->cd(0);
1145  c4->Update();
1146 
1147  sprintf(canvasname,"c_%s_page_%d",flag,page++);
1148  c_copy=(TCanvas*)c4->Clone();
1149  c_copy->SetName(canvasname);
1150  c_array->Add(c_copy);
1151 
1152  ps->NewPage();
1153 
1154  if(strcmp(flag,"ht2")==0) ps->Close();
1155 
1156  return retYield;
1157 }
1158 
1159 Float_t Pi0Analysis::getWeight(Float_t pT){
1160  //weight for neutrons (Levy function):
1161  if(pT<0.5) pT=0.5;
1162  Float_t weight=pT/pow((1.+(sqrt(pT*pT+1.) - 1.)/(0.215*12.55)),12.55);
1163  return weight;
1164 }
1165 
1166 void Pi0Analysis::storeCanvases(const Char_t *f_c){
1167  TFile *f_canvas=new TFile(f_c,"RECREATE");
1168  c_array->Write();
1169  f_canvas->Close();
1170 }
MyPoint
Definition: MyPoint.h:7
Generated by   doxygen 1.8.5