StRoot: StEmcPool/StPhotonAnalysis/Efficiency.cxx Source File

00001 
00002 #include <assert.h>
00003 
00004 #include <TColor.h>
00005 #include <TMultiGraph.h>
00006 #include <TGraphErrors.h>
00007 #include <TStyle.h>
00008 #include <TMinuit.h>
00009 #include <TRandom.h>
00010 #include <TTree.h>
00011 #include <TF1.h>
00012 #include <TFile.h>
00013 #include <TCanvas.h>
00014 #include <TVector3.h>
00015 #include <TLorentzVector.h>
00016 #include <TMath.h>
00017 #include <Riostream.h>
00018 
00019 #include <StEmcPool/StPhotonCommon/MyEvent.h> 
00020 #include <StEmcPool/StPhotonCommon/MyPoint.h> 
00021 #include <StEmcPool/StPhotonCommon/MyMcTrack.h> 
00022 
00023 #include "AnaCuts.h"
00024 #include "Pi0Analysis.h"
00025 
00026 #include "Efficiency.h"
00027 
00028 ClassImp(Efficiency)
00029 
00030 Efficiency::Efficiency(const char *input,const char *dir,const char *flag)
00031 {
00032   gStyle->SetOptDate(0);
00033 
00034   dirout=TString(dir);
00035 
00036   mFile=new TFile(input,"OPEN");
00037   myEventTree=(TTree*)mFile->Get("mEventTree");
00038   ev=new MyEvent();
00039   myEventTree->SetBranchAddress("branch",&ev);
00040 
00041   pseff=new TPostScript((dirout+"_eff.ps").Data(),-111);
00042   //init cuts
00043   cuts=new AnaCuts(flag);
00044   cuts->printCuts();
00045 
00046   USEWEIGHT=kTRUE;
00047   USEPYTHIAWEIGHT=kFALSE;
00048   USEBBCSPREAD=kTRUE;
00049 
00050   mFlag=flag;
00051   isMC=kFALSE;
00052   isPythia=kFALSE;
00053   isDAU=kFALSE;
00054   isPP05=kFALSE;
00055   if(strcmp(mFlag,"dAu")==0) isDAU=kTRUE;
00056   if(strcmp(mFlag,"pp05")==0) isPP05=kTRUE;
00057 
00058   NEUTRONS=kFALSE;
00059   ANTINEUTRONS=kFALSE;
00060   PIONS=kFALSE;
00061   ETAS=kFALSE;
00062   CHARGEDPIONS=kFALSE;
00063   PHOTONS=kFALSE;
00064   KZEROLONG=kFALSE;
00065 
00066   nNeutrons=0;
00067 
00068   htitle=TString("eff. for ");
00069 
00070   cout<<endl<<"Efficiency constructed!"<<endl<<endl;
00071 }
00072 
00073 Efficiency::~Efficiency()
00074 {
00075   cout<<endl<<"Efficiency destructed!"<<endl<<endl;
00076 }
00077 Int_t Efficiency::init()
00078 {
00079   //acceptance
00080   h_genMB=new TH2F("h_genMB","generated vs y and p_{T}",20,0.,20.,30,-0.5,1.5);
00081   h_accMB=new TH2F("h_accMB","acceptance vs y and p_{T}",20,0.,20.,30,-0.5,1.5);
00082 
00083   //conversions
00084   h_convGen=new TH1F("h_convGen","generated vs #eta",45,-0.2,1.3);
00085   h_convConv=new TH1F("h_convConv","converted vs #eta",45,-0.2,1.3);
00086   h_convConvSSD=new TH1F("h_convConvSSD","converted vs rapidity",45,-0.2,1.3);
00087   h_convConvSVT=new TH1F("h_convConvSVT","converted vs rapidity",45,-0.2,1.3);
00088   h_convConvIFC=new TH1F("h_convConvIFC","converted vs rapidity",45,-0.2,1.3);
00089 
00090   h_convConvNotCtb=new TH1F("h_convConvNotCtb","converted vs #eta before CTB",45,-0.2,1.3);
00091   h_stopRadius=new TH1F("h_stopRadius","conversion radius",1000,0.,450.);
00092 
00093 
00094   //ht
00095   h_HT1adc_id=new TH2F("h_HT1adc_id","ht-1 adc vs id",2400,.5,2400.5,800,200.5,1000.5);
00096   h_HT2adc_id=new TH2F("h_HT2adc_id","ht-2 adc vs id",2400,.5,2400.5,800,200.5,1000.5);
00097 
00098   //split clusters:
00099   h_splitClusAll=new TH1F("h_splitClusAll","h_splitClusAll",20,0.0,10.0);
00100   h_splitClus=new TH1F("h_splitClus","h_splitClus",20,0.0,10.0);
00101 
00102   //neutrons:
00103   h_nbarDet=new TH1F("h_nbarDet","h_nbarDet",cuts->nPtBinsEffMB,cuts->ptBinsEffMB.GetArray());
00104   h_nbarIn=new TH1F("h_nbarIn","h_nbarIn",cuts->nPtBinsEffMB,cuts->ptBinsEffMB.GetArray());
00105 
00106   //QA
00107   h_mcdist=new TH1F("h_mcdist","distance to mctrack",1000,0.,5.);
00108   h_mcdist2D=new TH2F("h_mcdist2D","distance to mctrack vs pT",20,0.,20.,100,0.,.25);
00109   h_dist=new TH1F("h_dist","distance to charged track",1000,0.,200.);
00110   h_dist2D=new TH2F("h_dist2D","distance to charged track vs pT",20,0.,20.,200,0.,50.);
00111   //eff:
00112   h_inputMB=new TH1F("h_inputMB","mc input;p_{T};dN/dp_{T}",cuts->nPtBinsEffMB,cuts->ptBinsEffMB.GetArray());
00113   h_inputHT1=new TH1F("h_inputHT1","mc input;p_{T};dN/dp_{T}",cuts->nPtBinsEffHT1,cuts->ptBinsEffHT1.GetArray());
00114   h_inputHT2=new TH1F("h_inputHT2","mc input;p_{T};dN/dp_{T}",cuts->nPtBinsEffHT2,cuts->ptBinsEffHT2.GetArray());
00115 
00116   h_inputDaughtersMB=new TH1F("h_inputDaughtersMB","mc inputDaughters;p_{T};dN/dp_{T}",cuts->nPtBinsEffMB,cuts->ptBinsEffMB.GetArray());
00117   h_inputDaughtersHT1=new TH1F("h_inputDaughtersHT1","mc inputDaughters;p_{T};dN/dp_{T}",cuts->nPtBinsEffHT1,cuts->ptBinsEffHT1.GetArray());
00118   h_inputDaughtersHT2=new TH1F("h_inputDaughtersHT2","mc inputDaughters;p_{T};dN/dp_{T}",cuts->nPtBinsEffHT2,cuts->ptBinsEffHT2.GetArray());
00119 
00120   h_recoMB=new TH1F("h_recoMB","reco.;p_{T};dN/dp_{T}",cuts->nPtBinsEffMB,cuts->ptBinsEffMB.GetArray());
00121   h_recoHT1=new TH1F("h_recoHT1","reco.;p_{T};dN/dp_{T}",cuts->nPtBinsEffHT1,cuts->ptBinsEffHT1.GetArray());
00122   h_recoHT2=new TH1F("h_recoHT2","reco.;p_{T};dN/dp_{T}",cuts->nPtBinsEffHT2,cuts->ptBinsEffHT2.GetArray());
00123 
00124   h_recoDaughtersMB=new TH1F("h_recoDaughtersMB","recoDaughters.;p_{T};dN/dp_{T}",cuts->nPtBinsEffMB,cuts->ptBinsEffMB.GetArray());
00125   h_recoDaughtersHT1=new TH1F("h_recoDaughtersHT1","recoDaughters.;p_{T};dN/dp_{T}",cuts->nPtBinsEffHT1,cuts->ptBinsEffHT1.GetArray());
00126   h_recoDaughtersHT2=new TH1F("h_recoDaughtersHT2","recoDaughters.;p_{T};dN/dp_{T}",cuts->nPtBinsEffHT2,cuts->ptBinsEffHT2.GetArray());
00127 
00128   h_minvMB=new TH2F("h_minvMB","reco.;p_{T};m_{inv}",
00129                     cuts->nPtBinsEffMB,cuts->ptBinsEffMB.GetArray(),cuts->nMinvBinsMB,cuts->mInvLowMB,cuts->mInvHighMB);
00130   h_minvHT1=new TH2F("h_minvHT1","reco.;p_{T};m_{inv}",
00131                      cuts->nPtBinsEffHT1,cuts->ptBinsEffHT1.GetArray(),cuts->nMinvBinsHT1,cuts->mInvLowHT1,cuts->mInvHighHT1);
00132   h_minvHT2=new TH2F("h_minvHT2","reco.;p_{T};m_{inv}",
00133                      cuts->nPtBinsEffHT2,cuts->ptBinsEffHT2.GetArray(),cuts->nMinvBinsHT2,cuts->mInvLowHT2,cuts->mInvHighHT2);
00134 
00135   h_matrixMB=new TH2F("h_matrixMB","pT gen. vs reco. MB",80,0.,20.,80,0.,20.);
00136 
00137   h_etaphi=new TH2F("h_etaphi","eta/phi neutral points",300,0.,1.,1800,-TMath::Pi(),TMath::Pi());
00138   h_pythiaPions=new TH1F("h_pythiaPions","pythia pions vs pT",200,0.,50.);
00139   h_pythiaPhotons=new TH1F("h_pythiaPhotons","pythia photons vs pT",200,0.,50.);
00140   h_pythiaPartonPt=new TH1F("h_pythiaPartonPt","pT of pythia process",200,0.,50);
00141 
00142   h_clusterWidth=new TH2F("h_clusterWidth","width BSMD eta+phi",20,0.,20.,100,0.,0.05);
00143   h_energyRatio=new TH2F("h_energyRatio","energy ratio BSMD/BTOW",20,0.,20.,160,0.,8.);
00144   h_towclusRatio=new TH2F("h_towclusRatio","hightow/cluster",20,0.,20.,50,0.,1.);
00145 
00146   h_vzMB=new TH1F("h_vzMB","z-vertex MB",480,-120.,120.);
00147   h_vzMB->Sumw2();
00148 
00149   h_asymmMB=new TH2F("h_asymmMB","asymmetry of inv mass pairs MB",cuts->nPtBinsMB,cuts->ptBinsMB.GetArray(),40,0.0,1.0);
00150   h_asymmMB->Sumw2();
00151   h_asymmHT1=new TH2F("h_asymmHT1","asymmetry of inv mass pairs HT1",cuts->nPtBinsHT1,cuts->ptBinsHT1.GetArray(),40,0.0,1.0);
00152   h_asymmHT1->Sumw2();
00153   h_asymmHT2=new TH2F("h_asymmHT2","asymmetry of inv mass pairs HT2",cuts->nPtBinsHT2,cuts->ptBinsHT2.GetArray(),40,0.0,1.0);
00154   h_asymmHT2->Sumw2();
00155 
00156   h_smdeta1=new TH2F("h_smdeta1","n eta strips in cluster vs p_{T}",20,0.,10.,10,-0.5,9.5);
00157   h_smdphi1=new TH2F("h_smdphi1","n phi strips in cluster vs p_{T}",20,0.,10.,10,-0.5,9.5);
00158   h_smdeta2=new TH2F("h_smdeta2","eta energy ratio in cluster vs p_{T}",20,0.,10.,100,0.0,10.);
00159   h_smdphi2=new TH2F("h_smdphi2","phi energy ratio in cluster vs p_{T}",20,0.,10.,100,0.0,10.);
00160 
00161   h_energyeta=new TH1F("h_energyeta","energy BSMDE",40,0.,10.);
00162   h_energyphi=new TH1F("h_energyeta","energy BSMDP",40,0.,10.);
00163 
00164   h_pionsVsEtaMB=new TH1F("h_pionsVsEta","dN_{#pi^{0}}/dy data",40,-0.5,1.5);
00165   h_pionsVsEtaMB->Sumw2();
00166 
00167 
00168   h_smdeta1->Sumw2();
00169   h_smdphi1->Sumw2();
00170   h_smdeta2->Sumw2();
00171   h_smdphi2->Sumw2();
00172 
00173 
00174   h_pythiaPartonPt->Sumw2();
00175   h_inputMB->Sumw2();
00176   h_inputHT1->Sumw2();
00177   h_inputHT2->Sumw2();
00178   h_inputDaughtersMB->Sumw2();
00179   h_inputDaughtersHT1->Sumw2();
00180   h_inputDaughtersHT2->Sumw2();
00181   h_recoMB->Sumw2();
00182   h_recoHT1->Sumw2();
00183   h_recoHT2->Sumw2();
00184   h_recoDaughtersMB->Sumw2();
00185   h_recoDaughtersHT1->Sumw2();
00186   h_recoDaughtersHT2->Sumw2();
00187   h_minvMB->Sumw2();
00188   h_minvHT1->Sumw2();
00189   h_minvHT2->Sumw2();
00190   h_matrixMB->Sumw2();
00191   h_pythiaPions->Sumw2();
00192   h_clusterWidth->Sumw2();
00193   h_energyRatio->Sumw2();
00194   h_towclusRatio->Sumw2();
00195 
00196 
00197 
00198   return 1; 
00199 }
00200 Int_t Efficiency::make(Int_t evmax)
00201 {
00202   Float_t WEIGHT=0.;
00203   Int_t pid=0;
00204   Int_t i=0;
00205   while(myEventTree->GetEntry(i))
00206     {
00207       if(evmax&&i>evmax)
00208         {
00209           cout<<"reached evmax,"<<endl;
00210           cout<<"abort loop!"<<endl;
00211           break;      
00212         }
00213       if(i%10000==0) cout<<"processing event: "<<i<<endl;
00214 
00215 
00216       MyMcTrack *mcTr=0;
00217       Float_t ptInput=0;
00218       Float_t rapInput=0;
00219       if(!isPythia){
00220         mcTr=ev->getMcTrack();
00221         ptInput=mcTr->momentum().Pt();
00222         rapInput=mcTr->momentum().PseudoRapidity();
00223       }
00224 
00225       if(i==0 && !isPythia){
00226         cout<<"not pythia, plain MC"<<endl;
00227         pid=mcTr->id();
00228         if(pid==111){PIONS=kTRUE;htitle=TString("pions");}
00229         else if(pid==221) {ETAS=kTRUE;htitle=TString("etas");}     
00230         else if(pid==2112) {NEUTRONS=kTRUE;htitle=TString("neutrons");}
00231         else if(pid==-2112) {ANTINEUTRONS=kTRUE;htitle=TString("antineutrons");}
00232         else if(pid==211) {CHARGEDPIONS=kTRUE;htitle=TString("piplus");}
00233         else if(pid==22) {PHOTONS=kTRUE;htitle=TString("single gamma");}
00234         else if(pid==130) {KZEROLONG=kTRUE;htitle=TString("single kzeroL");}
00235         else assert(0);
00236 
00237         htitle+=TString(" in ");
00238         htitle+=TString(mFlag);
00239       }
00240       else if(isPythia){
00241         htitle=TString("pions pythia");
00242       }
00243 
00244       if(NEUTRONS||ANTINEUTRONS){
00245         TLorentzVector lv(mcTr->momentum().X(),mcTr->momentum().Y(),mcTr->momentum().Z(),sqrt(mcTr->momentum().Mag2()+1.));
00246         rapInput=lv.Rapidity();
00247       }
00248       
00249       //get vertex
00250       TVector3 vPos=ev->vertex();
00251       //get vertex weight dAu + pp separately:
00252       WEIGHT=getWeightVertex(vPos.Z());
00253       h_vzMB->Fill(vPos.Z(),WEIGHT);
00254 
00255 
00256       //add spread from bbc:
00257       if(USEBBCSPREAD && isPP05 && gRandom->Rndm()<0.36){
00258         double bbc_spread=gRandom->Gaus(0.,43.5);
00259         vPos.SetZ(vPos.Z()+bbc_spread);
00260       }
00261 
00262       //set trigger from simu for dAu
00263       if(isDAU){
00264         Int_t tr=1;
00265         if(ev->highTowerAdc()>cuts->ht1AdcCUT) tr+=2;
00266         if(ev->highTowerAdc()>cuts->ht2AdcCUT) tr+=4;
00267         ev->setTrigger(tr);
00268       }
00269       else if(!isPP05){
00270         cout<<"wrong flag"<<endl;
00271         assert(0);
00272       }
00273 
00274       if(ev->trigger()&2) h_HT1adc_id->Fill(ev->highTowerId(),ev->highTowerAdc());
00275       if(ev->trigger()&4) h_HT2adc_id->Fill(ev->highTowerId(),ev->highTowerAdc());
00276 
00277       
00278       //event cuts
00279       ev->setMomentumInTpc(999999.);
00280 
00281       if(isDAU){
00282         if(!cuts->isEventOK(ev,"dAu")) {i++;continue;}
00283       }
00284       else if(isPP05){
00285         if(!cuts->isEventOK(ev,"pp05")) {i++;continue;}
00286       }
00287       else assert(0);
00288 
00289 
00290       //event ok, count neutrons:
00291       nNeutrons=nNeutrons+1.;
00292       float pt_rnd=10.*gRandom->Rndm();
00293       h_nbarIn->Fill(pt_rnd,exp(-0.5*pt_rnd));
00294       h_nbarDet->Fill(mcTr->momentum().Pt());
00295 
00296       if(ETAS||PIONS){
00297         if(ev->numberOfMcPhotons()!=2 && !isPythia){
00298           i++;
00299           continue;
00300         }
00301       }
00302       
00303       //fill input:
00304       if(isPythia && 0){
00305         Float_t pypT=ev->partonPt();
00306         if(pypT<15.){
00307           WEIGHT=1.;
00308         }
00309         else if(pypT<25.){
00310           WEIGHT=(216000/208000)*(0.0003895/0.002228);
00311         }
00312         else if(pypT>25.){
00313           WEIGHT=(216000/208000)*(1.016e-005/0.002228);
00314         }
00315       }
00316       if(isPythia){
00317         h_pythiaPartonPt->Fill(ev->partonPt(),WEIGHT);
00318         for(Int_t it=0;it<ev->getMcPionArray()->GetEntries();it++){
00319           MyMcTrack *pyt=(MyMcTrack*)ev->getMcPionArray()->At(it);
00320           if(pyt->momentum().PseudoRapidity()<cuts->rapPionMinCUT||
00321              pyt->momentum().PseudoRapidity()>cuts->rapPionMaxCUT) continue;
00322           h_pythiaPions->Fill(pyt->momentum().Pt(),WEIGHT);
00323           if(ev->trigger()&1){
00324             h_inputMB->Fill(pyt->momentum().Pt(),WEIGHT);
00325           }
00326           if(ev->trigger()&1){
00327             h_inputHT1->Fill(pyt->momentum().Pt(),WEIGHT);
00328           }
00329           if(ev->trigger()&1){
00330             h_inputHT2->Fill(pyt->momentum().Pt(),WEIGHT);
00331           }
00332         }
00333         for(Int_t it=0;it<ev->getMcPhotonArray()->GetEntries();it++){
00334           MyMcTrack *pyt=(MyMcTrack*)ev->getMcPhotonArray()->At(it);
00335           if(pyt->momentum().PseudoRapidity()<cuts->rapPionMinCUT||
00336              pyt->momentum().PseudoRapidity()>cuts->rapPionMaxCUT) continue;
00337           h_pythiaPhotons->Fill(pyt->momentum().Pt(),WEIGHT);
00338           if(ev->trigger()&1){
00339             h_inputDaughtersMB->Fill(pyt->momentum().Pt(),WEIGHT);
00340           }
00341           if(ev->trigger()&1){
00342             h_inputDaughtersHT1->Fill(pyt->momentum().Pt(),WEIGHT);
00343           }
00344           if(ev->trigger()&1){
00345             h_inputDaughtersHT2->Fill(pyt->momentum().Pt(),WEIGHT);
00346           }
00347         }
00348       }
00349 
00350 
00351 
00352       if(NEUTRONS){
00353         WEIGHT*=getWeightNeutrons(ptInput);     
00354         if(rapInput>cuts->rapidityMinCUT&&rapInput<cuts->rapidityMaxCUT){
00355           if(ev->trigger()&1) h_inputMB->Fill(mcTr->momentum().Pt(),WEIGHT);
00356           if(ev->trigger()&1) h_inputHT1->Fill(mcTr->momentum().Pt(),WEIGHT);
00357           if(ev->trigger()&1) h_inputHT2->Fill(mcTr->momentum().Pt(),WEIGHT);
00358         }
00359       }
00360 
00361       if(ANTINEUTRONS || KZEROLONG){
00362         WEIGHT*=getWeightAntiNeutrons(ptInput);
00363         if(rapInput>cuts->rapidityMinCUT&&rapInput<cuts->rapidityMaxCUT){
00364           if(ev->trigger()&1) h_inputMB->Fill(mcTr->momentum().Pt(),WEIGHT);
00365           if(ev->trigger()&1) h_inputHT1->Fill(mcTr->momentum().Pt(),WEIGHT);
00366           if(ev->trigger()&1) h_inputHT2->Fill(mcTr->momentum().Pt(),WEIGHT);
00367         }
00368       }
00369 
00370       if(PIONS||ETAS){
00371         if(PIONS) WEIGHT*=getWeightPions(ptInput);
00372         else if(ETAS) WEIGHT*=getWeightEtas(ptInput);
00373 
00374         //digamma
00375         if(rapInput>cuts->rapPionMinCUT&&rapInput<cuts->rapPionMaxCUT){
00376           if(ev->trigger()&1) h_inputMB->Fill(mcTr->momentum().Pt(),WEIGHT);
00377           if(ev->trigger()&1) h_inputHT1->Fill(mcTr->momentum().Pt(),WEIGHT);
00378           if(ev->trigger()&1) h_inputHT2->Fill(mcTr->momentum().Pt(),WEIGHT);
00379         }
00380         MyMcTrack *daughterA=(MyMcTrack*)ev->getMcPhotonArray()->At(0);
00381         MyMcTrack *daughterB=(MyMcTrack*)ev->getMcPhotonArray()->At(1);
00382 
00383         //if(TMath::Abs(daughterA->energy()-daughterB->energy())/(daughterA->energy()+daughterB->energy())<0.7){
00384         //  i++; 
00385         //  continue;
00386         //}
00387 
00388         //for acceptance
00389         h_genMB->Fill(mcTr->momentum().Pt(),mcTr->momentum().PseudoRapidity());
00390         if(daughterA->position().PseudoRapidity() > cuts->etaMinCUT &&
00391            daughterA->position().PseudoRapidity() < cuts->etaMaxCUT){
00392           if(daughterB->position().PseudoRapidity() > cuts->etaMinCUT &&
00393              daughterB->position().PseudoRapidity() < cuts->etaMaxCUT){
00394             h_accMB->Fill(mcTr->momentum().Pt(),mcTr->momentum().PseudoRapidity());
00395           }
00396         }
00397 
00398         if(daughterA->momentum().PseudoRapidity() > cuts->rapidityMinCUT &&
00399            daughterA->momentum().PseudoRapidity() < cuts->rapidityMaxCUT){
00400           if(ev->trigger()&1) h_inputDaughtersMB->Fill(daughterA->momentum().Pt(),WEIGHT);
00401           if(ev->trigger()&1) h_inputDaughtersHT1->Fill(daughterA->momentum().Pt(),WEIGHT);
00402           if(ev->trigger()&1) h_inputDaughtersHT2->Fill(daughterA->momentum().Pt(),WEIGHT);
00403         }
00404         if(daughterB->momentum().PseudoRapidity() > cuts->rapidityMinCUT &&
00405            daughterB->momentum().PseudoRapidity() < cuts->rapidityMaxCUT){
00406           if(ev->trigger()&1) h_inputDaughtersMB->Fill(daughterB->momentum().Pt(),WEIGHT);
00407           if(ev->trigger()&1) h_inputDaughtersHT1->Fill(daughterB->momentum().Pt(),WEIGHT);
00408           if(ev->trigger()&1) h_inputDaughtersHT2->Fill(daughterB->momentum().Pt(),WEIGHT);
00409         }
00410       }
00411       
00412       if(CHARGEDPIONS){
00413         WEIGHT*=getWeightPions(ptInput);
00414         if(rapInput>cuts->rapidityMinCUT&&rapInput<cuts->rapidityMaxCUT){
00415           if(ev->trigger()&1) h_inputMB->Fill(mcTr->momentum().Pt(),WEIGHT);
00416           if(ev->trigger()&1) h_inputHT1->Fill(mcTr->momentum().Pt(),WEIGHT);
00417           if(ev->trigger()&1) h_inputHT2->Fill(mcTr->momentum().Pt(),WEIGHT);
00418         }
00419       }
00420 
00421       if(PHOTONS){
00422         WEIGHT*=getWeightPhotons(ptInput);
00423         if(rapInput>cuts->rapidityMinCUT&&rapInput<cuts->rapidityMaxCUT){
00424           if(ev->trigger()&1) h_inputMB->Fill(mcTr->momentum().Pt(),WEIGHT);
00425           if(ev->trigger()&1) h_inputHT1->Fill(mcTr->momentum().Pt(),WEIGHT);
00426           if(ev->trigger()&1) h_inputHT2->Fill(mcTr->momentum().Pt(),WEIGHT);
00427 
00428           h_stopRadius->Fill(mcTr->stopRadius());
00429           h_convGen->Fill(mcTr->momentum().PseudoRapidity(),getWeightVertex(vPos.Z()));
00430           if(mcTr->stopRadius()>0&&mcTr->stopRadius()<223.5){
00431             h_convConv->Fill(mcTr->momentum().PseudoRapidity(),getWeightVertex(vPos.Z()));
00432           }
00433           if(mcTr->stopRadius()>0&&mcTr->stopRadius()<200.){
00434             h_convConvNotCtb->Fill(mcTr->momentum().PseudoRapidity(),getWeightVertex(vPos.Z()));
00435           }
00436           //SVT
00437           if(mcTr->stopRadius()>1.&&mcTr->stopRadius()<20.){
00438             h_convConvSVT->Fill(mcTr->momentum().PseudoRapidity(),getWeightVertex(vPos.Z()));
00439           }
00440           //SSD
00441           if(mcTr->stopRadius()>20.&&mcTr->stopRadius()<40.){
00442             h_convConvSSD->Fill(mcTr->momentum().PseudoRapidity(),getWeightVertex(vPos.Z()));
00443           }
00444           //IFC
00445           if(mcTr->stopRadius()>40.&&mcTr->stopRadius()<55.){
00446             h_convConvIFC->Fill(mcTr->momentum().PseudoRapidity(),getWeightVertex(vPos.Z()));
00447           }
00448         }
00449       }
00450 
00451       //check split clusters:
00452       int nclusters=0;
00453 
00454       MyPoint *p;
00455       MyPoint *pp;
00456       TClonesArray *clA=(TClonesArray*)ev->getPointArray();
00457       for(Int_t j=0;j<clA->GetEntries();j++)
00458         {
00459           p=(MyPoint*)clA->At(j);
00460           TVector3 pPos=p->position();
00461           TVector3 pMom=pPos-vPos;
00462           pMom.SetMag(p->energy());
00463 
00464           h_dist2D->Fill(pMom.Pt(),TMath::Abs(p->distanceToTrack()));
00465 
00466           if(!cuts->isPointOK(p,vPos)) continue;
00467 
00468           //chaeck split clusters:
00469           nclusters++;
00470           if(ev->trigger()&1 && pMom.Pt()>1.){
00471             h_etaphi->Fill(pPos.PseudoRapidity(),pPos.Phi(),WEIGHT);
00472           }
00473           Float_t PdistMC=9999.;
00474 
00475           h_dist->Fill(TMath::Abs(p->distanceToTrack()));
00476 
00477           MyMcTrack *closestTrack=0;
00478           if(PIONS||ETAS){
00479             PdistMC=9999.;//reset
00480             MyMcTrack *trA=(MyMcTrack*)ev->getMcPhotonArray()->At(0);
00481             MyMcTrack *trB=(MyMcTrack*)ev->getMcPhotonArray()->At(1);
00482             TVector3 mcPosA=trA->position();
00483             TVector3 mcPosB=trB->position();
00484             Float_t An=mcPosA.PseudoRapidity();
00485             Float_t Ap=mcPosA.Phi();
00486             Float_t Bn=mcPosB.PseudoRapidity();
00487             Float_t Bp=mcPosB.Phi();
00488             //point p:
00489             Float_t Pn=pPos.PseudoRapidity();
00490             Float_t Pp=pPos.Phi();
00491             Float_t dphiAP=TMath::Abs(Ap-Pp);
00492             while(dphiAP>2.*TMath::Pi()) dphiAP-=2.*TMath::Pi();
00493             Float_t dphiBP=TMath::Abs(Bp-Pp);
00494             while(dphiBP>2.*TMath::Pi()) dphiBP-=2.*TMath::Pi();
00495             Float_t dpA=sqrt(pow(An-Pn,2)+pow(dphiAP,2));
00496             Float_t dpB=sqrt(pow(Bn-Pn,2)+pow(dphiBP,2));
00497             //distance to mc track:
00498             PdistMC=dpA<dpB ? dpA : dpB;
00499             closestTrack=dpA<dpB ? trA : trB;
00500 
00501             h_mcdist->Fill(PdistMC);
00502             h_mcdist2D->Fill(pMom.Pt(),PdistMC);
00503           }
00504           
00505           //gamma:
00506           if(p->distanceToTrack()>cuts->photonCUT){
00507 
00508             if(pMom.PseudoRapidity()>cuts->rapidityMinCUT && pMom.PseudoRapidity()<cuts->rapidityMaxCUT){
00509               if(ETAS||PIONS||isPythia){
00510                 
00511                 if(ev->trigger()&1) h_recoDaughtersMB->Fill(pMom.Pt(),WEIGHT);
00512                 if(ev->trigger()&2) h_recoDaughtersHT1->Fill(pMom.Pt(),WEIGHT);
00513                 if(ev->trigger()&4) h_recoDaughtersHT2->Fill(pMom.Pt(),WEIGHT);
00514                 
00515               }
00516               else{
00517                 if(ev->trigger()&1) h_recoMB->Fill(pMom.Pt(),WEIGHT);
00518                 if(ev->trigger()&2) h_recoHT1->Fill(pMom.Pt(),WEIGHT);
00519                 if(ev->trigger()&4) h_recoHT2->Fill(pMom.Pt(),WEIGHT);
00520 
00521                 if(ev->trigger()&1) h_matrixMB->Fill(pMom.Pt(),ptInput);
00522               }
00523               
00524               if(ev->trigger()&1){
00525                 h_clusterWidth->Fill(pMom.Pt(),sqrt(p->widthEta()*p->widthEta()+p->widthPhi()*p->widthPhi()),WEIGHT);
00526                 h_energyRatio->Fill(pMom.Pt(),(p->energyEta()+p->energyPhi())/p->energy(),WEIGHT);
00527                 h_towclusRatio->Fill(pMom.Pt(),p->towerClusterEnergy(0)/p->energy(),WEIGHT);
00528                 
00529                 if(p->distanceToTrack()>30.){
00530                   h_smdeta1->Fill(pMom.Pt(),p->nHitsEta(),WEIGHT);
00531                   h_smdphi1->Fill(pMom.Pt(),p->nHitsPhi(),WEIGHT);
00532                   h_smdeta2->Fill(pMom.Pt(),p->energyEta()/p->energy(),WEIGHT);
00533                   h_smdphi2->Fill(pMom.Pt(),p->energyPhi()/p->energy(),WEIGHT);
00534                   
00535                   h_energyeta->Fill(p->energyEta());
00536                   h_energyphi->Fill(p->energyPhi());
00537                 }
00538               }
00539 
00540 
00541             }
00542          
00543           }
00544           
00545           for(Int_t jj=0;(ETAS||PIONS||PHOTONS||isPythia)&&jj<clA->GetEntries();jj++)
00546             {
00547               if(jj<=j) continue;
00548               pp=(MyPoint*)clA->At(jj);
00549               TVector3 ppPos=pp->position();
00550               TVector3 ppMom=ppPos-vPos;
00551               ppMom.SetMag(pp->energy());
00552               if(!cuts->isPointOK(pp,vPos)) continue;
00553 
00554               /*Float_t PPdistMC=9999.;
00555               if(!isPythia){
00556               MyMcTrack *trA=(MyMcTrack*)ev->getMcPhotonArray()->At(0);
00557               MyMcTrack *trB=(MyMcTrack*)ev->getMcPhotonArray()->At(1);
00558               TVector3 mcPosA=trA->position();
00559               TVector3 mcPosB=trB->position();
00560               Float_t An=mcPosA.PseudoRapidity();
00561               Float_t Ap=mcPosA.Phi();
00562               Float_t Bn=mcPosB.PseudoRapidity();
00563               Float_t Bp=mcPosB.Phi();  
00564               
00565               //point pp:
00566               Float_t PPn=ppPos.PseudoRapidity();
00567               Float_t PPp=ppPos.Phi();
00568               Float_t dphiAPP=TMath::Abs(Ap-PPp);
00569               while(dphiAPP>2.*TMath::Pi()) dphiAPP-=2.*TMath::Pi();
00570               Float_t dphiBPP=TMath::Abs(Bp-PPp);
00571               while(dphiBPP>2.*TMath::Pi()) dphiBPP-=2.*TMath::Pi();
00572               Float_t dppA=sqrt(pow(An-PPn,2)+pow(dphiAPP,2));
00573               Float_t dppB=sqrt(pow(Bn-PPn,2)+pow(dphiBPP,2));
00574               PPdistMC=dppA<dppB ? dppA : dppB;
00575               }*/
00576               
00577               //two neutrals:
00578               TVector3 pi0Mom=pMom+ppMom;
00579               Float_t angle=pMom.Angle(ppMom);
00580               Float_t minv=sqrt(2.*p->energy()*pp->energy()*(1. - cos(angle)));
00581               Float_t pTpion=pi0Mom.Pt();
00582               Float_t etapion=pi0Mom.PseudoRapidity();
00583               Float_t asymm=TMath::Abs(p->energy()-pp->energy())/(p->energy()+pp->energy());
00584 
00585               if(etapion<cuts->rapPionMinCUT||etapion>cuts->rapPionMaxCUT) continue;
00586               
00587               if(ev->trigger()&1){
00588                 if(asymm<=cuts->asymmetryCUT) h_minvMB->Fill(pTpion,minv,WEIGHT);
00589                 if(pTpion>1.){
00590                   if(minv>0.08 && minv<0.20){
00591                     h_asymmMB->Fill(pTpion,asymm,WEIGHT);
00592                     if(asymm<=cuts->asymmetryCUT){
00593                       h_pionsVsEtaMB->Fill(pi0Mom.PseudoRapidity(),WEIGHT);
00594                       h_matrixMB->Fill(pTpion,ptInput,WEIGHT);
00595                     }
00596                   }
00597                 }
00598               }
00599               if(ev->trigger()&2){
00600                 if(asymm<=cuts->asymmetryCUT) h_minvHT1->Fill(pTpion,minv,WEIGHT);
00601                 if(pTpion>4.){
00602                   if(minv>0.10 && minv<0.20){
00603                     h_asymmHT1->Fill(pTpion,asymm,WEIGHT);
00604                   }
00605                 }
00606               }
00607               if(ev->trigger()&4){
00608                 if(asymm<=cuts->asymmetryCUT) h_minvHT2->Fill(pTpion,minv,WEIGHT);
00609                 if(pTpion>6.){
00610                   if(minv>0.10 && minv<0.2){
00611                     h_asymmHT2->Fill(pTpion,asymm,WEIGHT);
00612                   }
00613                 }
00614               }
00615             
00616 
00617             }
00618         
00619 
00620         }
00621 
00622       if(PHOTONS && nclusters>0) h_splitClusAll->Fill(mcTr->momentum().Pt());
00623       if(PHOTONS && nclusters>1) h_splitClus->Fill(mcTr->momentum().Pt());
00624             
00625       i++;
00626     }
00627 
00628 
00629 
00630   if(ETAS||PIONS||isPythia){//extract yield with regular routine:
00631     Pi0Analysis *pi0=new Pi0Analysis((const char*)(dirout+"_effbins.ps").Data(),"/dev/null",mFlag);
00632     pi0->init("/dev/null");
00633     if(!isPythia) pi0->setMC(kTRUE);
00634     h_recoMB=new TH1F(*pi0->getYield(h_minvMB,"mb"));
00635     h_recoHT1=new TH1F(*pi0->getYield(h_minvHT1,"ht1"));
00636     h_recoHT2=new TH1F(*pi0->getYield(h_minvHT2,"ht2"));
00637     pi0->storeCanvases((dirout+"_canvases.root").Data());
00638     delete pi0;
00639     // !!!! calculate from dN/dpt to #/bin:
00640     for(Int_t i=1;i<=h_recoMB->GetNbinsX();i++)
00641       {
00642         h_recoMB->SetBinContent(i,h_recoMB->GetBinContent(i)*h_recoMB->GetBinWidth(i));
00643         h_recoMB->SetBinError(i,h_recoMB->GetBinError(i)*h_recoMB->GetBinWidth(i));
00644       }
00645     for(Int_t i=1;i<=h_recoHT1->GetNbinsX();i++)
00646       {
00647         h_recoHT1->SetBinContent(i,h_recoHT1->GetBinContent(i)*h_recoHT1->GetBinWidth(i));
00648         h_recoHT1->SetBinError(i,h_recoHT1->GetBinError(i)*h_recoHT1->GetBinWidth(i));
00649       }
00650     for(Int_t i=1;i<=h_recoHT2->GetNbinsX();i++)
00651       {
00652         h_recoHT2->SetBinContent(i,h_recoHT2->GetBinContent(i)*h_recoHT2->GetBinWidth(i));
00653         h_recoHT2->SetBinError(i,h_recoHT2->GetBinError(i)*h_recoHT2->GetBinWidth(i));
00654       }
00655   }
00656 
00657   gStyle->SetPalette(1);
00658   gStyle->SetStatStyle(0);
00659   
00660   pseff->On();
00661   TCanvas *c=new TCanvas("c","c",600,800);
00662   c->Divide(3,4);
00663   c->cd(1);
00664   gPad->SetLogy();
00665   h_inputMB->SetLineColor(1);
00666   h_inputMB->SetLineWidth(2);
00667   h_inputMB->SetTitle((htitle+TString(" MB")).Data());
00668   h_inputMB->Draw("hist");
00669   h_recoMB->SetLineColor(1);
00670   h_recoMB->SetLineWidth(2);
00671   h_recoMB->SetLineStyle(2);
00672   if(h_recoMB->GetMinimum()>0.) h_inputMB->SetMinimum(h_recoMB->GetMinimum()/10. );
00673   h_recoMB->Draw("histsame");
00674   c->cd(2);
00675   gPad->SetLogy();
00676   h_inputHT1->SetLineColor(TColor::GetColor(24,101,24));
00677   h_inputHT1->SetLineWidth(2);
00678   h_inputHT1->SetTitle((htitle+TString(" HT1")).Data());
00679   h_inputHT1->Draw("hist");
00680   h_recoHT1->SetLineColor(TColor::GetColor(24,101,24));
00681   h_recoHT1->SetLineWidth(2);
00682   h_recoHT1->SetLineStyle(2);
00683   if(h_recoHT1->GetMinimum()>0.) h_inputHT1->SetMinimum(h_recoHT1->GetMinimum()/10. );
00684   h_recoHT1->Draw("histsame");
00685   c->cd(3);
00686   gPad->SetLogy();
00687   h_inputHT2->SetLineColor(TColor::GetColor(24,28,174));
00688   h_inputHT2->SetLineWidth(2);
00689   h_inputHT2->SetTitle((htitle+TString(" HT2")).Data());
00690   h_inputHT2->Draw("hist");  
00691   h_recoHT2->SetLineColor(TColor::GetColor(24,28,174));
00692   h_recoHT2->SetLineWidth(2);
00693   h_recoHT2->SetLineStyle(2);
00694   if(h_recoHT2->GetMinimum()>0.) h_inputHT2->SetMinimum(h_recoHT2->GetMinimum()/10. );
00695   h_recoHT2->Draw("histsame");
00696 
00697   c->cd(4);
00698   if(NEUTRONS||ANTINEUTRONS||KZEROLONG){
00699     gPad->SetLogy();
00700   }
00701   h_effMB=new TH1F(*h_recoMB);
00702   h_effMB->SetNameTitle("h_effMB","efficiency MB;p_{T};#epsilon");
00703   h_effMB->SetLineStyle(1);
00704   h_effMB->Divide(h_inputMB);
00705   h_effMB->Draw("pe");
00706   c->cd(5);
00707   if(NEUTRONS||ANTINEUTRONS||KZEROLONG){
00708     gPad->SetLogy();
00709   }
00710   h_effHT1=new TH1F(*h_recoHT1);
00711   h_effHT1->SetNameTitle("h_effHT1","efficiency HT1;p_{T};#epsilon");
00712   h_effHT1->SetLineStyle(1);
00713   h_effHT1->Divide(h_inputHT1);
00714   h_effHT1->Draw("pe");
00715   c->cd(6);
00716   if(NEUTRONS||ANTINEUTRONS||KZEROLONG){
00717     gPad->SetLogy();
00718   }
00719   h_effHT2=new TH1F(*h_recoHT2);
00720   h_effHT2->SetNameTitle("h_effHT2","efficiency HT2;p_{T};#epsilon");
00721   h_effHT2->SetLineStyle(1);
00722   h_effHT2->Divide(h_inputHT2);
00723   h_effHT2->Draw("pe");
00724 
00725   c->cd(7);
00726   gPad->SetLogy();
00727   h_inputDaughtersMB->SetLineColor(1);
00728   h_inputDaughtersMB->SetLineWidth(2);
00729   h_inputDaughtersMB->Draw("hist");
00730   h_recoDaughtersMB->SetLineColor(1);
00731   h_recoDaughtersMB->SetLineWidth(2);
00732   h_recoDaughtersMB->SetLineStyle(2);
00733   h_recoDaughtersMB->Draw("histsame");
00734 
00735   c->cd(8);
00736   gPad->SetLogy();
00737   h_inputDaughtersHT1->SetLineColor(TColor::GetColor(24,101,24));
00738   h_inputDaughtersHT1->SetLineWidth(2);
00739   h_inputDaughtersHT1->Draw("hist");
00740   h_recoDaughtersHT1->SetLineColor(TColor::GetColor(24,101,24));
00741   h_recoDaughtersHT1->SetLineWidth(2);
00742   h_recoDaughtersHT1->SetLineStyle(2);
00743   h_recoDaughtersHT1->Draw("histsame");
00744   c->cd(9);
00745   gPad->SetLogy();
00746   h_inputDaughtersHT2->SetLineColor(TColor::GetColor(24,28,174));
00747   h_inputDaughtersHT2->SetLineWidth(2);
00748   h_inputDaughtersHT2->Draw("hist");
00749   h_recoDaughtersHT2->SetLineColor(TColor::GetColor(24,28,174));
00750   h_recoDaughtersHT2->SetLineWidth(2);
00751   h_recoDaughtersHT2->SetLineStyle(2);
00752   h_recoDaughtersHT2->Draw("histsame");
00753 
00754   c->cd(10);
00755   h_effDaughtersMB=new TH1F(*h_recoDaughtersMB);
00756   h_effDaughtersMB->SetNameTitle("h_effDaughtersMB","efficiency daughters MB;p_{T};#epsilon");
00757   h_effDaughtersMB->SetLineStyle(1);
00758   h_effDaughtersMB->Divide(h_inputDaughtersMB);
00759   h_effDaughtersMB->Draw("pe");
00760   c->cd(11);
00761   h_effDaughtersHT1=new TH1F(*h_recoDaughtersHT1);
00762   h_effDaughtersHT1->SetNameTitle("h_effDaughtersHT1","efficiency daughters HT1;p_{T};#epsilon");
00763   h_effDaughtersHT1->SetLineStyle(1);
00764   h_effDaughtersHT1->Divide(h_inputDaughtersHT1);
00765   h_effDaughtersHT1->Draw("pe");
00766   c->cd(12);
00767   h_effDaughtersHT2=new TH1F(*h_recoDaughtersHT2);
00768   h_effDaughtersHT2->SetNameTitle("h_effDaughtersHT2","efficiency daughters HT2;p_{T};#epsilon");
00769   h_effDaughtersHT2->SetLineStyle(1);
00770   h_effDaughtersHT2->Divide(h_inputDaughtersHT2);
00771   h_effDaughtersHT2->Draw("pe");
00772 
00773   
00774   c->cd(0);
00775   c->Update();
00776   c->SaveAs((dirout+"_effplots.root").Data());
00777   pseff->Close();
00778     
00779   
00780   return 1;
00781 }
00782 Int_t Efficiency::finish()
00783 {
00784   //save histograms
00785   mFileOut=new TFile((dirout+"_eff.root").Data(),"RECREATE");
00786   mFileOut->cd();
00787 
00788   h_genMB->Write();
00789   h_accMB->Write();
00790 
00791   h_HT1adc_id->Write();
00792   h_HT2adc_id->Write();
00793 
00794   h_splitClus->Write();
00795   h_splitClusAll->Write();
00796 
00797   h_convGen->Write();
00798   h_convConv->Write();
00799   h_convConvSSD->Write();
00800   h_convConvSVT->Write();
00801   h_convConvIFC->Write();
00802 
00803   h_convConvNotCtb->Write();
00804   h_stopRadius->Write();
00805 
00806   h_nbarDet->Write();
00807   h_nbarIn->Write();
00808 
00809   h_clusterWidth->Write();
00810   h_energyRatio->Write();
00811   h_towclusRatio->Write();
00812 
00813   h_smdeta1->Write();
00814   h_smdphi1->Write();
00815   h_smdeta2->Write();
00816   h_smdphi2->Write();
00817 
00818 
00819   h_pythiaPions->Write();
00820   h_pythiaPhotons->Write();
00821   h_pythiaPartonPt->Write();
00822   //
00823   h_mcdist->Write();
00824   h_mcdist2D->Write();
00825   h_dist->Write();
00826   h_dist2D->Write();
00827 
00828   h_effMB->Write();
00829   h_effHT1->Write();
00830   h_effHT2->Write();
00831 
00832   h_effDaughtersMB->Write();
00833   h_effDaughtersHT1->Write();
00834   h_effDaughtersHT2->Write();
00835 
00836   h_matrixMB->Write();
00837 
00838   h_etaphi->Write();
00839 
00840   h_minvMB->Write();
00841   h_minvHT1->Write();
00842   h_minvHT2->Write();
00843 
00844   h_asymmMB->Write();
00845   h_asymmHT1->Write();
00846   h_asymmHT2->Write();
00847 
00848   h_pionsVsEtaMB->Write();
00849 
00850   h_vzMB->Write();
00851 
00852   h_energyeta->Write();
00853   h_energyphi->Write();
00854 
00855   //
00856   mFileOut->Close();
00857   return 1;
00858 }
00859 Float_t Efficiency::getWeightPions(Float_t pT)
00860 {
00861   if(pT<0.2) pT=0.2;
00862   float weight=1.;
00863   if(isPP05){
00864     float p[]={592.,-9.31,5.03,-7.55,-1.7,6.06};
00865     float WS=1. - 1./(1.+TMath::Exp(TMath::Abs(p[4])*(pT-p[5])));
00866     weight=WS*p[2]*pow(pT,p[3]);
00867     weight+=(1.- WS)*p[0]*pow(1.+pT,p[1]);
00868     weight*=pT;
00869     weight=7.0e+05*pT*pow(1.+pT,-9.3);
00870   }
00871   else if(isDAU){
00872     float p[]={39.3,-8.58,0.828,-7.27,0.835,5.18};
00873     float WS=1. - 1./(1.+exp(p[4]*(pT-p[5])));
00874     weight=WS*p[2]*pow(1.+pT,p[3]);
00875     weight+=(1.- WS)*p[0]*pow(1.+pT,p[1]);
00876     weight*=pT;
00877     weight=7.0e+05*pT*pow(1.+pT,-9.3);
00878   }
00879 
00880   if(!USEWEIGHT) weight=1.;
00881   if(USEPYTHIAWEIGHT) weight=exp(-0.401*pT);
00882   return weight;
00883 }
00884 Float_t Efficiency::getWeightEtas(Float_t pT)
00885 {
00886   Float_t ME=0.548;//eta meson mass
00887   Float_t MP=0.135;//pi0 mass
00888   Float_t weight=0.45*pT/sqrt(pT*pT+ME*ME-MP*MP)*getWeightPions(sqrt(pT*pT+ME*ME-MP*MP));
00889   if(!USEWEIGHT) weight=1.;
00890   return weight;
00891 }
00892 Float_t Efficiency::getWeightAntiNeutrons(Float_t pT)
00893 {
00894   //hadron data for Levy function, nucl-ex/0601033:
00895   //-->  d^2N/(2*pi*pT*N*dpT*dy) = B/((1+((mT - m0)/nT))^n)
00896   // {p-dAu; pbar-dAu; p-pp; pbar-pp} and m0 = m_neutron = 1.0 GeV.
00897   //Double_t B[]={0.3,0.23,0.072,0.061};//->[0]
00898   //Double_t eB[]={0.01,0.01,0.005,0.005};
00899   Double_t T[]={0.205,0.215,0.179,0.173};//->[1]
00900   //Double_t eT[]={0.004,0.005,0.006,0.006};
00901   Double_t n[]={11.00,12.55,10.87,10.49};//->[2]
00902   //Double_t en[]={0.29,0.41,0.43,0.40};
00903   
00904   if(pT<0.5) pT=0.5;
00905   Float_t weight=0.;
00906   Int_t bin=-1;
00907   if(isDAU) bin=1;
00908   if(isPP05) bin=3;
00909   weight=pT/pow((1.+(sqrt(pT*pT+1.) - 1.)/(T[bin]*n[bin])),n[bin]);
00910   if(isPP05 && ANTINEUTRONS) weight*=exp(0.5*pT);
00911   if(KZEROLONG) weight=getWeightPions(pT);
00912   return weight;
00913 }
00914 Float_t Efficiency::getWeightNeutrons(Float_t pT)
00915 {
00916   //hadron data for Levy function, nucl-ex/0601033:
00917   //-->  d^2N/(2*pi*pT*N*dpT*dy) = B/((1+((mT - m0)/nT))^n)
00918   // {p-dAu; pbar-dAu; p-pp; pbar-pp} and m0 = m_neutron = 1.0 GeV.
00919   //Double_t B[]={0.3,0.23,0.072,0.061};//->[0]
00920   //Double_t eB[]={0.01,0.01,0.005,0.005};
00921   Double_t T[]={0.205,0.215,0.179,0.173};//->[1]
00922   //Double_t eT[]={0.004,0.005,0.006,0.006};
00923   Double_t n[]={11.00,12.55,10.87,10.49};//->[2]
00924   //Double_t en[]={0.29,0.41,0.43,0.40};
00925 
00926   if(pT<0.5) pT=0.5;
00927   Float_t weight=0.;
00928   Int_t bin=-1;
00929   if(isDAU) bin=0;
00930   if(isPP05) bin=2;
00931   weight=pT/pow((1.+(sqrt(pT*pT+1.) - 1.)/(T[bin]*n[bin])),n[bin]);
00932   return weight;
00933 }
00934 Float_t Efficiency::getWeightPhotons(Float_t pT)
00935 {
00936   if(pT<0.2) pT=0.2;
00937   return getWeightPions(pT);  
00938 }
00939 Float_t Efficiency::getWeightVertex(Float_t z)
00940 {
00941   Float_t weight=1.;
00942   if(isDAU){
00943     weight=1.04 + 0.00547*z - 0.00016*z*z - 2.8e-06*z*z*z;
00944     weight=weight + 4.7e-08*z*z*z*z + 4.5e-10*z*z*z*z*z;
00945   }
00946   if(isPP05){
00947     float sigma_data2=50.*50.;
00948     float sigma_mc2=60.*60.;
00949     float factor=(1./sigma_data2)-(1./sigma_mc2);
00950     weight=exp(-0.5*z*z*factor);    
00951   }
00952   return weight;
00953 }