StRoot: StEmcPool/StPhotonAnalysis/Pi0Analysis.cxx Source File

00001 #include <stdio.h>
00002 #include <TMultiGraph.h>
00003 #include <TStyle.h>
00004 #include <TMinuit.h>
00005 #include <TTree.h>
00006 #include <TF1.h>
00007 #include <TFile.h>
00008 #include <TCanvas.h>
00009 #include <TVector3.h>
00010 #include <Riostream.h>
00011 #include <TMath.h>
00012 
00013 #include <StEmcPool/StPhotonCommon/MyEvent.h> 
00014 #include <StEmcPool/StPhotonCommon/MyPoint.h> 
00015 #include <StEmcPool/StPhotonCommon/MyMcTrack.h> 
00016 
00017 #include "AnaCuts.h"
00018 #include "EventMixer.h"
00019 
00020 #include "Pi0Analysis.h"
00021 
00022 //temp solution
00023 Bool_t MINBIAS=kTRUE;
00024 Bool_t EXCLUDERANGE=kTRUE;
00025 Bool_t EXCLUDEETA=kFALSE;
00026 Bool_t doNEUTRONS=kFALSE;
00027 Float_t rightPI=0.22;
00028 Float_t leftPI=0.08;//was 0.12
00029 Double_t CONSTANT=0.;
00030 Double_t combFit(Double_t *x,Double_t *par)
00031 {
00032   //old settings:
00033   //leftPI=0.08;
00034   //rightPI=0.24;
00035 
00036   if( (x[0]>leftPI&&x[0]<rightPI) && EXCLUDERANGE )
00037     {
00038       TF1::RejectPoint();
00039       return 0;
00040     }
00041   if( (x[0]>0.5&&x[0]<0.65) && EXCLUDEETA )//thesis result: 0.48
00042     {
00043       TF1::RejectPoint();
00044       return 0;
00045     }
00046  
00047   
00048   Double_t ret=par[0]*x[0]+par[1]*x[0]*x[0]+par[2]*x[0]*x[0]*x[0]+CONSTANT;
00049   
00050   return ret;
00051 }
00052 Double_t sumFit(Double_t *x,Double_t *par)
00053 {
00054   Double_t ret=par[0]*x[0]+par[1]*x[0]*x[0]+par[2]*x[0]*x[0]*x[0];
00055   Double_t w=x[0]-par[4];
00056   Double_t s=par[5];
00057   Double_t ret2=par[3]*TMath::Exp(-0.5*(w/s)*(w/s));
00058   return ret+ret2;
00059 }
00060 
00061 ClassImp(Pi0Analysis)
00062 
00063 Pi0Analysis::Pi0Analysis(const Char_t *psout, const Char_t *psout2, const Char_t *flag)
00064 {
00065   gStyle->SetOptDate(0);
00066   ps=new TPostScript(psout,-111);
00067   ps2=new TPostScript(psout2,-111);
00068 
00069   c_array=new TObjArray(100);
00070 
00071 
00072   cout<<"DEFAULT CONSTRUCTOR FOR PI0ANALYSIS!!!!"<<endl;
00073   cout<<"storing ps in: "<<psout<<endl;
00074 
00075   //init cuts
00076   cuts=new AnaCuts(flag);
00077   //cuts->printCuts();
00078 
00079   //eventmixer
00080   mixer=new EventMixer(flag);
00081   
00082   runPrev=0;
00083   startdatePrev=0;
00084   starttimePrev=0;
00085   numberOfMB=0;
00086   nMBinRun=0;
00087   numberOfHT1=0;
00088   nHT1inRun=0;
00089   numberOfHT2=0;
00090   nHT2inRun=0;
00091   psMB=0;
00092   psHT1=0;
00093   psHT2=0;
00094 
00095   psHT1_eff=0;
00096   psHT1_eff2=0;
00097   nMB_eff=0;
00098 
00099   isMC=kFALSE;
00100   isPythia=kFALSE;
00101   isDAU=kFALSE;
00102   isPP05=kFALSE;
00103   isAUAU200=kFALSE;
00104   isHIJING=kFALSE;
00105   if(strcmp(flag,"dAu")==0){
00106     isDAU=kTRUE;
00107   }
00108   if(strcmp(flag,"pp05")==0){
00109     isPP05=kTRUE;
00110   }
00111   if(strcmp(flag,"auau200")==0){
00112     isAUAU200=kTRUE;
00113   }
00114 
00115   noMINBIAS=kFALSE;
00116 
00117   WEIGHT=1.;
00118   iev_0=0;
00119   iev_1=0;
00120   iev_2=0;
00121 
00122 }
00123 Pi0Analysis::~Pi0Analysis()
00124 {
00125   cout<<endl<<"Pi0Analysis destructed!"<<endl<<endl;
00126 }
00127 Int_t Pi0Analysis::init(const Char_t *output)
00128 {
00129 
00130   mFileOut=new TFile(output,"RECREATE");
00131 
00132   h_bbcVsTpc=new TH2F("h_bbcVsTpc","BBC vs TPC vertex z",481,-240.5,240.5,480,-240.,240.);
00133   h_bbcVsTpcCorr=new TH2F("h_bbcVsTpcCorr","corrected BBC vs TPC vertex z",481,-240.5,240.5,480,-240.,240.);
00134   h_bbcRes=new TH1F("h_bbcres","bbc resolution",100,-2.,2.);
00135 
00136   h_smdeta1=new TH2F("h_smdeta1","n eta strips in cluster vs p_{T}",20,0.,10.,10,0.5,10.5);
00137   h_smdphi1=new TH2F("h_smdphi1","n phi strips in cluster vs p_{T}",20,0.,10.,10,0.5,10.5);
00138   h_smdeta2=new TH2F("h_smdeta2","eta energy ratio in cluster vs p_{T}",20,0.,10.,100,0.0,10.);
00139   h_smdphi2=new TH2F("h_smdphi2","phi energy ratio in cluster vs p_{T}",20,0.,10.,100,0.0,10.);
00140 
00141   h_ratioMB=new TH1F("ratioMB","E neutral over E total",60,-.2,1.);
00142   h_ratioMB->Sumw2();
00143   h_ratioHT1=new TH1F("ratioHT1","E neutral over E total",60,-.2,1.);
00144   h_ratioHT1->Sumw2();
00145   h_ratioHT2=new TH1F("ratioHT2","E neutral over E total",60,-.2,1.);
00146   h_ratioHT2->Sumw2();
00147   h_vzMB=new TH1F("h_vzMB","z-vertex MB",480,-120.,120.);
00148   h_vzMB->Sumw2();
00149   h_vzHT1=new TH1F("h_vzHT1","z-vertex HT1",480,-120.,120.);
00150   h_vzHT1->Sumw2();
00151   h_vzHT2=new TH1F("h_vzHT2","z-vertex HT2",480,-120.,120.);
00152   h_vzHT2->Sumw2();
00153   h_trigidHT1=new TH1F("h_trigidHT1","id of trigger HT1",4800,0.5,4800.5);
00154   h_trigidHT1->Sumw2();
00155   h_trigidHT2=new TH1F("h_trigidHT2","id of trigger HT2",4800,0.5,4800.5);
00156   h_trigidHT2->Sumw2();
00157   h_trigidHT1aftercut=new TH1F("h_trigidHT1aftercut","id of trigger HT1 after cut",4800,0.5,4800.5);
00158   h_trigidHT1aftercut->Sumw2();
00159   h_trigidHT2aftercut=new TH1F("h_trigidHT2aftercut","id of trigger HT2 after cut",4800,0.5,4800.5);
00160   h_trigidHT2aftercut->Sumw2();
00161   h_events=new TH1F("h_events","trigger counts: mb=+1,ht1=+2,ht2=+4,sn=+8",21,-0.5,20.5);
00162   h_events->Sumw2();
00163 
00164   h_HT1adc_id=new TH2F("h_HT1adc_id","ht-1 adc vs id",2400,.5,2400.5,800,200.5,1000.5);
00165   h_HT2adc_id=new TH2F("h_HT2adc_id","ht-2 adc vs id",2400,.5,2400.5,800,200.5,1000.5);
00166 
00167   h_etaphi=new TH2F("h_etaphi","eta/phi of neutral points",300,0.0,1.0,1800,-3.14,3.14);
00168   h_etaphi->Sumw2();
00169   h_rapphi=new TH2F("h_rapphi","rap/phi of neutral points",300,-0.2,1.2,1800,-3.14,3.14);
00170   h_rapphi->Sumw2();
00171   h_dist=new TH1F("h_dist","distance of point to track",1000,0.0,200.0);
00172   h_dist->Sumw2();
00173 
00174   h_dist2DMB=new TH2F("h_dist2DMB","dist. vs pt MB",cuts->nPtBinsMB,cuts->ptBinsMB.GetArray(),1000,0.0,1000.0);
00175   h_dist2DHT1=new TH2F("h_dist2DHT1","dist. vs pt HT1",cuts->nPtBinsHT1,cuts->ptBinsHT1.GetArray(),1000,0.0,1000.0);
00176   h_dist2DHT2=new TH2F("h_dist2DHT2","dist. vs pt HT2",cuts->nPtBinsHT2,cuts->ptBinsHT2.GetArray(),1000,0.0,1000.0);
00177 
00178   h_dist2DMBpions=new TH2F("h_dist2DMBpions","dist. vs pt MB",cuts->nPtBinsMB,cuts->ptBinsMB.GetArray(),1000,0.0,1000.0);
00179   h_dist2DHT1pions=new TH2F("h_dist2DHT1pions","dist. vs pt HT1",cuts->nPtBinsHT1,cuts->ptBinsHT1.GetArray(),1000,0.0,1000.0);
00180   h_dist2DHT2pions=new TH2F("h_dist2DHT2pions","dist. vs pt HT2",cuts->nPtBinsHT2,cuts->ptBinsHT2.GetArray(),1000,0.0,1000.0);
00181   
00182   h_asymmMB=new TH2F("h_asymmMB","asymmetry of inv mass pairs MB",cuts->nPtBinsMB,cuts->ptBinsMB.GetArray(),40,0.0,1.0);
00183   h_asymmMB->Sumw2();
00184   h_asymmHT1=new TH2F("h_asymmHT1","asymmetry of inv mass pairs HT1",cuts->nPtBinsHT1,cuts->ptBinsHT1.GetArray(),40,0.0,1.0);
00185   h_asymmHT1->Sumw2();
00186   h_asymmHT2=new TH2F("h_asymmHT2","asymmetry of inv mass pairs HT2",cuts->nPtBinsHT2,cuts->ptBinsHT2.GetArray(),40,0.0,1.0);
00187   h_asymmHT2->Sumw2();
00188   h_asymmMBbg=new TH2F("h_asymmMBbg","asymmetry of inv mass pairs MB (bg)",cuts->nPtBinsMB,cuts->ptBinsMB.GetArray(),40,0.0,1.0);
00189   h_asymmMBbg->Sumw2();
00190   h_asymmHT1bg=new TH2F("h_asymmHT1bg","asymmetry of inv mass pairs HT1 (bg)",cuts->nPtBinsHT1,cuts->ptBinsHT1.GetArray(),40,0.0,1.0);
00191   h_asymmHT1bg->Sumw2();
00192   h_asymmHT2bg=new TH2F("h_asymmHT2bg","asymmetry of inv mass pairs HT2 (bg)",cuts->nPtBinsHT2,cuts->ptBinsHT2.GetArray(),40,0.0,1.0);
00193   h_asymmHT2bg->Sumw2();
00194 
00195   h_minvMB=new TH2F("h_minvMB","invariant mass vs pT MB",cuts->nPtBinsMB,cuts->ptBinsMB.GetArray(),cuts->nMinvBinsMB,cuts->mInvLowMB,cuts->mInvHighMB);
00196   h_minvMB->Sumw2();
00197   h_minvHT1=new TH2F("h_minvHT1","invariant mass vs pT HT1",cuts->nPtBinsHT1,cuts->ptBinsHT1.GetArray(),cuts->nMinvBinsHT1,cuts->mInvLowHT1,cuts->mInvHighHT1);
00198   h_minvHT1->Sumw2();
00199   h_minvHT2=new TH2F("h_minvHT2","invariant mass vs pT HT2",cuts->nPtBinsHT2,cuts->ptBinsHT2.GetArray(),cuts->nMinvBinsHT2,cuts->mInvLowHT2,cuts->mInvHighHT2);
00200   h_minvHT2->Sumw2();
00201 
00202   h_pionsVsEtaMB=new TH1F("h_pionsVsEta","dN_{#pi^{0}}/dy data",40,-0.5,1.5);
00203 
00204   h_gammaMB=new TH1F("h_gammaMB","photon yield vs pT MB",cuts->nPtBinsMB,cuts->ptBinsMB.GetArray());
00205   h_gammaMB->Sumw2();
00206   h_gammaHT1=new TH1F("h_gammaHT1","photon yield vs pT HT1",cuts->nPtBinsHT1,cuts->ptBinsHT1.GetArray());
00207   h_gammaHT1->Sumw2();
00208   h_gammaHT2=new TH1F("h_gammaHT2","photon yield vs pT HT2",cuts->nPtBinsHT2,cuts->ptBinsHT2.GetArray());
00209   h_gammaHT2->Sumw2();
00210 
00211   h_pythiaPartonPt=new TH1F("h_pythiaPartonPt","pT of pythia process",200,0.,50);  
00212   h_pythiaPartonPt->Sumw2();
00213   h_pythiaPions=new TH1F("h_pythiaPions","pythia pions vs pT",200,0.,50.);
00214   h_pythiaPions->Sumw2();
00215   h_pythiaPionsMB=new TH1F("h_pythiaPionsMB","pythia pions vs pT MB",cuts->nPtBinsMB,cuts->ptBinsMB.GetArray());
00216   h_pythiaPionsMB->Sumw2();
00217   h_pythiaPionsHT1=new TH1F("h_pythiaPionsHT1","pythia pions vs pT HT1",cuts->nPtBinsHT1,cuts->ptBinsHT1.GetArray());
00218   h_pythiaPionsHT1->Sumw2();
00219   h_pythiaPionsHT2=new TH1F("h_pythiaPionsHT2","pythia pions vs pT HT2",cuts->nPtBinsHT2,cuts->ptBinsHT2.GetArray());
00220   h_pythiaPionsHT2->Sumw2();
00221 
00222   isHot=(TArrayI)cuts->isHot;
00223 
00224   h_adcHT1=new TH1F("h_adcHT1","highest adc en>0. for HT1 after cuts",1000,0.,1000.);
00225   h_adcHT2=new TH1F("h_adcHT2","highest adc en>0. for HT2 after cuts",1000,0.,1000.);
00226 
00227   h_mcneutronsMB=new TH1F("h_mcneutronsMB","gen. neutrons vs pT MB",cuts->nPtBinsEffMB,cuts->ptBinsEffMB.GetArray());
00228   h_mcneutronsHT1=new TH1F("h_mcneutronsHT1","gen. neutrons vs pT HT1",cuts->nPtBinsEff,cuts->ptBinsEff.GetArray());
00229   h_mcneutronsHT2=new TH1F("h_mcneutronsHT2","gen. neutrons vs pT HT2",cuts->nPtBinsEff,cuts->ptBinsEff.GetArray());
00230   h_mcneutronsWeightMB=new TH1F("h_mcneutronsWeightMB","gen. neutrons vs pT (weighted) MB",cuts->nPtBinsEff,cuts->ptBinsEff.GetArray());
00231   h_mcneutronsWeightHT1=new TH1F("h_mcneutronsWeightHT1","gen. neutrons vs pT (weighted) HT1",cuts->nPtBinsEff,cuts->ptBinsEff.GetArray());
00232   h_mcneutronsWeightHT2=new TH1F("h_mcneutronsWeightHT2","gen. neutrons vs pT (weighted) HT2",cuts->nPtBinsEff,cuts->ptBinsEff.GetArray());
00233   h_neutronsMB=new TH1F("h_neutronsMB","reco. neutrons vs pT MB",cuts->nPtBinsEff,cuts->ptBinsEff.GetArray());
00234   h_neutronsHT1=new TH1F("h_neutronsHT1","reco. neutrons vs pT HT1",cuts->nPtBinsEff,cuts->ptBinsEff.GetArray());
00235   h_neutronsHT2=new TH1F("h_neutronsHT2","reco. neutrons vs pT HT2",cuts->nPtBinsEff,cuts->ptBinsEff.GetArray());
00236 
00237   h_EvsE=new TH2F("h_EvsE","E neutral vs E in TPC",80,0.,40,80,0.,40.);
00238 
00239   h_nstripsETA=new TH2F("h_nstripsETA","n strips vs pt bsmde",15,0.,15.,6,0.,6.);
00240   h_nstripsPHI=new TH2F("h_nstripsPHI","n strips vs pt bsmdp",15,0.,15.,6,0.,6.);
00241 
00242   h_clusterWidth=new TH2F("h_clusterWidth","width BSMD eta+phi",20,0.,20.,100,0.,0.05);
00243   h_clusterWidth->Sumw2();
00244   h_energyRatio=new TH2F("h_energyRatio","energy ratio BSMD/BTOW",20,0.,20.,160,0.,8.);
00245   h_energyRatio->Sumw2();
00246 
00247   if(isPP05){
00248     fout_mb.open("timestamps/pp_timestamps_mb.list");
00249     fout_ht1.open("timestamps/pp_timestamps_ht1.list");
00250     fout_ht2.open("timestamps/pp_timestamps_ht2.list");
00251   }
00252   if(isDAU){
00253     fout_mb.open("timestamps/dau_timestamps_mb.list");
00254     fout_ht1.open("timestamps/dau_timestamps_ht1.list");
00255     fout_ht2.open("timestamps/dau_timestamps_ht2.list");
00256   }
00257 
00258   return 1; 
00259 }
00260 Int_t Pi0Analysis::make(Int_t evmax, const Char_t *input)
00261 {
00262   mFile=new TFile(input,"OPEN");
00263   myEventTree=(TTree*)mFile->Get("mEventTree");
00264   ev=new MyEvent();
00265   myEventTree->SetBranchAddress("branch",&ev);
00266   
00267   Int_t i=0;
00268   while(myEventTree->GetEntry(i)){
00269     if(evmax&&i>evmax){
00270       cout<<"reached evmax,"<<endl;
00271       cout<<"abort loop!"<<endl;
00272       break;      
00273     }
00274     if(i%500000==0) cout<<"processing "<<input<<" evt: "<<i<<endl;
00275 
00276     if(i==0){
00277       startdatePrev=ev->date();
00278       starttimePrev=ev->time();
00279     }
00280 
00281     //assign trigger for MC event:
00282     Int_t mcTr=1;//mb
00283     if(isMC&&!isPythia){
00284       if(ev->highTowerAdc()>cuts->ht1AdcCUT) mcTr+=2;
00285       if(ev->highTowerAdc()>cuts->ht2AdcCUT) mcTr+=4;
00286       ev->setTrigger(mcTr);
00287     }
00288 
00289     TVector3 vPos=ev->vertex();
00290 
00291     float bbc_vertz=(ev->bbcVertexZ()-10.77)/2.82;
00292     //now bbc_vert is west-east, then from correlation:
00293     bbc_vertz+=4.04;
00294     bbc_vertz/=0.354;
00295     ev->setBbcVertexZ(bbc_vertz);
00296 
00297     if(ev->trigger()&1 && TMath::Abs(vPos.Z())>0.00000001){
00298       h_bbcVsTpc->Fill(vPos.Z(),ev->bbcVertexZ());
00299     }
00300    
00301     if(ev->trigger()&1) {h_vzMB->Fill(vPos.Z());}
00302     if(ev->trigger()&2) {h_vzHT1->Fill(vPos.Z());h_trigidHT1->Fill(ev->highTowerId());}
00303     if(ev->trigger()&4) {h_vzHT2->Fill(vPos.Z());h_trigidHT2->Fill(ev->highTowerId());}
00304 
00305     if(ev->trigger()&4 && TMath::Abs(vPos.Z())<60.){h_EvsE->Fill(ev->momentumInTpc(),ev->energyInBarrel());}
00306     
00307     //set tpc momentum for simu:
00308     if(isMC) ev->setMomentumInTpc(99999.);
00309 
00310     //exclude depleted minbias
00311     if(noMINBIAS && ev->trigger()&1){
00312       int trigg=ev->trigger();
00313       ev->setTrigger(trigg-1);
00314     }
00315 
00316     //event cuts + hot tower:
00317     if(isDAU && !cuts->isEventOK(ev,"dAu")) {i++;continue;}
00318     if(isPP05 && !cuts->isEventOK(ev,"pp05")) {i++;continue;}    
00319 
00320 
00321 
00322     //fill mixer
00323     //if(isDAU || isPP05) mixer->addEvent(ev);
00324 
00325     //check resolution bbc:
00326     if(isPP05 && TMath::Abs(vPos.Z())>0.){
00327       float bbcres=(ev->bbcVertexZ()-vPos.Z())/vPos.Z();
00328       h_bbcRes->Fill(bbcres);
00329     }
00330 
00331     Float_t ratioE=TMath::Abs(ev->energyInBarrel()+ev->momentumInTpc())>0. ? ev->energyInBarrel()/(ev->energyInBarrel()+ev->momentumInTpc()) : -.1;
00332 
00333     if(i>0 && ev->runId()!=runPrev){
00334       psHT1_eff+=psHT1*nHT1inRun;
00335       nMB_eff+=psMB*nMBinRun;
00336 
00337       fout_mb<<runPrev<<" "<<startdatePrev<<" "<<starttimePrev<<" "<<nMBinRun<<endl;
00338       fout_ht1<<runPrev<<" "<<startdatePrev<<" "<<starttimePrev<<" "<<nHT1inRun<<endl;
00339       fout_ht2<<runPrev<<" "<<startdatePrev<<" "<<starttimePrev<<" "<<nHT2inRun<<endl;
00340   
00341       nMBinRun=0;
00342       nHT1inRun=0;
00343       nHT2inRun=0;
00344       startdatePrev=ev->date();
00345       starttimePrev=ev->time();
00346     }
00347     if(ev->trigger()&1) {nMBinRun++;}
00348     if(ev->trigger()&2) {nHT1inRun++;}
00349     if(ev->trigger()&4) {nHT2inRun++;}
00350     psMB=(ev->prescale(0) > ev->prescale(1)) ? (Float_t)ev->prescale(0) : (Float_t)ev->prescale(1);
00351     psHT1=(Float_t)ev->prescale(2);
00352     psHT2=(Float_t)ev->prescale(3);
00353     
00354     //get neutrons:
00355     for(Int_t j=0;j<ev->numberOfMcPions();j++){
00356       MyMcTrack *pyt=(MyMcTrack*)ev->getMcPionArray()->At(j);
00357       //check if neutron in rapidity bite:
00358       //change this to rapidity...
00359       if(pyt->momentum().PseudoRapidity()<cuts->rapidityMinCUT||pyt->momentum().PseudoRapidity()>cuts->rapidityMaxCUT) continue;
00360       h_mcneutronsMB->Fill(pyt->momentum().Pt());
00361       h_mcneutronsHT1->Fill(pyt->momentum().Pt());
00362       h_mcneutronsHT2->Fill(pyt->momentum().Pt());
00363       Float_t weight=getWeight(pyt->momentum().Pt());
00364       h_mcneutronsWeightMB->Fill(pyt->momentum().Pt(),weight);
00365       h_mcneutronsWeightHT1->Fill(pyt->momentum().Pt(),weight);
00366       h_mcneutronsWeightHT2->Fill(pyt->momentum().Pt(),weight);
00367     }
00368 
00369     if(ev->trigger()&1){
00370       h_ratioMB->Fill(ratioE);
00371       numberOfMB++;
00372     }
00373     if(ev->trigger()&2){
00374       h_ratioHT1->Fill(ratioE);
00375       numberOfHT1++;
00376       h_trigidHT1aftercut->Fill(ev->highTowerId());
00377 
00378       h_adcHT1->Fill(ev->highTowerAdc());
00379       h_HT1adc_id->Fill(ev->highTowerId(),ev->highTowerAdc());
00380       }
00381     if(ev->trigger()&4){
00382       h_ratioHT2->Fill(ratioE);
00383       numberOfHT2++;
00384       h_trigidHT2aftercut->Fill(ev->highTowerId());
00385 
00386       h_adcHT2->Fill(ev->highTowerAdc());
00387       h_HT2adc_id->Fill(ev->highTowerId(),ev->highTowerAdc());
00388     }
00389 
00390     h_events->Fill(ev->trigger());
00391 
00392     if(isPythia){
00393       Float_t pypT=ev->partonPt();
00394       if(pypT<15.){
00395         WEIGHT=1.;
00396         iev_0++;
00397       }
00398       else if(pypT<25.){
00399         WEIGHT=(216000/208000)*(0.0003895/0.002228);
00400         iev_1++;
00401       }
00402       else if(pypT>25.){
00403         iev_2++;
00404         WEIGHT=(216000/208000)*(1.016e-005/0.002228);
00405       }
00406       h_pythiaPartonPt->Fill(ev->partonPt(),WEIGHT);
00407       for(Int_t it=0;it<ev->getMcPionArray()->GetEntries();it++){
00408         MyMcTrack *pyt=(MyMcTrack*)ev->getMcPionArray()->At(it);
00409         if(pyt->momentum().PseudoRapidity()<cuts->rapPionMinCUT||pyt->momentum().PseudoRapidity()>cuts->rapPionMaxCUT) continue;
00410         h_pythiaPions->Fill(pyt->momentum().Pt(),WEIGHT);
00411         if(ev->trigger()&1){
00412           h_pythiaPionsMB->Fill(pyt->momentum().Pt(),WEIGHT);
00413         }
00414         if(ev->trigger()&2){
00415           h_pythiaPionsHT1->Fill(pyt->momentum().Pt(),WEIGHT);
00416         }
00417         if(ev->trigger()&4){
00418           h_pythiaPionsHT2->Fill(pyt->momentum().Pt(),WEIGHT);
00419         }
00420       }
00421     }
00422 
00423     runPrev=ev->runId();
00424     //loop on points
00425     MyPoint *p;
00426     MyPoint *pp;
00427     TClonesArray *clA=(TClonesArray*)ev->getPointArray();
00428     for(Int_t j=0;j<clA->GetEntries();j++){
00429       p=(MyPoint*)clA->At(j);
00430       TVector3 pPos=p->position();
00431       TVector3 pMom=pPos - vPos;
00432       pMom.SetMag(p->energy());
00433 
00434       if(!cuts->isPointOK(p,vPos)) continue;
00435       
00436       h_dist->Fill(p->distanceToTrack());
00437       if(ev->trigger()&1) h_dist2DMB->Fill(pMom.Pt(),p->distanceToTrack());
00438       if(ev->trigger()&2) h_dist2DHT1->Fill(pMom.Pt(),p->distanceToTrack());
00439       if(ev->trigger()&4) h_dist2DHT2->Fill(pMom.Pt(),p->distanceToTrack());
00440 
00441       if(ev->trigger()&1 && pMom.Pt()>1.){
00442         h_etaphi->Fill(pPos.PseudoRapidity(),pPos.Phi());
00443         h_rapphi->Fill(pMom.PseudoRapidity(),pMom.Phi());
00444       }
00445 
00446       if(pMom.PseudoRapidity()>cuts->rapidityMinCUT && pMom.PseudoRapidity()<cuts->rapidityMaxCUT){
00447         //first neutrons:
00448         if(doNEUTRONS){
00449           Float_t D=1000000.;
00450           Float_t pt4weight=10.;
00451           for(Int_t j=0;j<ev->numberOfMcPions();j++){
00452             MyMcTrack *pyt=(MyMcTrack*)ev->getMcPionArray()->At(j);
00453             TVector3 mcPos=pyt->position();
00454             Float_t n=mcPos.PseudoRapidity();
00455             Float_t p=mcPos.Phi();
00456             Float_t nn=pPos.PseudoRapidity();
00457             Float_t pp=pPos.Phi();
00458             Float_t d=sqrt(pow(p-pp,2)+pow(n-nn,2));
00459             if(d<D){
00460               D=d;
00461               pt4weight=pyt->momentum().Pt();
00462             }
00463           }
00464           if(D<0.025){
00465             if(ev->trigger()&1) h_neutronsMB->Fill(pMom.Pt(),getWeight(pt4weight));
00466             if(ev->trigger()&2) h_neutronsHT1->Fill(pMom.Pt(),getWeight(pt4weight));
00467             if(ev->trigger()&4) h_neutronsHT2->Fill(pMom.Pt(),getWeight(pt4weight));
00468           }
00469         }
00470         
00471         //fill photons:
00472         if(p->distanceToTrack()>cuts->photonCUT){
00473           if(ev->trigger()&1) h_gammaMB->Fill(pMom.Pt(),WEIGHT);
00474           if(ev->trigger()&2) h_gammaHT1->Fill(pMom.Pt(),WEIGHT);
00475           if(ev->trigger()&4) h_gammaHT2->Fill(pMom.Pt(),WEIGHT);
00476           
00477           if(ev->trigger()&1){
00478             h_clusterWidth->Fill(pMom.Pt(),sqrt(p->widthEta()*p->widthEta()+p->widthPhi()*p->widthPhi()));
00479             h_energyRatio->Fill(pMom.Pt(),(p->energyEta()+p->energyPhi())/p->energy());
00480             h_nstripsETA->Fill(pMom.Pt(),p->nHitsEta(),WEIGHT);
00481             h_nstripsPHI->Fill(pMom.Pt(),p->nHitsPhi(),WEIGHT);
00482             //super isolated:
00483             if(p->distanceToTrack()>30.){
00484               h_smdeta1->Fill(pMom.Pt(),p->nHitsEta());
00485               h_smdphi1->Fill(pMom.Pt(),p->nHitsPhi());
00486               h_smdeta2->Fill(pMom.Pt(),p->energyEta()/p->energy());
00487               h_smdphi2->Fill(pMom.Pt(),p->energyPhi()/p->energy());
00488             }
00489           }
00490         }
00491         
00492       }
00493       
00494       for(Int_t jj=0;jj<clA->GetEntries();jj++){
00495         if(jj<=j) continue;
00496         pp=(MyPoint*)clA->At(jj);
00497         TVector3 ppPos=pp->position();
00498         TVector3 ppMom=ppPos - vPos;
00499         ppMom.SetMag(pp->energy());
00500         if(!cuts->isPointOK(pp,vPos)) continue;
00501 
00502         if(p->distanceToTrack()<pp->distanceToTrack()){
00503           if(ev->trigger()&1) h_dist2DMBpions->Fill(pMom.Pt(),p->distanceToTrack());
00504           if(ev->trigger()&2) h_dist2DHT1pions->Fill(pMom.Pt(),p->distanceToTrack());
00505           if(ev->trigger()&4) h_dist2DHT2pions->Fill(pMom.Pt(),p->distanceToTrack());
00506         }
00507         else{
00508           if(ev->trigger()&1) h_dist2DMBpions->Fill(pMom.Pt(),pp->distanceToTrack());
00509           if(ev->trigger()&2) h_dist2DHT1pions->Fill(pMom.Pt(),pp->distanceToTrack());
00510           if(ev->trigger()&4) h_dist2DHT2pions->Fill(pMom.Pt(),pp->distanceToTrack());
00511         }
00512            
00513         //two neutrals
00514         TVector3 pi0Mom=pMom+ppMom;
00515         Float_t angle=pMom.Angle(ppMom);
00516         Float_t Epion=p->energy()+pp->energy();
00517         Float_t asymm=TMath::Abs(p->energy()-pp->energy())/Epion;
00518         Float_t minv=sqrt(2.*p->energy()*pp->energy()*(1. - cos(angle)));
00519         Float_t pTpion=pi0Mom.Pt();
00520         
00521         //select rapidity bite / asymm cut for pions:
00522         if(pi0Mom.PseudoRapidity()<cuts->rapPionMinCUT||pi0Mom.PseudoRapidity()>cuts->rapPionMaxCUT) continue;
00523                 
00524         if(ev->trigger()&1){
00525           if(asymm<=cuts->asymmetryCUT) h_minvMB->Fill(pTpion,minv,WEIGHT);
00526           if(pTpion>1.){
00527             if(minv>0.08 && minv<0.20){
00528             h_asymmMB->Fill(pTpion,asymm);
00529             if(asymm<=cuts->asymmetryCUT) h_pionsVsEtaMB->Fill(pi0Mom.PseudoRapidity());
00530             }
00531             else h_asymmMBbg->Fill(pTpion,asymm);
00532           }
00533         }
00534         if(ev->trigger()&2){
00535           if(asymm<=cuts->asymmetryCUT) h_minvHT1->Fill(pTpion,minv,WEIGHT);
00536           if(pTpion>4.){
00537             if(minv>0.10 && minv<0.20){
00538               h_asymmHT1->Fill(pTpion,asymm);
00539             }
00540             else h_asymmHT1bg->Fill(pTpion,asymm);
00541           }
00542         }
00543         if(ev->trigger()&4){
00544           if(asymm<=cuts->asymmetryCUT) h_minvHT2->Fill(pTpion,minv,WEIGHT);
00545           if(pTpion>6.){
00546             if(minv>0.10 && minv<0.2){
00547               h_asymmHT2->Fill(pTpion,asymm);
00548             }
00549             else h_asymmHT2bg->Fill(pTpion,asymm);
00550           }
00551         }
00552       
00553       }
00554 
00555 
00556     }
00557     
00558     
00559     i++;
00560   }
00561   
00562   return 1;
00563 }
00564 
00565 void Pi0Analysis::printPrescales()
00566 {
00567   cout<<"--------------------------------"<<endl;
00568   cout<<"      stats processed:          "<<endl;
00569   cout<<"                                "<<endl;
00570   cout<<" number of mb:  "<<numberOfMB<<endl;
00571   cout<<" number of ht1: "<<numberOfHT1<<endl;
00572   cout<<" number of ht2: "<<numberOfHT2<<endl;
00573   cout<<"--------------------------------"<<endl;
00574   cout<<" prescale corrections:          "<<endl;
00575   cout<<"                                "<<endl;
00576   cout<<" psHT1_eff "<<psHT1_eff/(1.*numberOfHT1)<<endl;
00577   cout<<endl;
00578   cout<<" N_mb*ps_mb: "<<nMB_eff<<endl;
00579   cout<<"--------------------------------"<<endl;
00580 }
00581 Int_t Pi0Analysis::finish()
00582 {
00583   h_minvMB_mixed=new TH2F(*(TH2F*)mixer->getMinvMB());
00584   mFileOut->Add(h_minvMB_mixed);
00585   mFileOut->Write();
00586 
00587   return 1;
00588 }
00589 
00590 void Pi0Analysis::getYield()
00591 {
00592   h_yieldMB=new TH1F(*getYield(NULL,"mb"));
00593   h_yieldHT1=new TH1F(*getYield(NULL,"ht1"));
00594   h_yieldHT2=new TH1F(*getYield(NULL,"ht2"));
00595 }
00596 
00597 //new routine:
00598 TH1F *Pi0Analysis::getYield(TH2F *h, const Char_t *flag)
00599 {
00600   ps->On();
00601   if(!h){
00602     cout<<"zero pointer"<<endl;
00603     if(strcmp(flag,"mb")==0) h=new TH2F(*h_minvMB);
00604     if(strcmp(flag,"ht1")==0) h=new TH2F(*h_minvHT1);
00605     if(strcmp(flag,"ht2")==0) h=new TH2F(*h_minvHT2);
00606   }
00607   TString retname(h->GetName());
00608   retname.Append("_yield");
00609   TString rettitle(h->GetTitle());
00610   rettitle.Append("_yield");
00611   TH1F *retYield=new TH1F(*((TH1F*)h->ProjectionX()));
00612   retYield->Reset();
00613   retYield->Sumw2();
00614   retYield->SetNameTitle(retname.Data(),rettitle.Data());
00615 
00616   Char_t canvasname[100];
00617   int page=0;
00618   TCanvas *c_copy;
00619 
00620   Float_t NSIGMAHI=3.;
00621   Float_t NSIGMALO=3.;
00622 
00623   const int nx=h->GetNbinsX();
00624   Float_t *x=new Float_t[nx];
00625   Float_t *ex=new Float_t[nx];
00626   Float_t *y=new Float_t[nx];
00627   Float_t *ey=new Float_t[nx];
00628   Float_t *s=new Float_t[nx];
00629   Float_t *es=new Float_t[nx];
00630   Float_t *m=new Float_t[nx];
00631   Float_t *em=new Float_t[nx];
00632   
00633   cout<<"------------- "<<flag<<" --------------"<<endl;
00634 
00635   Int_t color=1;
00636   Int_t cbgcolor=1;
00637   Int_t peakcolor=1;
00638   Int_t fillcolor=16;
00639   MINBIAS=kTRUE;
00640   if(strcmp(flag,"ht1")==0) {color=4; cbgcolor=2; peakcolor=1; fillcolor=38; MINBIAS=kFALSE;}
00641   if(strcmp(flag,"ht2")==0) {color=2; cbgcolor=4; peakcolor=1; fillcolor=46; MINBIAS=kFALSE;}
00642 
00643   TCanvas *c=new TCanvas(flag,flag,600,600);
00644   c->Divide(3,3);
00645 
00646   TF1 *combbg[50];
00647   Char_t combname[100];
00648 
00649   TF1 *pi0peak[50];
00650   Char_t peakname[100];
00651   TF1 *sumfit[50];
00652   Char_t sumname[100];
00653 
00654   TH1D *ptslice[50];
00655   Char_t name[100];
00656 
00657   TH1D *ptslice_bg[50];
00658 
00659   TH1D *ptslice_bgsub[50];
00660   Char_t subname[100];
00661   Char_t title[100];
00662   
00663   for(Int_t i=1;i<=nx;i++){
00664     if((i-1)%9==0){
00665       if(i>1){
00666         c->Update();
00667         sprintf(canvasname,"c_%s_page_%d",flag,page++);
00668         c_copy=(TCanvas*)c->Clone();
00669         c_copy->SetName(canvasname);
00670         c_array->Add(c_copy);
00671       }
00672       ps->NewPage();
00673     }
00674     c->cd((i-1)%9 + 1);
00675       
00676     sprintf(combname,"comb_%d_%s",i,flag);
00677    
00678     sprintf(peakname,"gaus_%d_%s",i,flag);
00679     sprintf(sumname,"sum_%d_%s",i,flag);
00680     sprintf(name,"ptslice_%d_%s",i,flag);
00681     sprintf(subname,"subslice_%d_%s",i,flag);
00682     sprintf(title,"%.2f < p_{T} < %.2f",h->GetXaxis()->GetBinLowEdge(i),h->GetXaxis()->GetBinUpEdge(i));
00683     
00684     //get bin
00685     ptslice[i]=h->ProjectionY(name,i,i,"e");
00686 
00687     //for next iter.
00688     ptslice_bg[i]=new TH1D(*ptslice[i]);
00689 
00690 
00691     //set fitfunction
00692     combbg[i]=new TF1(combname,combFit,0.0,1.0,3);
00693 
00694 
00695     sumfit[i]=new TF1(sumname,sumFit,0.0,1.0,6);
00696    
00697 
00698     float firstbin=0;
00699     for(int ib=1;ib<=ptslice[i]->GetNbinsX();ib++){
00700       if(ptslice[i]->GetBinContent(ib)>0.){
00701         firstbin=(float)ptslice[i]->GetXaxis()->GetBinCenter(ib);
00702         break;
00703       }
00704     }
00705 
00706     combbg[i]->SetParameters(1.0,-1.0,1.0);
00707     
00708     EXCLUDERANGE=kTRUE;//leave out pi0 range, interval:
00709     EXCLUDEETA=kTRUE;
00710     if(strcmp(flag,"mb")==0) rightPI=isDAU ? 0.2 : 0.2;//dau was 0.26 : 0.22
00711     if(strcmp(flag,"ht1")==0) rightPI=isDAU ? 0.281 : 0.281;
00712     if(strcmp(flag,"ht2")==0) rightPI=isDAU ? 0.26 : 0.26;
00713     
00714     if(isDAU && strcmp(flag,"ht1")==0){
00715       if(h->GetXaxis()->GetBinCenter(i)>9.){
00716         combbg[i]->FixParameter(2,0.);//was 9.
00717         //rightPI=0.3;
00718       }
00719     }
00720     if(isDAU && strcmp(flag,"ht2")==0){
00721       if(h->GetXaxis()->GetBinCenter(i)>11.){//was 11.
00722         combbg[i]->FixParameter(2,0.);//was 9.
00723         //rightPI=0.3;
00724       }
00725     }
00726 
00727     if(isPP05 && h->GetXaxis()->GetBinCenter(i)>9.) combbg[i]->FixParameter(2,0.);
00728 
00729     combbg[i]->SetRange(0.0,0.8);
00730     if(isPP05){
00731       combbg[i]->SetRange(0.28,0.8);//was 0.0-0.8 (24/01/08)
00732       if(h->GetXaxis()->GetBinCenter(i)<7. && strcmp(flag,"ht2")==0){
00733       combbg[i]->SetRange(0.3,0.8);
00734       }
00735     }
00736     if(isDAU) combbg[i]->SetRange(0.28,0.8);    
00737 
00738     sumfit[i]->SetRange(0.0,0.8);//was 0.0-0.8, this to show background
00739     
00740     
00741     //MARKER    
00742     if(isPP05 && strcmp(flag,"mb")==0){
00743       combbg[i]->SetRange(0.0,0.4);//0.8
00744       if(h->GetXaxis()->GetBinCenter(i)<2.){
00745         combbg[i]->SetRange(0.05,0.4);
00746         leftPI=0.081;
00747         rightPI=0.19;
00748       }
00749       if(h->GetXaxis()->GetBinCenter(i)>2.){
00750         combbg[i]->FixParameter(2,0.);
00751         combbg[i]->SetRange(0.05,0.4);
00752         leftPI=0.10;
00753         rightPI=0.22;
00754       }
00755       if(h->GetXaxis()->GetBinCenter(i)>3.){
00756         leftPI=0.06;
00757         rightPI=0.22;
00758       }
00759       sumfit[i]->SetRange(0.0,0.4);
00760     }
00761     if(isDAU && strcmp(flag,"mb")==0){
00762       combbg[i]->SetRange(0.0,0.4);//0.8
00763       if(h->GetXaxis()->GetBinCenter(i)<2.){
00764         combbg[i]->SetRange(0.0,0.4);
00765         leftPI=0.0;
00766         rightPI=0.2;
00767       }
00768       if(h->GetXaxis()->GetBinCenter(i)>2.){
00769         leftPI=0.0;
00770         rightPI=0.22;
00771       }
00772       if(h->GetXaxis()->GetBinCenter(i)>4.){
00773         combbg[i]->FixParameter(2,0.);
00774         leftPI=0.0;
00775         rightPI=0.22;
00776       }
00777       sumfit[i]->SetRange(0.0,0.4);
00778     }
00779 
00780 
00781     if(!isMC||isPythia) ptslice[i]->Fit(combname,"QR0");
00782     EXCLUDERANGE=kFALSE;
00783     EXCLUDEETA=kFALSE;
00784     combbg[i]->SetRange(0.0,0.8);
00785     if(strcmp(flag,"mb")==0) combbg[i]->SetRange(0.0,0.4);
00786     //cout<<gMinuit->fCstatu.Data()<<endl;
00787     
00788     ptslice[i]->SetMarkerStyle(4);
00789     ptslice[i]->SetMarkerSize(.7);
00790     ptslice[i]->SetMarkerColor(color);
00791     ptslice[i]->SetLineColor(color);
00792     ptslice[i]->Draw();
00793     combbg[i]->SetLineStyle(2);
00794     combbg[i]->SetLineWidth(1);
00795     combbg[i]->SetLineColor(cbgcolor);
00796     sumfit[i]->SetLineWidth(1);
00797     sumfit[i]->SetLineColor(cbgcolor);
00798 
00799 
00800     ptslice_bgsub[i]=new TH1D(*ptslice[i]);
00801     ptslice_bgsub[i]->SetName(subname);
00802     if(!isMC||isPythia){
00803       ptslice_bgsub[i]->Add(combbg[i],-1.0);
00804       combbg[i]->Draw("same");
00805     }
00806     
00807     //fit peak
00808     pi0peak[i]=new TF1(peakname,"gaus",0.05,0.2);
00809     pi0peak[i]->SetParLimits(1,.1,.2);
00810     
00811     if(!isMC||isPythia) ptslice_bgsub[i]->Fit(peakname,"QR0");
00812     else ptslice[i]->Fit(peakname,"QR0");
00813     pi0peak[i]->SetLineColor(peakcolor); 
00814     pi0peak[i]->SetLineStyle(3);
00815     pi0peak[i]->SetLineWidth(2);
00816     pi0peak[i]->SetRange(0.0,0.8);
00817 
00818     //draw sum of comb. + sign.
00819     Double_t params[6];
00820     combbg[i]->GetParameters(&params[0]);
00821     pi0peak[i]->GetParameters(&params[3]);
00822     sumfit[i]->SetParameters(params);
00823     if(!isMC||isPythia){
00824       sumfit[i]->Draw("same");
00825     }
00826 
00827     ptslice[i]->SetTitle(title);
00828     TAxis *ax=ptslice[i]->GetXaxis();
00829     ax->SetRangeUser(0.0,.8);
00830     if(strcmp(flag,"mb")==0) ax->SetRangeUser(0.0,0.4);
00831 
00832     if(!isMC||isPythia){
00833       Double_t MINIMUM=ptslice_bgsub[i]->GetMinimum();
00834       ptslice[i]->SetMinimum(MINIMUM);
00835     }
00836     
00837     //get 1st yield
00838     Float_t mean=pi0peak[i]->GetParameter(1);
00839     Float_t meanerr=pi0peak[i]->GetParError(1);
00840     Float_t sigma=pi0peak[i]->GetParameter(2);
00841     Float_t sigmaerr=pi0peak[i]->GetParError(2);
00842     
00843     Float_t yield=pi0peak[i]->Integral(mean-2.*sigma,mean+2.*sigma);
00844     yield/=ptslice_bgsub[i]->GetBinWidth(1);
00845     Float_t erryield=0.;
00846     
00847     if(yield>0.) erryield=sqrt(yield);
00848     else {yield=0.;erryield=0.;}
00849     
00850     x[i-1]=h->GetXaxis()->GetBinCenter(i);
00851     ex[i-1]=0.;
00852     y[i-1]=yield/h->GetXaxis()->GetBinWidth(i);// divide dpt
00853     ey[i-1]=erryield/h->GetXaxis()->GetBinWidth(i);
00854     m[i-1]=mean;
00855     em[i-1]=meanerr;
00856     s[i-1]=sigma;
00857     es[i-1]=sigmaerr;
00858     
00859   }
00860   c->Update();
00861   sprintf(canvasname,"c_%s_page_%d",flag,page++);
00862   c_copy=(TCanvas*)c->Clone();
00863   c_copy->SetName(canvasname);
00864   c_array->Add(c_copy);
00865   ps->NewPage();
00866   
00867   Float_t pTmin=1.;
00868   Float_t pTmax=7.;
00869   Float_t thres=0;
00870   if(strcmp(flag,"ht1")==0) {pTmin=3.;pTmax=11.;thres=2.5;}
00871   if(strcmp(flag,"ht2")==0) {pTmin=5.;pTmax=15.;thres=4.5;}
00872   //fit mean
00873   TF1 *meanfit=new TF1("meanfit","[0]*(1.0 - exp(-[1]*(x-[2])))");
00874   meanfit->SetParameter(0,.145);
00875   meanfit->SetParameter(1,2.);
00876   meanfit->SetParameter(2,thres);
00877   meanfit->SetParLimits(2,thres-2.0,10.);
00878   meanfit->SetRange(pTmin,pTmax);
00879   meanfit->SetLineColor(color);
00880   meanfit->SetLineStyle(2);
00881   //fit sigma
00882   TF1 *sigmafit=new TF1("sigmafit","[0]+[1]*exp(-[2]*(x-[3]))");
00883   sigmafit->SetParameter(0,0.02);
00884   sigmafit->SetParameter(1,0.05);
00885   sigmafit->SetParameter(2,1.0);
00886   sigmafit->SetParameter(3,thres);
00887   sigmafit->SetParLimits(3,thres-2.0,10.);
00888   sigmafit->SetRange(pTmin,pTmax);
00889   sigmafit->SetLineColor(color);
00890   sigmafit->SetLineStyle(2);
00891 
00892   TCanvas *c2=new TCanvas(TString("c2_")+TString(flag),TString("c2_")+TString(flag),600,600);
00893   c2->Divide(2,2);
00894   c2->cd(1);
00895 
00896   //draw mean,sigma,yield vs pT
00897   TGraphErrors *gry=new TGraphErrors(nx,x,y,ex,ey);
00898   gry->SetName((TString(flag)+"_yieldfit").Data());
00899   gPad->SetLogy();
00900   gry->SetMarkerStyle(25);
00901   gry->SetMarkerSize(.9);
00902   gry->SetMarkerColor(color);
00903   gry->GetXaxis()->SetRangeUser(pTmin,pTmax);
00904   gry->Draw("ap");
00905   gry->SetMaximum(10.*gry->GetHistogram()->GetMaximum());
00906   Float_t minimum=.1*gry->GetHistogram()->GetMinimum();
00907   gry->SetMinimum(minimum > 1. ? minimum : 1.);
00908 
00909   c2->cd(2);
00910   TGraphErrors *grm=new TGraphErrors(nx,x,m,ex,em);
00911   grm->SetName((TString(flag)+"_meanfit").Data());
00912   grm->SetMarkerStyle(25);
00913   grm->SetMarkerSize(.9);
00914   grm->SetMarkerColor(color);
00915   grm->GetXaxis()->SetRangeUser(pTmin,pTmax);
00916   grm->Draw("ap");
00917   grm->Fit("meanfit","QR");
00918   meanfit->Draw("same");
00919   c2->cd(3);
00920   TGraphErrors *grs=new TGraphErrors(nx,x,s,ex,es);
00921   grs->SetName((TString(flag)+"_sigmafit").Data());
00922   grs->SetMarkerStyle(25);
00923   grs->SetMarkerSize(.9);
00924   grs->SetMarkerColor(color);
00925   grs->SetMaximum(0.15);
00926   grs->SetMinimum(0.);
00927   grs->GetXaxis()->SetRangeUser(pTmin,pTmax);
00928   grs->Draw("ap");
00929   grs->Fit("sigmafit","QR");
00930   c2->cd(0);
00931 
00932   //c2->Update();
00933   //ps->NewPage();
00934   delete c2;
00935 
00936   TCanvas *cr=new TCanvas(TString("cr_")+TString(flag),TString("cr_")+TString(flag),1000,1000);
00937   cr->Divide(4,5);
00938 
00939   //re-integrate peak region (with correct error..):
00940   Float_t *x2=new Float_t[nx];
00941   Float_t *ex2=new Float_t[nx];
00942   Float_t *y2=new Float_t[nx];
00943   Float_t *ey2=new Float_t[nx];
00944   Float_t *m2=new Float_t[nx];
00945   Float_t *em2=new Float_t[nx];
00946   Float_t *s2=new Float_t[nx];
00947   Float_t *es2=new Float_t[nx];
00948   Float_t *bgy=new Float_t[nx];
00949 
00950 
00951   TF1 *finalpeak[100];
00952   TH1D *intregion[100];
00953   Char_t newname[100];
00954   Char_t newname2[100];
00955   Char_t newtitle[100];
00956   Char_t newfit[100];
00957   for(Int_t g=1;g<=nx;g++)
00958     {
00959       sprintf(newname,"ptbin_%d_%s",g,flag);
00960       sprintf(newname,"ptbin2_%d_%s",g,flag);
00961       sprintf(newtitle,"%.2f < p_{T} < %.2f",h->GetXaxis()->GetBinLowEdge(g),h->GetXaxis()->GetBinUpEdge(g));
00962       sprintf(newfit,"fit_peak_%d_%s",g,flag);
00963       finalpeak[g]=new TF1(newfit,"gaus");
00964       //get pT
00965       Float_t xval=h->GetXaxis()->GetBinCenter(g);
00966       //limit peak to +/- 2 sigma from extrapolation
00967       Float_t left=meanfit->Eval(xval)-2.*sigmafit->Eval(xval);
00968       Float_t right=meanfit->Eval(xval)+2.*sigmafit->Eval(xval);
00969       
00970       finalpeak[g]->SetParLimits(1,left,right);
00971       //define fit range +/- 3*sigma
00972       left=meanfit->Eval(xval)-3*sigmafit->Eval(xval);
00973       right=meanfit->Eval(xval)+3*sigmafit->Eval(xval);
00974 
00975       finalpeak[g]->SetRange(left,right);
00976       finalpeak[g]->SetLineColor(color);
00977       finalpeak[g]->SetLineStyle(2);
00978       finalpeak[g]->SetLineWidth(1);
00979       cr->cd(g);
00980       ptslice_bgsub[g]->SetName(newname);
00981       ptslice_bgsub[g]->SetTitle(newtitle);
00982       ptslice_bgsub[g]->GetXaxis()->SetTitle("M_{inv} (GeV/c^{2})");
00983       ptslice_bgsub[g]->GetYaxis()->SetTitle("counts (10 MeV/c^{2})^{-1}");
00984       if(strcmp(flag,"mb")!=0){
00985         ptslice_bgsub[g]->GetYaxis()->SetTitle("counts (20 MeV/c^{2})^{-1}");
00986       }
00987       ptslice_bgsub[g]->SetFillColor(0);
00988       ptslice_bgsub[g]->SetLineColor(color);
00989       ptslice_bgsub[g]->Draw("hist");
00990       ptslice_bgsub[g]->GetXaxis()->SetRangeUser(0.0,.8);
00991       if(strcmp(flag,"mb")==0){
00992         ptslice_bgsub[g]->GetXaxis()->SetRangeUser(0.0,0.4);
00993         //fit residual:
00994       }
00995 
00996       ptslice_bgsub[g]->Fit(newfit,"QR0");
00997       finalpeak[g]->SetRange(0.0,1.0);
00998       finalpeak[g]->Draw("same");
00999 
01000 
01001       //get integration region
01002       intregion[g]=new TH1D(*ptslice_bgsub[g]);
01003       intregion[g]->SetName(newname2);
01004       intregion[g]->SetFillColor(fillcolor);
01005       Float_t LEFT=finalpeak[g]->GetParameter(1)-NSIGMALO*finalpeak[g]->GetParameter(2);
01006       Float_t RIGHT=finalpeak[g]->GetParameter(1)+(NSIGMAHI)*finalpeak[g]->GetParameter(2);
01007       if(isMC) RIGHT=finalpeak[g]->GetParameter(1)+(NSIGMAHI+3)*finalpeak[g]->GetParameter(2);
01008       if(!isMC){
01009         if(isDAU){
01010           if(strcmp(flag,"mb")==0){
01011             if(RIGHT>0.2) RIGHT=0.2;
01012             if(LEFT<0.08) LEFT=0.081;
01013           }
01014           if(strcmp(flag,"ht1")==0&&xval>9.){
01015             LEFT=0.081;
01016             RIGHT=0.26;
01017           }
01018           if(strcmp(flag,"ht2")==0&&xval>9.){
01019             LEFT=0.081;
01020             if(xval>12.) RIGHT=0.3;
01021           }
01022         }
01023         
01024         if(isPP05){
01025           if(strcmp(flag,"mb")==0){
01026             //if(RIGHT>0.2) RIGHT=0.19;
01027             if(xval>2.) RIGHT=0.199;
01028             if(LEFT<0.08) LEFT=0.081;
01029           }
01030           if(strcmp(flag,"ht1")==0&&xval>9.){
01031             LEFT=0.081;
01032             RIGHT=0.241;
01033           }
01034           if(strcmp(flag,"ht2")==0&&xval>10.){
01035             LEFT=0.081;
01036           }
01037 
01038 
01039         }
01040       
01041       }
01042       
01043       
01044       x2[g-1]=xval;
01045       ex2[g-1]=0.;
01046             
01047       Float_t dpt=h->GetXaxis()->GetBinWidth(g);
01048       intregion[g]->GetXaxis()->SetRangeUser(LEFT,RIGHT);
01049     
01050 
01051       ptslice_bg[g]->GetXaxis()->SetRangeUser(LEFT,RIGHT);
01052       
01053       y2[g-1]=intregion[g]->Integral();
01054       bgy[g-1]=ptslice_bg[g]->Integral() - intregion[g]->Integral();
01055       if(!isMC) ey2[g-1]=sqrt(y2[g-1]+2.*bgy[g-1]);
01056       else{
01057         //calculate montecarlo error:
01058         ey2[g-1]=0;
01059         for(Int_t i=1;i<=intregion[g]->GetNbinsX();i++){
01060           ey2[g-1]+=pow(intregion[g]->GetBinError(i),2);
01061         }
01062         ey2[g-1]=sqrt(ey2[g-1]);
01063       } 
01064       y2[g-1]/=dpt;
01065       ey2[g-1]/=dpt;
01066       //cout<<xval<<"\t"<<y2[g-1]<<" +/- "<<ey2[g-1]<<endl;
01067       if(strcmp(flag,"mb")==0) cout<<xval<<"\t"<<dpt*y2[g-1]/bgy[g-1]<<"\t"<<ey2[g-1]/y2[g-1]<<endl;
01068 
01069       intregion[g]->Draw("histsame");
01070 
01071       //get mean from final gaus fit
01072       m2[g-1]=finalpeak[g]->GetParameter(1);
01073       em2[g-1]=finalpeak[g]->GetParError(1);
01074       s2[g-1]=finalpeak[g]->GetParameter(2);
01075       es2[g-1]=finalpeak[g]->GetParError(2);
01076 
01077       retYield->SetBinContent(g,y2[g-1]);
01078       retYield->SetBinError(g,ey2[g-1]);
01079 
01080     }
01081   cr->cd(0);
01082   cr->Update();
01083   sprintf(canvasname,"c_%s_page_%d",flag,page++);
01084   c_copy=(TCanvas*)cr->Clone();
01085   c_copy->SetName(canvasname);
01086   c_array->Add(c_copy);
01087   ps->NewPage();  
01088   
01089 
01090   TCanvas *c4=new TCanvas(TString("c4_")+TString(flag),TString("c4_")+TString(flag),600,900);
01091   c4->Divide(2,3);
01092   c4->cd(1);
01093   TGraphErrors *regry=new TGraphErrors(nx,x2,y2,ex2,ey2);
01094   regry->SetName((TString(flag)+"_yield").Data());
01095   gPad->SetLogy();
01096   regry->SetMarkerStyle(21);
01097   regry->SetMarkerSize(.9);
01098   regry->SetMarkerColor(color);
01099   regry->GetXaxis()->SetRangeUser(pTmin,pTmax);
01100   gry->GetXaxis()->SetRangeUser(pTmin,pTmax);
01101   Float_t maxy=gry->GetHistogram()->GetMaximum();
01102   Float_t miny=gry->GetHistogram()->GetMinimum() + 1.;
01103   if(isMC) miny=gry->GetHistogram()->GetMinimum() + 0.001;
01104   if(isPythia) miny=gry->GetHistogram()->GetMinimum();
01105   regry->SetMaximum(maxy);
01106   regry->SetMinimum(miny);
01107   gry->SetMaximum(maxy);
01108   gry->SetMinimum(miny);
01109   regry->Draw("ap");
01110   c4->cd(2);
01111   gPad->SetLogy();
01112   gry->Draw("ap");
01113   c4->cd(3);
01114   TGraphErrors *regrm=new TGraphErrors(nx,x2,m2,ex2,em2);
01115   regrm->SetName((TString(flag)+"_mean").Data());
01116   regrm->SetMarkerStyle(21);
01117   regrm->SetMarkerSize(.9);
01118   regrm->SetMarkerColor(color);
01119   regrm->GetXaxis()->SetRangeUser(pTmin,pTmax);
01120   grm->GetXaxis()->SetRangeUser(pTmin,pTmax);
01121   regrm->SetMaximum(.2);
01122   regrm->SetMinimum(.1);
01123   grm->SetMaximum(.2);
01124   grm->SetMinimum(.1);
01125   regrm->Draw("ap");
01126   c4->cd(4);
01127   grm->Draw("ap");
01128   c4->cd(5);
01129   TGraphErrors *regrs=new TGraphErrors(nx,x2,s2,ex2,es2);
01130   regrs->SetName((TString(flag)+"_sigma").Data());
01131   regrs->SetMarkerStyle(21);
01132   regrs->SetMarkerSize(.9);
01133   regrs->SetMarkerColor(color);
01134   regrs->GetXaxis()->SetRangeUser(pTmin,pTmax);
01135   grs->GetXaxis()->SetRangeUser(pTmin,pTmax);
01136   regrs->SetMaximum(.05);
01137   regrs->SetMinimum(0.);
01138   grs->SetMaximum(.05);
01139   grs->SetMinimum(0.);
01140   regrs->Draw("ap");
01141   c4->cd(6);
01142   grs->Draw("ap");
01143   c4->cd(0);
01144   c4->Update();
01145 
01146   sprintf(canvasname,"c_%s_page_%d",flag,page++);
01147   c_copy=(TCanvas*)c4->Clone();
01148   c_copy->SetName(canvasname);
01149   c_array->Add(c_copy);
01150 
01151   ps->NewPage();
01152 
01153   if(strcmp(flag,"ht2")==0) ps->Close();
01154 
01155   return retYield;
01156 }
01157 
01158 Float_t Pi0Analysis::getWeight(Float_t pT){
01159   //weight for neutrons (Levy function):
01160   if(pT<0.5) pT=0.5;
01161   Float_t weight=pT/pow((1.+(sqrt(pT*pT+1.) - 1.)/(0.215*12.55)),12.55);
01162   return weight;
01163 }
01164 
01165 void Pi0Analysis::storeCanvases(const Char_t *f_c){
01166   TFile *f_canvas=new TFile(f_c,"RECREATE");
01167   c_array->Write();
01168   f_canvas->Close();
01169 }