+
class AliSpectraBothHistoManager;
class AliSpectraBothEventCuts;
class AliSpectraBothTrackCuts;
Int_t Marker[6]={20,21,22,24,25,26};
Double_t Range[3]={0.3,0.3,0.5}; // LowPt range for pi k p
+
+Double_t FitRange[2]={-2.0,2.0};
+Double_t CutRange[2]={-3.0,3.0};
+Double_t minptformaterial[6]={0.0,0.2,0.0,0.0,0.2,0.0};
+Double_t maxptformaterial[6]={0.0,0.6,1.3,0.0,0.6,0.0};
+Double_t minptforWD[6]={0.2,100.0,0.3,0.2,100.0,0.3};
+Double_t maxptforWD[6]={1.5,-100.0,2.0,1.5,-100.0,2.0};
+Double_t minRanges[3]={0.3,0.3,0.45};
+Double_t maxRanges[3]={1.5,1.2,2.2};
+Double_t fMaxContaminationPIDMC=0.2;
+
+
enum ECharge_t {
kPositive,
kNegative,
kNCharges
};
enum {
- kdodca=0x1,
- kgeantflukaKaon=0x2,
- kgeantflukaProton=0x4,
- knormalizationtoeventspassingPhySel=0x8,
- kveretxcorrectionandbadchunkscorr=0x10,
- kmcisusedasdata=0x20
+ kdodca=0x1, //dca fits are made
+ kgeantflukaKaon=0x2,// geant fluka correction is used for kaons
+ kgeantflukaProton=0x4, // geant fluka correction is used for protons and antiprotons
+ knormalizationtoeventspassingPhySel=0x8,// spectra are divided by number of events passing physic selection
+ kveretxcorrectionandbadchunkscorr=0x10, // correction for difference in z vertex distribution in data and MC and correction for bad chunks is applied
+ kmcisusedasdata=0x20, // the result of the looping over MC is used as data input
+ kdonotusedcacuts=0x40, // allows to use the constant dca cut for all pt bins not the pt dependet defined in stardrad track cuts 2011
+ kuseprimaryPIDcont=0x80, //pid contamination is calculated using only primiary particle in this case K should use dca fits
+ knormalizationwithbin0integralsdata=0x100, // the normalization factor is calcualte using integral over z vertex distributions (in this case reconstructed vertex disitrbution uses z vertex for data)
+ knormalizationwithbin0integralsMC=0x200, //in this case reconstructed vertex disitrbution uses z vertex for data, those to options will be use only if knormalizationtoeventspassingPhySel is not set
+ kuserangeonfigfile=0x400, // use of config file for dca fit settings
+ kskipconcutonspectra=0x800, //do not use conPID<02 cut useful for syst. studies
+ kuseTOFmatchingcorrection=0x1000 // if set tof matching correction is applied.
+
};
Bool_t OpenFile(TString dirname, TString outputname, Bool_t mcflag,Bool_t mcasdata=false);
-void AnalysisBoth (UInt_t options=0xF,TString outdate, TString outnamedata, TString outnamemc="" )
+void AnalysisBoth (UInt_t options=0xF,TString outdate, TString outnamedata, TString outnamemc="",TString configfile="" )
{
+ gStyle->SetOptStat(0);
TH1::AddDirectory(kFALSE);
gSystem->Load("libCore.so");
gSystem->Load("libPhysics.so");
mass[1] = TDatabasePDG::Instance()->GetParticle("K+")->Mass();
mass[2] = TDatabasePDG::Instance()->GetParticle("proton")->Mass();
- AliESDtrackCuts* esdtrackcuts= AliESDtrackCuts::GetStandardITSTPCTrackCuts2011(kTRUE);
- TString formulastring(esdtrackcuts->GetMaxDCAToVertexXYPtDep());
- formulastring.ReplaceAll("pt","x");
- TFormula* dcacutxy=new TFormula("dcacutxy",formulastring.Data());
-
+ TFormula* dcacutxy=0x0;
+ if(!(options&kdonotusedcacuts))
+ {
+
+ AliESDtrackCuts* esdtrackcuts= AliESDtrackCuts::GetStandardITSTPCTrackCuts2011(kTRUE);
+ TString formulastring(esdtrackcuts->GetMaxDCAToVertexXYPtDep());
+ formulastring.ReplaceAll("pt","x");
+ dcacutxy=new TFormula("dcacutxy",formulastring.Data());
+ }
+ if(options&kuserangeonfigfile)
+ if(!ReadConfigFile(configfile))
+ return;
TList* lout=new TList();
-
+
TString indirname=Form("/output/train%s",outdate.Data());
//TString indirname("/output/train24072012");
TH1F* contWD[6];
TH1F* contMat[6];
TH1F* confinal[6];
+ TH1F* contPIDpri[6];
+ TH1F* contSecMC[6];
TH1F* contfit[12];
TH1F* contWDfit[12];
GetPtHistFromPtDCAhisto("hHistPtRecTrue","conPID",managermc,contPID,dcacutxy);
GetPtHistFromPtDCAhisto("hHistPtRecSigmaSecondaryWeakDecay","conWD",managermc,contWD,dcacutxy);
GetPtHistFromPtDCAhisto("hHistPtRecSigmaSecondaryMaterial","conMat",managermc,contMat,dcacutxy);
+ GetPtHistFromPtDCAhisto("hHistPtRecSigmaPrimary","conPIDprimary",managermc,contPIDpri,dcacutxy);
+
-
-// Double_t neventsdata = ecutsdata->NumberOfEvents();
Double_t neventsmcall = 1 ; //if loop over MC is done after or befor events cuts this will be changed
Double_t neventsdata = 1;
Double_t neventsmc = 1;
+ //Normaliztion of MCtruth depends if the loop was done after of before ESD event cuts.
+ //In currect code this cannot be check on the level of macro.
+ //If the loop was done before MC should be done to all processed events (NumberOfProcessedEvents())
+ //If loop was done after MC should be normalized to all accepted events (NumberOfEvents())
+ // The option one will be alaways use.
+
+ neventsmcall= ecutsmc->NumberOfProcessedEvents();
if(options&knormalizationtoeventspassingPhySel)
{
- neventsmcall= ecutsmc->NumberOfProcessedEvents();
+ //neventsmcall= ecutsmc->NumberOfProcessedEvents();
neventsdata=ecutsdata->NumberOfPhysSelEvents();
neventsmc=ecutsmc->NumberOfPhysSelEvents();
}
+ else if ((options&knormalizationwithbin0integralsdata)||(options&knormalizationwithbin0integralsMC))
+ {
+ neventsdata=Normaliztionwithbin0integrals(options);
+ neventsmc=ecutsmc->NumberOfPhysSelEvents();
+ }
else
{
- neventsdata=ecutsdata->NumberOfEvents();
+ neventsdata=ecutsdata->NumberOfEvents(); //number of accepted events
neventsmc=ecutsmc->NumberOfEvents();
neventsmcall= ecutsmc->NumberOfEvents();
-
}
GetMCTruth(MCTruth);
-
+ cout<<neventsdata<<" Events"<<endl;
TH1F* allgen=((TH1F*)managermc->GetPtHistogram1D("hHistPtGen",1,1))->Clone();
TH1F* spectraall=(TH1F*)allrecdata->Clone("recNch");
+ spectraall->SetTitle("recNch");
TH1F* contall=(TH1F*)allrecMC->Clone("contall");
- contall->Add(alleff,-1);
+ contall->SetTitle("contall");
+ //contall->Add(alleff,-1);
+ SubHistWithFullCorr(contall,alleff);
alleff->Divide(alleff,allgen,1,1,"B");
contall->Divide(contall,allrecMC,1,1,"B");
lout->Add(alleff);
lout->Add(allrecdata);
lout->Add(spectraall);
- lout->Add(contall);
+ lout->Add(contall);
+
for (int i=0;i<6;i++)
{
TString tmpname(rawspectramc[i]->GetTitle());
tmpname.ReplaceAll("SpectraMC","%s");
- contallMC[i]=(TH1F*)rawspectramc[i]->Clone(Form(tmpname.Data(),"contallMC"));
+ contallMC[i]=(TH1F*)rawspectramc[i]->Clone(Form(tmpname.Data(),"contallMC"));
+ contallMC[i]->SetTitle(Form(tmpname.Data(),"contallMC"));
contfit[i]=(TH1F*)rawspectramc[i]->Clone(Form(tmpname.Data(),"contfit"));
+ contfit[i]->SetTitle(Form(tmpname.Data(),"contfit"));
contWDfit[i]=(TH1F*)rawspectramc[i]->Clone(Form(tmpname.Data(),"contWDfit"));
+ contWDfit[i]->SetTitle(Form(tmpname.Data(),"contWDfit"));
contMatfit[i]=(TH1F*)rawspectramc[i]->Clone(Form(tmpname.Data(),"contMatfit"));
+ contMatfit[i]->SetTitle(Form(tmpname.Data(),"contMatfit"));
primaryfit[i]=(TH1F*)rawspectramc[i]->Clone(Form(tmpname.Data(),"primaryfit"));
-
+ primaryfit[i]->SetTitle(Form(tmpname.Data(),"primaryfit"));
contfit[i+6]=(TH1F*)rawspectramc[i]->Clone(Form(tmpname.Data(),"contfitonMC"));
+ contfit[i+6]->SetTitle(Form(tmpname.Data(),"contfitonMC"));
contWDfit[i+6]=(TH1F*)rawspectramc[i]->Clone(Form(tmpname.Data(),"contWDfitonMC"));
+ contWDfit[i+6]->SetTitle(Form(tmpname.Data(),"contWDfitonMC"));
contMatfit[i+6]=(TH1F*)rawspectramc[i]->Clone(Form(tmpname.Data(),"contMatfitonMC"));
+ contMatfit[i+6]->SetTitle(Form(tmpname.Data(),"contMatfitonMC"));
primaryfit[i+6]=(TH1F*)rawspectramc[i]->Clone(Form(tmpname.Data(),"primaryfitMC"));
+ primaryfit[i+6]->SetTitle(Form(tmpname.Data(),"primaryfitMC"));
+
+
contfit[i]->Reset();
contWDfit[i]->Reset();
contMatfit[i]->Reset();
SetBintoOne(primaryfit[i]);
SetBintoOne(primaryfit[i+6]);
spectra[i]=(TH1F*)rawspectradata[i]->Clone(Form(tmpname.Data(),"SpectraFinal"));
-
+ spectra[i]->SetTitle(Form(tmpname.Data(),"SpectraFinal"));
spectraLeonardo[i]=(TH1F*)rawspectradata[i]->Clone(Form(tmpname.Data(),"SpectraFinalLeonardo"));
+ spectraLeonardo[i]->SetTitle(Form(tmpname.Data(),"SpectraFinalLeonardo"));
+
corrLeonardo[i]=(TH1F*)MCTruth[i]->Clone(Form(tmpname.Data(),"CorrFactLeonardo"));
-
+ corrLeonardo[i]->SetTitle(Form(tmpname.Data(),"CorrFactLeonardo"));
corrLeonardo[i]->Divide(corrLeonardo[i],rawspectramc[i],1,1,"B");
- contallMC[i]->Add(eff[i],-1.0);
- RecomputeErrors(contallMC[i]);
+ //contallMC[i]->Add(eff[i],-1.0);
+ SubHistWithFullCorr(contallMC[i],eff[i]);
+ //RecomputeErrors(contallMC[i]);
contallMC[i]->Sumw2();
contallMC[i]->Divide(contallMC[i],rawspectramc[i],1,1,"B");
-
+
+ // contamintaion from PID but only primaries
+ //contPIDpri[i]->Add(eff[i],-1.0);
+ SubHistWithFullCorr(contPIDpri[i],eff[i]);
+
+ //RecomputeErrors(contPIDpri[i]);
+ contPIDpri[i]->Divide(contPIDpri[i],rawspectramc[i],1,1,"B");
+
eff[i]->Divide(eff[i],MCTruth[i],1,1,"B");
contPID[i]->Sumw2();
rawspectramc[i]->Sumw2();
- contPID[i]->Add(contPID[i],rawspectramc[i],-1,1);
- RecomputeErrors(contPID[i]);
+ //contPID[i]->Add(contPID[i],rawspectramc[i],-1,1);
+ SubHistWithFullCorr(contPID[i],rawspectramc[i]);
+ contPID[i]->Scale(-1.0);
+
+ //RecomputeErrors(contPID[i]);
contPID[i]->ResetStats();
contPID[i]->Sumw2();
contPID[i]->Divide(contPID[i],rawspectramc[i],1,1,"B");
-
- confinal[i]=(TH1F*)contPID[i]->Clone(Form(tmpname.Data(),"confinal"));
-
+
+ if(options&kuseprimaryPIDcont)
+ confinal[i]=(TH1F*)contPIDpri[i]->Clone(Form(tmpname.Data(),"confinal"));
+ else
+ confinal[i]=(TH1F*)contPID[i]->Clone(Form(tmpname.Data(),"confinal"));
+ confinal[i]->SetTitle(Form(tmpname.Data(),"confinal"));
+ contSecMC[i]=(TH1F*)contWD[i]->Clone(Form(tmpname.Data(),"conSecMC"));
+ contSecMC[i]->Add(contMat[i]);
contWD[i]->Divide(contWD[i],rawspectramc[i],1,1,"B");
contMat[i]->Divide(contMat[i],rawspectramc[i],1,1,"B");
-
+ contSecMC[i]->Divide(contSecMC[i],rawspectramc[i],1,1,"B");
rawspectradata[i]->Scale(1./neventsdata,"width");
lout->Add(contfit[i]);
lout->Add(contWDfit[i]);
lout->Add(contMatfit[i]);
+ lout->Add(primaryfit[i]);
lout->Add(contfit[i+6]);
lout->Add(contWDfit[i+6]);
lout->Add(contMatfit[i+6]);
+ lout->Add(primaryfit[i+6]);
lout->Add(spectra[i]);
lout->Add(spectraLeonardo[i]);
lout->Add(confinal[i]);
+ lout->Add(contPIDpri[i]);
+ lout->Add(contSecMC[i]);
}
- TFile* fout=new TFile(Form("./results/ResMY_%s_%s.root",outnamemc.Data(),outdate.Data()),"RECREATE");
+ outdate.ReplaceAll("/","_");
+ configfile.ReplaceAll(".","_");
+ TFile* fout=0x0;
+ if(configfile.Length()>0&&(options&kuserangeonfigfile))
+ fout=new TFile(Form("./results/ResMY_%s_%s_%#X_%s.root",outnamemc.Data(),outdate.Data(),options,configfile.Data()),"RECREATE");
+ else
+ fout=new TFile(Form("./results/ResMY_%s_%s_%#X.root",outnamemc.Data(),outdate.Data(),options),"RECREATE");
if (options&kdodca)
DCACorrectionMarek(managerdata,managermc,dcacutxy,fout,contfit,contWDfit,contMatfit,primaryfit);
for (int i=0;i<6;i++)
{
if(options&kdodca)
{
- confinal[i]->Add(contfit[i]);
+ if((options&kuseprimaryPIDcont)||(i!=1&&i!=4)) //if we do not use cont PId only for primary and this is a kaon that do not use fit
+ confinal[i]->Add(contfit[i]);
GetCorrectedSpectra(spectra[i],rawspectradata[i],eff[i],confinal[i]);
}
else
GetCorrectedSpectra(spectra[i],rawspectradata[i],eff[i],contallMC[i]);
}
GetCorrectedSpectraLeonardo(spectraLeonardo[i],corrLeonardo[i],primaryfit[i],primaryfit[i+6]);
- CleanHisto(spectra[i],-1,100,contPID[i]);
- CleanHisto(spectraLeonardo[i],-1,100,contPID[i]);
+ if(options&kskipconcutonspectra)
+ continue;
+ if(options&kuseprimaryPIDcont)
+ {
+ CleanHisto(spectra[i],-1,100,contPIDpri[i]);
+ CleanHisto(spectraLeonardo[i],-1,100,contPIDpri[i]);
+ }
+ else
+ {
+ CleanHisto(spectra[i],-1,100,contPID[i]);
+ CleanHisto(spectraLeonardo[i],-1,100,contPID[i]);
+ }
}
GFCorrection(spectra,tcutsdata->GetPtTOFMatching(),options);
GFCorrection(spectraLeonardo,tcutsdata->GetPtTOFMatching(),options);
- MatchingTOFEff(spectra,lout);
- MatchingTOFEff(spectraLeonardo);
- TH1F* allch=GetSumAllCh(spectra,mass);
+ TH1F* allch=GetSumAllCh(spectra,mass);
lout->Add(allch);
-
+ if(options&kuseTOFmatchingcorrection)
+ {
+ MatchingTOFEff(spectra,lout);
+ MatchingTOFEff(spectraLeonardo);
+ TOFMatchingForNch(spectraall);
+ }
// lout->ls();
fout->cd();
TList* listqa=new TList();
TList* canvaslist=new TList();
Float_t vertexcorrection=1.0;
- Float_t corrbadchunksvtx=QAPlotsBoth(managerdata,managermc,ecutsdata,ecutsmc,tcutsdata,tcutsmc,listqa,canvaslist);
+ Float_t corrbadchunksvtx=1.0;
+ if ((options&knormalizationwithbin0integralsdata)||(options&knormalizationwithbin0integralsMC))
+ corrbadchunksvtx=QAPlotsBoth(managerdata,managermc,ecutsdata,ecutsmc,tcutsdata,tcutsmc,listqa,canvaslist,0);
+ else
+ corrbadchunksvtx=QAPlotsBoth(managerdata,managermc,ecutsdata,ecutsmc,tcutsdata,tcutsmc,listqa,canvaslist,1);
if (options&kveretxcorrectionandbadchunkscorr)
vertexcorrection=corrbadchunksvtx;
cout<<" VTX corr="<<vertexcorrection<<endl;
canvaslist->Write("outputcanvas",TObject::kSingleKey);
fout->Close();
+ //Normaliztionwithbin0integrals();
}
tcutsmc->PrintCuts();
if(!managermc||!ecutsmc||!tcutsmc)
return false;
+ if(managermc->GetGenMulvsRawMulHistogram("hHistGenMulvsRawMul")->GetEntries()!=ecutsmc->GetHistoCuts()->GetBinContent(3))
+ cout<<"Please check MC file possible problem with merging"<<" "<<managermc->GetGenMulvsRawMulHistogram("hHistGenMulvsRawMul")->GetEntries()<<" "<<ecutsmc->GetHistoCuts()->GetBinContent(3)<<endl;
}
else
{
tcutsdata->PrintCuts();
if(!managerdata||!ecutsdata||!tcutsdata)
return false;
+ if(managerdata->GetGenMulvsRawMulHistogram("hHistGenMulvsRawMul")->GetEntries()!=ecutsdata->GetHistoCuts()->GetBinContent(3))
+ cout<<"Please check DATA file possible problem with merging"<<" "<<anagerdata->GetGenMulvsRawMulHistogram("hHistGenMulvsRawMul")->GetEntries()<<" "<<ecutsdata->GetHistoCuts()->GetBinContent(3)<<endl;
+
}
return true;
}
{
Double_t lowedge=histo->GetBinLowEdge(ibin);
Float_t cut=dcacutxy->Eval(lowedge);
- TH1F* dcahist=(TH1F*)hman->GetDCAHistogram1D(name.Data(),lowedge,lowedge));
- Float_t inyield=dcahist->Integral(dcahist->GetXaxis()->FindBin(-1.0*cut),dcahist->GetXaxis()->FindBin(cut));
- cout<<"corr data "<<histo->GetBinContent(ibin)<<" "<<inyield<<" "<<dcahist->Integral()<<" "<<hnameout.Data()<<endl;
- cout<<"test dca "<<lowedge<<" "<<dcacutxy->Eval(lowedge)<<" "<<dcacutxy->Eval(histo->GetXaxis()->GetBinUpEdge(ibin))<<" "<<dcahist->GetBinLowEdge(dcahist->GetXaxis()->FindBin(-1.0*cut))<<" "<<dcahist->GetXaxis()->GetBinUpEdge(dcahist->GetXaxis()->FindBin(-1.0*cut))<<endl;
- histo->SetBinContent(ibin,inyield);
- histo->SetBinError(ibin,TMath::Sqrt(inyield));
+ TH1F* dcahist=(TH1F*)hman->GetDCAHistogram1D(name.Data(),lowedge,lowedge);
+// Float_t inyield=dcahist->Integral(dcahist->GetXaxis()->FindBin(-1.0*cut),dcahist->GetXaxis()->FindBin(cut));
+ Float_t testyield=0.0;
+ Float_t testerror=0.0;
+ for (int itest=dcahist->GetXaxis()->FindBin(-1.0*cut);itest<=dcahist->GetXaxis()->FindBin(cut);itest++)
+ {
+ testyield+=dcahist->GetBinContent(itest);
+ testerror+=dcahist->GetBinError(itest)*dcahist->GetBinError(itest);
+ }
+// cout<<"corr data "<<histo->GetBinContent(ibin)<<" "<<inyield<<" "<<dcahist->Integral()<<" "<<hnameout.Data()<<endl;
+// cout<<"test dca "<<lowedge<<" "<<dcacutxy->Eval(lowedge)<<" "<<dcacutxy->Eval(histo->GetXaxis()->GetBinUpEdge(ibin))<<" "<<dcahist->GetBinLowEdge(dcahist->GetXaxis()->FindBin(-1.0*cut))<<" "<<dcahist->GetXaxis()->GetBinUpEdge(dcahist->GetXaxis()->FindBin(-1.0*cut))<<endl;
+
+// cout<<testyield<<" "<<TMath::Sqrt(testerror)<<" error2 "<<inyield<<" "<<TMath::Sqrt(inyield)<<endl;
+
+ //histo->SetBinContent(ibin,inyield);
+ //histo->SetBinError(ibin,TMath::Sqrt(inyield));
+ histo->SetBinContent(ibin,testyield);
+ histo->SetBinError(ibin,TMath::Sqrt(testerror));
+
}
}
histo->Sumw2();
void GetPtHistFromPtDCAhisto(TString hnamein, TString hnameout, AliSpectraBothHistoManager* hman,TH1F** histo,TFormula* dcacutxy)
{
- Float_t min[3]={0.3,0.3,0.4};
- Float_t max[3]={1.5,1.2,2.2};
+ //Float_t min[3]={0.3,0.3,0.45};
+ //Float_t max[3]={1.5,1.2,2.2};
for(Int_t icharge=0;icharge<2;icharge++)
{
for(Int_t ipart=0;ipart<3;ipart++)
{
Int_t index=ipart+3*icharge;
- //TString hname=Form("hHistPtRecSigma%s%s",Particle[ipart].Data(),Sign[icharge].Data());
- printf("Getting %s",hnamein.Data());
+ Printf("Getting %s",hnamein.Data());
TString nameinfinal=Form("%s%s%s",hnamein.Data(),Particle[ipart].Data(),Sign[icharge].Data());
TString nameoutfinal=Form("%s%s%s",hnameout.Data(),Particle[ipart].Data(),Sign[icharge].Data());
histo[index]=GetOneHistFromPtDCAhisto(nameinfinal,nameoutfinal,hman,dcacutxy);
- /*histo[index] =(TH1F*)((TH1F*) hman->GetPtHistogram1D(Form("%s%s%s",hnamein.Data(),Particle[ipart].Data(),Sign[icharge].Data()),-1,-1))->Clone();
- histo[index]->SetName(Form("%s%s%s",hnameout.Data(),Particle[ipart].Data(),Sign[icharge].Data()));
- histo[index]->SetTitle(Form("%s%s%s",hnameout.Data(),Particle[ipart].Data(),Sign[icharge].Data()));
-
- if(dcacutxy)
- {
- for(int ibin=1;ibin<histo[index]->GetNbinsX();ibin++)
- {
- Double_t lowedge=histo[index]->GetBinLowEdge(ibin);
- Float_t cut=dcacutxy->Eval(lowedge);
- TH1F* dcahist=(TH1F*)hman->GetDCAHistogram1D(Form("%s%s%s",hnamein.Data(),Particle[ipart].Data(),Sign[icharge].Data()),lowedge,lowedge));
- Float_t inyield=dcahist->Integral(dcahist->GetXaxis()->FindBin(-1.0*cut),dcahist->GetXaxis()->FindBin(cut));
- cout<<"corr data "<<histo[index]->GetBinContent(ibin)<<" "<<inyield<<" "<<dcahist->Integral()<<endl;
- histo[index]->SetBinContent(ibin,inyield);
- histo[index]->SetBinError(ibin,TMath::Sqrt(inyield));
- }
- }*/
- CleanHisto(histo[index],min[ipart],max[ipart]);
- // histo[index]->Sumw2();
+ CleanHisto(histo[index],minRanges[ipart],maxRanges[ipart]);
}
}
}
}
if(contpid)
{
- if(contpid->GetBinContent(i)>0.2)
+ if(contpid->GetBinContent(i)>fMaxContaminationPIDMC)
{
h->SetBinContent(i,0);
h->SetBinError(i,0);
void DCACorrectionMarek(AliSpectraBothHistoManager* hman_data, AliSpectraBothHistoManager* hman_mc,TFormula* fun,TFile *fout,TH1F** hcon,TH1F** hconWD,TH1F** hconMat,TH1F** hprimary)
{
printf("\n\n-> DCA Correction");
- Double_t FitRange[2]={-2.0,2.0};
- Double_t CutRange[2]={-1.0,1.0};
- Double_t minptformaterial[6]={0.0,100.0,0.0,0.0,100.0,0.0};
- Double_t maxptformaterial[6]={0.0,-100.0,1.3,0.0,-100.0,0.0};
- Bool_t usefit[3]={true,false,true};
- Double_t maxbinforfit[6]={1.5,0,2.0,1.5,0,2.0};
+
+
Printf("\DCACorr");
- // TH1F *hcorrection[2];
- /* TCanvas *ccorrection=new TCanvas("DCAcorrection","DCAcorrection",700,500);
- TCanvas *cRatiocorrection=new TCanvas("DCARatiocorrection","DCARatiocorrection",700,500);
- cRatiocorrection->Divide(2,1);
- ccorrection->Divide(2,1);*/
TString sample[2]={"data","mc"};
ofstream debug("debugDCA.txt");
TList* listofdcafits=new TList();
for(Int_t isample=0;isample<2;isample++)
{
- /*hcorrection[isample]=(TH1F*)Spectra[index]->Clone();
- hcorrection[isample]->Reset("all");*/
for(Int_t ibin_data=7;ibin_data<40;ibin_data++)
{
Double_t lowedge=hcon[index]->GetBinLowEdge(ibin_data);
Double_t binwidth=hcon[index]->GetBinWidth(ibin_data);
- debug<<"NEW "<<Particle[ipart].Data()<<" "<<Sign[icharge].Data()<<" "<<lowedge<<endl;
- CutRange[1]=fun->Eval(lowedge);
- CutRange[0]=-1.0*CutRange[1];
+ debug<<"NEW "<<Particle[ipart].Data()<<" "<<Sign[icharge].Data()<<" "<<lowedge<<endl;
+ if(fun)
+ {
+ CutRange[1]=fun->Eval(lowedge);
+ CutRange[0]=-1.0*CutRange[1];
+ }
debug<<"cut "<<CutRange[1]<<" "<<CutRange[0]<<endl;
- Bool_t useMaterial=kFALSE;
- cout<<"try to fit "<<lowedge<<" "<<maxbinforfit[index]<<endl;
- if(lowedge>maxbinforfit[index])
- continue;
- if(lowedge>minptformaterial[index]&&lowedge<maxptformaterial[index])
- useMaterial=kTRUE;
-
+ Short_t fitsettings=DCAfitsettings(lowedge+0.5*binwidth,index);
+ debug<<"settings "<< fitsettings<<endl;
+ if(fitsettings==0)
+ continue;
+
TCanvas *cDCA=new TCanvas(Form("cDCA%d%s%s%sbin%d",index,sample[isample].Data(),Particle[ipart].Data(),Sign[icharge].Data(),ibin_data),Form("cDCA%d%s%s%sbin%d",index,sample[isample].Data(),Particle[ipart].Data(),Sign[icharge].Data(),ibin_data),1700,1500);
+ TLegend* Leg1 = new TLegend(0.6,0.6,0.85,0.85,"","NDC");
+ Leg1->SetFillStyle(kFALSE);
+ Leg1->SetLineColor(kWhite);
+ Leg1->SetBorderSize(0);
+
+
if(isample==0)
TH1F *hToFit =(TH1F*) ((TH1F*)hman_data->GetDCAHistogram1D(Form("hHistPtRecSigma%s%s",Particle[ipart].Data(),Sign[icharge].Data()),lowedge,lowedge))->Clone();
if(isample==1)
TH1F *hmc2=(TH1F*) ((TH1F*)hman_mc->GetDCAHistogram1D(Form("hHistPtRecSigmaSecondaryWeakDecay%s%s",Particle[ipart].Data(),Sign[icharge].Data()),lowedge,lowedge))->Clone();
TH1F *hmc3=(TH1F*) ((TH1F*)hman_mc->GetDCAHistogram1D(Form("hHistPtRecSigmaSecondaryMaterial%s%s",Particle[ipart].Data(),Sign[icharge].Data()),lowedge,lowedge))->Clone();
Double_t minentries=1;
- debug<<" Entries "<<isample<<" "<<hToFit->GetEntries()<<" "<<hmc1->GetEntries()<<" "<<hmc2->GetEntries()<<" "<<hmc3->GetEntries()<<endl;
- debug<<"2 Entries "<<isample<<" "<<hToFit->Integral(hmc1->GetXaxis()->FindBin(CutRange[0]),hmc1->GetXaxis()->FindBin(CutRange[1]))<<" "<<hmc1->Integral(hmc1->GetXaxis()->FindBin(CutRange[0]),hmc1->GetXaxis()->FindBin(CutRange[1]))<<" "<<hmc2->Integral(hmc1->GetXaxis()->FindBin(CutRange[0]),hmc1->GetXaxis()->FindBin(CutRange[1]))<<" "<<hmc3->Integral(hmc1->GetXaxis()->FindBin(CutRange[0]),hmc1->GetXaxis()->FindBin(CutRange[1]))<<endl;
- debug<< CutRange[0]<<" "<<CutRange[1]<<" "<<lowedge<<endl;
- if(hToFit->GetEntries()<=minentries || hmc1->GetEntries()<=minentries || hmc2->GetEntries()<=minentries || hmc3->GetEntries()<=minentries)
+ debug<<hToFit->GetEntries()<<" "<<hmc1->GetEntries()<<" "<<hmc2->GetEntries()<<" "<<hmc3->GetEntries()<<endl;
+ debug<<((fitsettings&0x1)&&hmc2->GetEntries()<=minentries)<<" "<<((fitsettings&0x2)&&hmc3->GetEntries()<=minentries)<<endl;
+ if(hToFit->GetEntries()<=minentries || hmc1->GetEntries()<=minentries || ((fitsettings&0x1)&&hmc2->GetEntries()<=minentries) || ((fitsettings&0x2)&&hmc3->GetEntries()<=minentries))
continue;
- //Data and MC can have different stat
hToFit->Sumw2();
hmc1->Sumw2();
hmc2->Sumw2();
hmc3->Sumw2();
Float_t corrforrebinning[4]={1.0,1.0,1.0,1.0};
+ Int_t binCutRange[]={hmc1->GetXaxis()->FindBin(CutRange[0]),hmc1->GetXaxis()->FindBin(CutRange[1])};
+
-
if(hmc3->GetNbinsX()>300)
{
- corrforrebinning[0]=hToFit->Integral(hmc1->GetXaxis()->FindBin(CutRange[0]),hmc1->GetXaxis()->FindBin(CutRange[1]));
- corrforrebinning[1]=hmc1->Integral(hmc1->GetXaxis()->FindBin(CutRange[0]),hmc1->GetXaxis()->FindBin(CutRange[1]));
- corrforrebinning[2]=hmc2->Integral(hmc1->GetXaxis()->FindBin(CutRange[0]),hmc1->GetXaxis()->FindBin(CutRange[1]));
- corrforrebinning[3]=hmc3->Integral(hmc1->GetXaxis()->FindBin(CutRange[0]),hmc1->GetXaxis()->FindBin(CutRange[1]));
+ corrforrebinning[0]=hToFit->Integral(binCutRange[0],binCutRange[1]);
+ corrforrebinning[1]=hmc1->Integral(binCutRange[0],binCutRange[1]);
+ corrforrebinning[2]=hmc2->Integral(binCutRange[0],binCutRange[1]);
+ corrforrebinning[3]=hmc3->Integral(binCutRange[0],binCutRange[1]);
hToFit->Rebin(30);
hmc1->Rebin(30);
hmc2->Rebin(30);
hmc3->Rebin(30);
- if(hToFit->Integral(hmc1->GetXaxis()->FindBin(CutRange[0]),hmc1->GetXaxis()->FindBin(CutRange[1]))>0.0)
- corrforrebinning[0]=corrforrebinning[0]/hToFit->Integral(hmc1->GetXaxis()->FindBin(CutRange[0]),hmc1->GetXaxis()->FindBin(CutRange[1]));
+
+ binCutRange[0]=hmc1->GetXaxis()->FindBin(CutRange[0]);
+ binCutRange[1]=hmc1->GetXaxis()->FindBin(CutRange[1]);
+
+
+ //after rebbing we lose resolution of the dca this correction also us to do obtain inside used dca
+
+ if(hToFit->Integral(binCutRange[0],binCutRange[1])>0.0)
+ corrforrebinning[0]=corrforrebinning[0]/hToFit->Integral(binCutRange[0],binCutRange[1]);
else
corrforrebinning[0]=1.0;
- if(hmc1->Integral(hmc1->GetXaxis()->FindBin(CutRange[0]),hmc1->GetXaxis()->FindBin(CutRange[1]))>0.0)
- corrforrebinning[1]=corrforrebinning[1]/hmc1->Integral(hmc1->GetXaxis()->FindBin(CutRange[0]),hmc1->GetXaxis()->FindBin(CutRange[1]));
+ if(hmc1->Integral(binCutRange[0],binCutRange[1])>0.0)
+ corrforrebinning[1]=corrforrebinning[1]/hmc1->Integral(binCutRange[0],binCutRange[1]);
else
corrforrebinning[1]=1.0;
- if(hmc2->Integral(hmc1->GetXaxis()->FindBin(CutRange[0]),hmc1->GetXaxis()->FindBin(CutRange[1]))>0.0)
- corrforrebinning[2]=corrforrebinning[2]/hmc2->Integral(hmc1->GetXaxis()->FindBin(CutRange[0]),hmc1->GetXaxis()->FindBin(CutRange[1]));
+ if(hmc2->Integral(binCutRange[0],binCutRange[1])>0.0)
+ corrforrebinning[2]=corrforrebinning[2]/hmc2->Integral(binCutRange[0],binCutRange[1]);
else
corrforrebinning[2]=1.0;
- if(hmc3->Integral(hmc1->GetXaxis()->FindBin(CutRange[0]),hmc1->GetXaxis()->FindBin(CutRange[1]))>0.0)
- corrforrebinning[3]=corrforrebinning[3]/hmc3->Integral(hmc1->GetXaxis()->FindBin(CutRange[0]),hmc1->GetXaxis()->FindBin(CutRange[1]));
+ if(hmc3->Integral(binCutRange[0],binCutRange[1])>0.0)
+ corrforrebinning[3]=corrforrebinning[3]/hmc3->Integral(binCutRange[0],binCutRange[1]);
else
corrforrebinning[3]=1.0;
- debug<<" cor bin "<<corrforrebinning[0]<<" "<<corrforrebinning[1]<<" "<<corrforrebinning[2]<<" "<<corrforrebinning[3]<<endl;
}
gPad->SetLogy();
TObjArray *mc=0x0;
- if(useMaterial)
+ Int_t Npar=3;
+ if(fitsettings==3)
mc = new TObjArray(3); // MC histograms are put in this array
else
+ {
mc = new TObjArray(2);
+ Npar=2;
+ }
+
mc->Add(hmc1);
- mc->Add(hmc2);
- if(useMaterial)
+ if(fitsettings&0x1)
+ mc->Add(hmc2);
+ if(fitsettings&0x2)
mc->Add(hmc3);
TFractionFitter* fit = new TFractionFitter(hToFit,mc); // initialise
fit->Constrain(0,0.0,1.0); // constrain fraction 1 to be between 0 and 1
- fit->Constrain(1,0.0,1.0); // constrain fraction 1 to be between 0 and 1
- if(useMaterial)
- fit->Constrain(2,0.0,1.0); // constrain fraction 1 to be between 0 and 1
- fit->SetRangeX(hToFit->GetXaxis()->FindBin(FitRange[0]),hToFit->GetXaxis()->FindBin(FitRange[1]));
- hToFit->GetXaxis()->SetRange(hToFit->GetXaxis()->FindBin(FitRange[0]),hToFit->GetXaxis()->FindBin(FitRange[1]));
+ Double_t defaultStep = 0.01;
+ for (int iparm = 0; iparm < Npar; ++iparm)
+ {
+ TString name("frac"); name += iparm;
+ if(iparm==0)
+ fit->GetFitter()->SetParameter(iparm,name.Data(),0.85,0.01,0.0,1.0);
+ else if (Npar==2)
+ fit->GetFitter()->SetParameter(iparm,name.Data(),0.15,0.01,0.0,1.0);
+ else
+ fit->GetFitter()->SetParameter(iparm,name.Data(),0.075,0.01,0.0,1.0);
+
+ }
+ if(fitsettings&0x1)
+ fit->Constrain(1,0.0,1.0); // constrain fraction 1 to be between 0 and 1
+ if(fitsettings&0x2)
+ fit->Constrain(1+(fitsettings&0x1),0.0,1.0); // constrain fraction 1 to be between 0 and 1
+
+ Int_t binFitRange[]={hmc1->GetXaxis()->FindBin(FitRange[0]),hmc1->GetXaxis()->FindBin(FitRange[1])};
+ fit->SetRangeX(binFitRange[0],binFitRange[1]);
+ hToFit->GetXaxis()->SetRange(binFitRange[0],binFitRange[1]);
hToFit->SetTitle(Form("DCA distr - %s %s %s %lf",sample[isample].Data(),Particle[ipart].Data(),Sign[icharge].Data(),lowedge));
Int_t status = fit->Fit(); // perform the fit
cout << "fit status: " << status << endl;
debug<<"fit status: " << status << endl;
- if (status == 0 && usefit[ipart])
- { // check on fit status
+ if (status == 0)
+ {
+ Double_t v1=0.0,v2=0.0,v3=0.0;
+ Double_t ev1=0.0,ev2=0.0,ev3=0.0;
+ Double_t cov=0.0;
+
+ // check on fit status
TH1F* result = (TH1F*) fit->GetPlot();
TH1F* PrimMCPred=(TH1F*)fit->GetMCPrediction(0);
- TH1F* secStMCPred=(TH1F*)fit->GetMCPrediction(1);
+
+ TH1F* secStMCPred=0X0;
TH1F* secMCPred=0x0;
- if(useMaterial)
- secMCPred=(TH1F*)fit->GetMCPrediction(2);
- Double_t v1=0,v2=0,v3=0;
- Double_t ev1=0,ev2=0,ev3=0;
- Double_t cov=0.0;
//first method, use directly the fit result
fit->GetResult(0,v1,ev1);
- fit->GetResult(1,v2,ev2);
- if(useMaterial)
+
+ if(fitsettings&0x1)
{
- fit->GetResult(2,v3,ev3);
- fit->GetFitter()->GetCovarianceMatrixElement(1,2);
+ fit->GetResult(1,v2,ev2);
+ secStMCPred=(TH1F*)fit->GetMCPrediction(1);
+
+ }
+ if(fitsettings&0x2)
+ {
+ fit->GetResult(1+(fitsettings&0x1),v3,ev3);
+ secMCPred=(TH1F*)fit->GetMCPrediction(1+(fitsettings&0x1));
+ if(fitsettings&0x1)
+ cov=fit->GetFitter()->GetCovarianceMatrixElement(1,2);
+
+
}
- debug<<v1<<" "<<ev1<<" "<<v2<<" "<<ev2<<" "<<v3<<" "<<ev3<<" "<<endl;
+
-
-
- Float_t normalizationdata=hToFit->Integral(hToFit->GetXaxis()->FindBin(CutRange[0]),hToFit->GetXaxis()->FindBin(CutRange[1]))/hToFit->Integral(hToFit->GetXaxis()->FindBin(FitRange[0]),hToFit->GetXaxis()->FindBin(FitRange[1]));
+ debug<<v1<<" "<<ev1<<" "<<v2<<" "<<ev2<<" "<<v3<<" "<<ev3<<" "<<" "<<cov<<endl;
+
+ // becuase dca cut range is not a fit range the results from TFractionFitter should be rescale
+ // This can be done in two ways or use input histograms or output histograms
+ // The difference between those two methods should be on the level of statistical error
+ // I use output histograms
- Float_t normalizationmc1=hmc1->Integral(hmc1->GetXaxis()->FindBin(CutRange[0]),hmc1->GetXaxis()->FindBin(CutRange[1]))/hmc1->Integral(hmc1->GetXaxis()->FindBin(FitRange[0]),hmc1->GetXaxis()->FindBin(FitRange[1]))/normalizationdata;
- Float_t normalizationmc2=hmc2->Integral(hmc2->GetXaxis()->FindBin(CutRange[0]),hmc2->GetXaxis()->FindBin(CutRange[1]))/hmc2->Integral(hmc1->GetXaxis()->FindBin(FitRange[0]),hmc2->GetXaxis()->FindBin(FitRange[1]))/normalizationdata;
- Float_t normalizationmc3=hmc3->Integral(hmc3->GetXaxis()->FindBin(CutRange[0]),hmc3->GetXaxis()->FindBin(CutRange[1]))/hmc3->Integral(hmc1->GetXaxis()->FindBin(FitRange[0]),hmc3->GetXaxis()->FindBin(FitRange[1]))/normalizationdata;
+ // Method 1 input histo
+
+ Float_t normalizationdata=hToFit->Integral(hToFit->GetXaxis()->FindBin(CutRange[0]),hToFit->GetXaxis()->FindBin(CutRange[1]))/hToFit->Integral(binFitRange[0],binFitRange[1]);
+ normalizationdata*=corrforrebinning[0];
+
+ Float_t normalizationmc1=(hmc1->Integral(binCutRange[0],binCutRange[1])/hmc1->Integral(binFitRange[0],binFitRange[1]))/normalizationdata;
+ Float_t normalizationmc2=0.0;
+ if(fitsettings&0x1)
+ normalizationmc2=(hmc2->Integral(binCutRange[0],binCutRange[1])/hmc2->Integral(binFitRange[0],binFitRange[1]))/normalizationdata;
+ Float_t normalizationmc3=0.0;
+ if(fitsettings&0x2)
+ normalizationmc3=(hmc3->Integral(binCutRange[0],binCutRange[1])/hmc3->Integral(binFitRange[0],binFitRange[1]))/normalizationdata;
+
+ normalizationmc1*=corrforrebinning[1];
+ normalizationmc2*=corrforrebinning[2];
+ normalizationmc3*=corrforrebinning[3];
+
debug<<"After Nor"<<endl;
debug<<v1*normalizationmc1<<" "<<ev1*normalizationmc1<<" "<<v2*normalizationmc2<<" "<<ev2*normalizationmc2<<" "<<v3*normalizationmc3<<" "<<ev3*normalizationmc3<<" "<<endl;
debug<<1.0-v1*normalizationmc1<<" "<<ev1*normalizationmc1<<" "<<v2*normalizationmc2+v3*normalizationmc3<<" "<<TMath::Sqrt(ev2*ev2*normalizationmc2*normalizationmc2+ev3*ev3*normalizationmc3*normalizationmc3+cov*normalizationmc3*normalizationmc2)<<endl;
- Float_t normalizationdata1=result->Integral(result->GetXaxis()->FindBin(CutRange[0]),result->GetXaxis()->FindBin(CutRange[1]))/result->Integral(result->GetXaxis()->FindBin(FitRange[0]),result->GetXaxis()->FindBin(FitRange[1]));
+ debug<<"addtional info"<<endl;
+
+
+ //Method 2 output histo
+
+ Float_t normalizationdata1=result->Integral(binCutRange[0],binCutRange[1])/result->Integral(binFitRange[0],binFitRange[1]);
+ // if the cut range is bigger the fit range we should calculate the normalization factor for data using the data histogram
+ // because result histogram has entries only in fits range
+ if(FitRange[0]>CutRange[0]||FitRange[1]<CutRange[1])
+ normalizationdata1=normalizationdata;
+
normalizationdata1*=corrforrebinning[0];
- Float_t normalizationmc11=PrimMCPred->Integral(PrimMCPred->GetXaxis()->FindBin(CutRange[0]),PrimMCPred->GetXaxis()->FindBin(CutRange[1]))/PrimMCPred->Integral(PrimMCPred->GetXaxis()->FindBin(FitRange[0]),PrimMCPred->GetXaxis()->FindBin(FitRange[1]))/normalizationdata1;
- Float_t normalizationmc21=secStMCPred->Integral(secStMCPred->GetXaxis()->FindBin(CutRange[0]),secStMCPred->GetXaxis()->FindBin(CutRange[1]))/secStMCPred->Integral(secStMCPred->GetXaxis()->FindBin(FitRange[0]),secStMCPred->GetXaxis()->FindBin(FitRange[1]))/normalizationdata1;
- Float_t normalizationmc31=0;
- if(useMaterial)
- normalizationmc31=secMCPred->Integral(secMCPred->GetXaxis()->FindBin(CutRange[0]),secMCPred->GetXaxis()->FindBin(CutRange[1]))/secMCPred->Integral(secMCPred->GetXaxis()->FindBin(FitRange[0]),secMCPred->GetXaxis()->FindBin(FitRange[1]))/normalizationdata1;
+ Float_t normalizationmc11=(PrimMCPred->Integral(binCutRange[0],binCutRange[1])/PrimMCPred->Integral(binFitRange[0],binFitRange[1]))/normalizationdata1;
+ Float_t normalizationmc21=0.0;
+ if(fitsettings&0x1)
+ normalizationmc21=(secStMCPred->Integral(binCutRange[0],binCutRange[1])/secStMCPred->Integral(binFitRange[0],binFitRange[1]))/normalizationdata1;
+ Float_t normalizationmc31=0.0;
+ if(fitsettings&0x2)
+ normalizationmc31=(secMCPred->Integral(binCutRange[0],binCutRange[1])/secMCPred->Integral(binFitRange[0],binFitRange[1]))/normalizationdata1;
normalizationmc11*=corrforrebinning[1];
normalizationmc21*=corrforrebinning[2];
normalizationmc31*=corrforrebinning[3];
- debug<<"After Nor 2"<<endl;
+ debug<<"After Nor 2"<<endl;
debug<<v1*normalizationmc11<<" "<<ev1*normalizationmc11<<" "<<v2*normalizationmc21<<" "<<ev2*normalizationmc21<<" "<<v3*normalizationmc31<<" "<<ev3*normalizationmc31<<endl;
debug<<1.0-v1*normalizationmc11<<" "<<ev1*normalizationmc11<<" "<<v2*normalizationmc21+v3*normalizationmc31<<" "<<TMath::Sqrt(ev2*ev2*normalizationmc21*normalizationmc21+ev3*ev3*normalizationmc31*normalizationmc31+cov*normalizationmc31*normalizationmc21)<<endl;
- debug<<CutRange[0]<<" "<<CutRange[1]<<endl;
- debug<<" Entries "<<isample<<" "<<hToFit->GetEntries()<<" "<<hmc1->GetEntries()<<" "<<hmc2->GetEntries()<<" "<<hmc3->GetEntries()<<endl;
- debug<<"2 Entries "<<isample<<" "<<hToFit->Integral(result->GetXaxis()->FindBin(CutRange[0]),result->GetXaxis()->FindBin(CutRange[1]))<<" "<<hmc1->Integral(result->GetXaxis()->FindBin(CutRange[0]),result->GetXaxis()->FindBin(CutRange[1]))<<" "<<hmc2->Integral(result->GetXaxis()->FindBin(CutRange[0]),result->GetXaxis()->FindBin(CutRange[1]))<<" "<<hmc3->Integral(result->GetXaxis()->FindBin(CutRange[0]),result->GetXaxis()->FindBin(CutRange[1]))<<endl;
-
hconWD[index+6*isample]->SetBinContent(ibin_data,v2*normalizationmc21);
hconWD[index+6*isample]->SetBinError(ibin_data,ev2*normalizationmc21);
hconMat[index+6*isample]->SetBinContent(ibin_data,v3*normalizationmc31);
hconMat[index+6*isample]->SetBinError(ibin_data,ev3*normalizationmc31);
hprimary[index+6*isample]->SetBinContent(ibin_data,v1*normalizationmc11);
- hprimary[index+6*isample]->SetBinContent(ibin_data,v1*normalizationmc11);
- if(useMaterial)
- {
- hcon[index+6*isample]->SetBinContent(ibin_data,v2*normalizationmc21+v3*normalizationmc31);
- hcon[index+6*isample]->SetBinError(ibin_data,TMath::Sqrt(ev2*ev2*normalizationmc21*normalizationmc21+ev3*ev3*normalizationmc31*normalizationmc31+cov*normalizationmc31*normalizationmc21));
- }
- else
- {
- hcon[index+6*isample]->SetBinContent(ibin_data,v2*normalizationmc21);
- hcon[index+6*isample]->SetBinError(ibin_data,ev2*normalizationmc21);
- }
+ hprimary[index+6*isample]->SetBinError(ibin_data,ev1*normalizationmc11);
+ hcon[index+6*isample]->SetBinContent(ibin_data,v2*normalizationmc21+v3*normalizationmc31);
+ hcon[index+6*isample]->SetBinError(ibin_data,TMath::Sqrt(ev2*ev2*normalizationmc21*normalizationmc21+ev3*ev3*normalizationmc31*normalizationmc31+cov*normalizationmc31*normalizationmc21));
result->Scale(1.0/result->Integral(result->GetXaxis()->FindBin(FitRange[0]),result->GetXaxis()->FindBin(FitRange[1])));
hToFit->Scale(1.0/hToFit->Integral(hToFit->GetXaxis()->FindBin(FitRange[0]),hToFit->GetXaxis()->FindBin(FitRange[1])));
PrimMCPred->Scale(v1/PrimMCPred->Integral(PrimMCPred->GetXaxis()->FindBin(FitRange[0]),PrimMCPred->GetXaxis()->FindBin(FitRange[1])));
- secStMCPred->Scale(v2/secStMCPred->Integral(secStMCPred->GetXaxis()->FindBin(FitRange[0]),secStMCPred->GetXaxis()->FindBin(FitRange[1])));
- if(useMaterial)
- secMCPred->Scale(v3/secMCPred->Integral(secMCPred->GetXaxis()->FindBin(FitRange[0]),secMCPred->GetXaxis()->FindBin(FitRange[1])));
-
- result->SetLineColor(kBlack);
- PrimMCPred->SetLineColor(kGreen+2);
- secStMCPred->SetLineColor(kRed);
-
+
hToFit->SetMinimum(0.0001);
hToFit->DrawClone("E1x0");
result->SetTitle("Fit result");
+ result->SetLineColor(kBlack);
+ Leg1->AddEntry(result,"result","lp");
result->DrawClone("lhistsame");
+
+ PrimMCPred->SetLineColor(kGreen+2);
+ PrimMCPred->SetLineStyle(2);
+ PrimMCPred->SetLineWidth(3.0);
+ Leg1->AddEntry(PrimMCPred,"Prmi.","l");
PrimMCPred->DrawClone("lhistsame");
- secStMCPred->DrawClone("lhistsame");
- if(useMaterial)
+ if(fitsettings&0x1)
{
- secMCPred->SetLineColor(kBlue);
- secMCPred->DrawClone("lhistsame");
+
+ secStMCPred->Scale(v2/secStMCPred->Integral(secStMCPred->GetXaxis()->FindBin(FitRange[0]),secStMCPred->GetXaxis()->FindBin(FitRange[1])));
+ secStMCPred->SetLineColor(kRed);
+ secStMCPred->SetLineWidth(3.0);
+
+ secStMCPred->SetLineStyle(3);
+ Leg1->AddEntry(secStMCPred,"Sec.WD","l");
+ secStMCPred->DrawClone("lhistsame");
+
}
- }
+ if(fitsettings&0x2)
+ {
+
+ secMCPred->Scale(v3/secMCPred->Integral(secMCPred->GetXaxis()->FindBin(FitRange[0]),secMCPred->GetXaxis()->FindBin(FitRange[1])));
+ secMCPred->SetLineColor(kBlue);
+ secMCPred->SetLineWidth(3.0);
+
+ secMCPred->SetLineStyle(4);
+ Leg1->AddEntry(secMCPred,"Sec.Mat","l");
+ secMCPred->DrawClone("lhistsame");
+
+ }
+ }
else
{
hconWD[index+6*isample]->SetBinContent(ibin_data,0.0);
hprimary[index+6*isample]->SetBinContent(ibin_data,1.0);
hprimary[index+6*isample]->SetBinError(ibin_data,0.0);
}
+ Leg1->Draw();
listofdcafits->Add(cDCA);
//cDCA->Write();
void RecomputeErrors(TH1* h)
{
for (int i=0; i<=h->GetXaxis()->GetNbins(); i++)
+ {
+ cout<<h->GetBinContent(i)<<" "<<h->GetBinError(i)<<" error "<<TMath::Sqrt(h->GetBinContent(i))<<endl;
h->SetBinError(i,TMath::Sqrt(h->GetBinContent(i)));
+ }
h->Sumw2();
}
void SetBintoOne(TH1* h)
gGFCorrectionKaonMinus->SetTitle("gGFCorrectionKaonMinus");
TString fnameGeanFlukaK="GFCorrection/correctionForCrossSection.321.root";
TFile *fGeanFlukaK= new TFile(fnameGeanFlukaK.Data());
- if (!fGeanFlukaK)
- return;
+ if (!fGeanFlukaK)
+ {
+ fnameGeanFlukaK="$ALICE_ROOT/PWGLF/SPECTRA/PiKaPr/TestAOD/GFCorrection/correctionForCrossSection.321.root";
+ fGeanFlukaK= new TFile(fnameGeanFlukaK.Data());
+ if (!fGeanFlukaK)
+ return;
+ }
TH1F *hGeantFlukaKPos=(TH1F*)fGeanFlukaK->Get("gHistCorrectionForCrossSectionParticles");
TH1F *hGeantFlukaKNeg=(TH1F*)fGeanFlukaK->Get("gHistCorrectionForCrossSectionAntiParticles");
//getting GF func for Kaons with TOF
const Int_t kNCharge=2;
Int_t kPos=0;
Int_t kNeg=1;
- TFile* fGFProtons = new TFile ("GFCorrection/correctionForCrossSection.root");
+ TString fnameGFProtons= "GFCorrection/correctionForCrossSection.root";
+ TFile* fGFProtons = new TFile (fnameGFProtons.Data());
+ if (!fGFProtons)
+ {
+ fnameGFProtons=="$ALICE_ROOT/PWGLF/SPECTRA/PiKaPr/TestAOD/GFCorrection/correctionForCrossSection.root";
+ TFile* fGFProtons = new TFile (fnameGFProtons.Data());
+ if (!fGFProtons)
+ return;
+ }
+
+
+
TH2D * hCorrFluka[kNCharge];
TH2D * hCorrFluka[2];
hCorrFluka[kPos] = (TH2D*)fGFProtons->Get("gHistCorrectionForCrossSectionProtons");
{
if(TOFMatchingScalling[0]<0.0&&TOFMatchingScalling[1]<0.0)
{
- TH1F *hMatcEffPos_data=(TH1F*)tcutsdata->GetHistoNMatchedPos();
- hMatcEffPos_data->Divide((TH1F*)tcutsdata->GetHistoNSelectedPos());
+ TH1F *hMatcEffPos_data=(TH1F*)tcutsdata->GetHistoNMatchedPos();
+ hMatcEffPos_data->Sumw2();
+ //hMatcEffPos_data->Divide((TH1F*)tcutsdata->GetHistoNSelectedPos());
+ hMatcEffPos_data->Divide(hMatcEffPos_data,(TH1F*)tcutsdata->GetHistoNSelectedPos(),1,1,"B");
hMatcEffPos_data->SetTitle("Matching Eff Pos - data");
TH1F *hMatcEffNeg_data=(TH1F*)tcutsdata->GetHistoNMatchedNeg();
- hMatcEffNeg_data->Divide((TH1F*)tcutsdata->GetHistoNSelectedNeg());
+ hMatcEffNeg_data->Sumw2();
+ //hMatcEffNeg_data->Divide((TH1F*)tcutsdata->GetHistoNSelectedNeg());
+ hMatcEffNeg_data->Divide(hMatcEffNeg_data,(TH1F*)tcutsdata->GetHistoNSelectedNeg(),1,1,"B");
hMatcEffNeg_data->SetTitle("Matching Eff Neg - data");
TH1F *hMatcEffPos_mc=(TH1F*)tcutsmc->GetHistoNMatchedPos();
- hMatcEffPos_mc->Divide((TH1F*)tcutsmc->GetHistoNSelectedPos());
+ hMatcEffPos_mc->Sumw2();
+ //hMatcEffPos_mc->Divide((TH1F*)tcutsmc->GetHistoNSelectedPos());
+ hMatcEffPos_mc->Divide(hMatcEffPos_mc,(TH1F*)tcutsmc->GetHistoNSelectedPos(),1,1,"B");
hMatcEffPos_mc->SetTitle("Matching Eff Pos - mc");
TH1F *hMatcEffNeg_mc=(TH1F*)tcutsmc->GetHistoNMatchedNeg();
- hMatcEffNeg_mc->Divide((TH1F*)tcutsmc->GetHistoNSelectedNeg());
+ hMatcEffNeg_mc->Sumw2();
+ //hMatcEffNeg_mc->Divide((TH1F*)tcutsmc->GetHistoNSelectedNeg());
+ hMatcEffNeg_mc->Divide(hMatcEffNeg_mc,(TH1F*)tcutsmc->GetHistoNSelectedNeg(),1,1,"B");
hMatcEffNeg_mc->SetTitle("Matching Eff Neg - mc");
TOFMatchingScalling[0]=pol0MatchPos_data->GetParameter(0);
TOFMatchingScalling[1]=pol0MatchNeg_data->GetParameter(0);
+
}
//Correction spectra for matching efficiency
//For the moment I'm using the inclusive correction
Float_t maxy=eta2y(pt,masstmp,0.8);
Float_t conver=maxy*(TMath::Sqrt(1-masstmp*masstmp/(mt2*TMath::CosH(maxy)*TMath::CosH(maxy)))+TMath::Sqrt(1-masstmp*masstmp/(mt2*TMath::CosH(0.0)*TMath::CosH(0.0))));
conver=conver/1.6;
- cout<<maxy<<" "<<conver<<" "<<masstmp<<""<<spectra[i]->GetName()<<endl;
+ //cout<<maxy<<" "<<conver<<" "<<masstmp<<""<<spectra[i]->GetName()<<endl;
Float_t bincontent=allch->GetBinContent(j);
Float_t binerror=allch->GetBinError(j);
bincontent+=conver*value;
}
}
+
+Short_t DCAfitsettings (Float_t pt, Int_t type)
+{
+ Short_t value=0x0;
+ if (pt<maxptformaterial[type]&&pt>minptformaterial[type])
+ value=value+2;
+ if (pt<maxptforWD[type]&&pt>minptforWD[type])
+ value=value+1;
+ return value;
+
+}
+
+Float_t Normaliztionwithbin0integrals(UInt_t options)
+{
+
+ TH1F* bin0mcRec=(TH1F*)ecutsmc->GetHistoVtxGenerated()->Clone("Bin0_rec");
+ TH1F* bin0mcMC=(TH1F*)ecutsmc->GetHistoVtxGenerated()->Clone("Bin0_MC");
+
+ TH1F* vertexmc=ecutsmc->GetHistoVtxAftSelwithoutZvertexCut();
+ TH1F* vertexmcMCz=ecutsmc->GetHistoVtxAftSelwithoutZvertexCutusingMCz();
+ TH1F* vertexdata=ecutsdata->GetHistoVtxAftSelwithoutZvertexCut();
+
+ TH1I* histodata=ecutsdata->GetHistoCuts();
+ TH1I* histomc=ecutsmc->GetHistoCuts();
+
+ Float_t dataevents=(Float_t)histodata->GetBinContent(3);
+ //cout<<histodata->GetBinContent(2)<<endl;
+ Float_t databin0events=((Float_t)histodata->GetBinContent(2))-((Float_t)histodata->GetBinContent(4));
+
+ bin0mcRec->Sumw2();
+ bin0mcMC->Sumw2();
+
+ bin0mcRec->Add(vertexmc,-1);
+ bin0mcMC->Add(vertexmcMCz,-1);
+
+ bin0mcRec->Divide(vertexmc);
+ bin0mcMC->Divide(vertexmcMCz);
+
+ bin0mcRec->Multiply(vertexdata);
+ bin0mcMC->Multiply(vertexdata);
+
+ Float_t bin0mcRecN=0.0;
+ Float_t bin0mcMCN=0.0;
+
+ for (int i=0;i<=bin0mcRec->GetXaxis()->GetNbins();i++)
+ {
+ bin0mcRecN+=bin0mcRec->GetBinContent(i);
+ bin0mcMCN+=bin0mcMC->GetBinContent(i);
+
+ }
+ bin0mcRec->Scale(databin0events/bin0mcRecN);
+ bin0mcMC->Scale(databin0events/bin0mcMCN);
+
+ Int_t binmin=bin0mcRec->GetXaxis()->FindBin(-10);
+ Int_t binmax=bin0mcRec->GetXaxis()->FindBin(10)-1;
+ cout<< bin0mcRec->GetXaxis()->GetBinLowEdge(binmin)<<" "<<bin0mcRec->GetXaxis()->GetBinUpEdge(binmax)<<endl;
+ cout<<bin0mcRecN<<" "<<bin0mcMCN<<" "<<databin0events<<endl;
+ cout<<dataevents<<" normalization "<<dataevents+bin0mcRec->Integral(binmin,binmax)<<" "<<dataevents+bin0mcMC->Integral(binmin,binmax)<<endl;
+ cout<<histodata->GetBinContent(2)<<" "<<histodata->GetBinContent(4)<<endl;
+ if ((options&knormalizationwithbin0integralsdata)==knormalizationwithbin0integralsdata)
+ return dataevents+bin0mcRec->Integral(binmin,binmax);
+ else if ((options&kknormalizationwithbin0integralsMC)==knormalizationwithbin0integralsMC)
+ return dataevents+bin0mcMC->Integral(binmin,binmax) ;
+ else
+ return 1;
+}
+
+
+Bool_t ReadConfigFile(TString configfile)
+{
+ ifstream infile(configfile.Data());
+ if(infile.is_open()==false)
+ return false;
+ TString namesofSetting[35]={"CutRangeMin","CutRangeMax","FitRangeMin","FitRangeMax","MinMatPionPlus","MaxMatPionPlus","MinMatKaonPlus","MaxMatKaonPlus","MinMatProtonPlus","MaxMatProtonPlus","MinMatPionMinus","MaxMatPionMinus","MinMatKaonMinus","MaxMatKaonMinus","MinMatProtonMinus","MaxMatProtonMinus","MinWDPionPlus","MaxWDPionPlus","MinWDKaonPlus","MaxWDKaonPlus","MinWDProtonPlus","MaxWDProtonPlus","MinWDPionMinus","MaxWDPionMinus","MinWDKaonMinus","MaxWDKaonMinus","MinWDProtonMinus","MaxWDProtonMinus","MaxContaminationPIDMC","MinPions","MaxPions","MinKaons","MaxKaons","MinProtons","MaxProtons"};
+
+ char buffer[256];
+ while (infile.eof()==false)
+ {
+ buffer[0]='#';
+ while (buffer[0]=='#'&&infile.eof()==false)
+ infile.getline(buffer,256);
+ TString tmpstring(buffer);
+ cout<<buffer<<endl;
+ if(tmpstring.Contains(namesofSetting[0]))
+ CutRange[0]=(tmpstring.Remove(0,namesofSetting[0].Length()+1)).Atof();
+ else if (tmpstring.Contains(namesofSetting[1]))
+ CutRange[1]=(tmpstring.Remove(0,namesofSetting[1].Length()+1)).Atof();
+ else if (tmpstring.Contains(namesofSetting[2]))
+ FitRange[0]=(tmpstring.Remove(0,namesofSetting[2].Length()+1)).Atof();
+ else if (tmpstring.Contains(namesofSetting[3]))
+ FitRange[1]=(tmpstring.Remove(0,namesofSetting[3].Length()+1)).Atof();
+ else if (tmpstring.Contains(namesofSetting[4]))
+ minptformaterial[0]=(tmpstring.Remove(0,namesofSetting[4].Length()+1)).Atof();
+ else if (tmpstring.Contains(namesofSetting[5]))
+ maxptformaterial[0]=(tmpstring.Remove(0,namesofSetting[5].Length()+1)).Atof();
+ else if (tmpstring.Contains(namesofSetting[6]))
+ minptformaterial[1]=(tmpstring.Remove(0,namesofSetting[6].Length()+1)).Atof();
+ else if (tmpstring.Contains(namesofSetting[7]))
+ maxptformaterial[1]=(tmpstring.Remove(0,namesofSetting[7].Length()+1)).Atof();
+ else if (tmpstring.Contains(namesofSetting[8]))
+ minptformaterial[2]=(tmpstring.Remove(0,namesofSetting[8].Length()+1)).Atof();
+ else if (tmpstring.Contains(namesofSetting[9]))
+ maxptformaterial[2]=(tmpstring.Remove(0,namesofSetting[9].Length()+1)).Atof();
+ else if (tmpstring.Contains(namesofSetting[10]))
+ minptformaterial[3]=(tmpstring.Remove(0,namesofSetting[10].Length()+1)).Atof();
+ else if (tmpstring.Contains(namesofSetting[11]))
+ maxptformaterial[3]=(tmpstring.Remove(0,namesofSetting[11].Length()+1)).Atof();
+ else if (tmpstring.Contains(namesofSetting[12]))
+ minptformaterial[4]=(tmpstring.Remove(0,namesofSetting[12].Length()+1)).Atof();
+ else if (tmpstring.Contains(namesofSetting[13]))
+ maxptformaterial[4]=(tmpstring.Remove(0,namesofSetting[13].Length()+1)).Atof();
+ else if (tmpstring.Contains(namesofSetting[14]))
+ minptformaterial[5]=(tmpstring.Remove(0,namesofSetting[14].Length()+1)).Atof();
+ else if (tmpstring.Contains(namesofSetting[15]))
+ maxptformaterial[5]=(tmpstring.Remove(0,namesofSetting[15].Length()+1)).Atof();
+ else if (tmpstring.Contains(namesofSetting[16]))
+ minptforWD[0]=(tmpstring.Remove(0,namesofSetting[16].Length()+1)).Atof();
+ else if (tmpstring.Contains(namesofSetting[17]))
+ maxptforWD[0]=(tmpstring.Remove(0,namesofSetting[17].Length()+1)).Atof();
+ else if (tmpstring.Contains(namesofSetting[18]))
+ minptforWD[1]=(tmpstring.Remove(0,namesofSetting[18].Length()+1)).Atof();
+ else if (tmpstring.Contains(namesofSetting[19]))
+ maxptforWD[1]=(tmpstring.Remove(0,namesofSetting[19].Length()+1)).Atof();
+ else if (tmpstring.Contains(namesofSetting[20]))
+ minptforWD[2]=(tmpstring.Remove(0,namesofSetting[20].Length()+1)).Atof();
+ else if (tmpstring.Contains(namesofSetting[21]))
+ maxptforWD[2]=(tmpstring.Remove(0,namesofSetting[21].Length()+1)).Atof();
+ else if (tmpstring.Contains(namesofSetting[22]))
+ minptforWD[3]=(tmpstring.Remove(0,namesofSetting[22].Length()+1)).Atof();
+ else if (tmpstring.Contains(namesofSetting[23]))
+ maxptforWD[3]=(tmpstring.Remove(0,namesofSetting[23].Length()+1)).Atof();
+ else if (tmpstring.Contains(namesofSetting[24]))
+ minptforWD[4]=(tmpstring.Remove(0,namesofSetting[24].Length()+1)).Atof();
+ else if (tmpstring.Contains(namesofSetting[25]))
+ maxptforWD[4]=(tmpstring.Remove(0,namesofSetting[25].Length()+1)).Atof();
+ else if (tmpstring.Contains(namesofSetting[26]))
+ minptforWD[5]=(tmpstring.Remove(0,namesofSetting[26].Length()+1)).Atof();
+ else if (tmpstring.Contains(namesofSetting[27]))
+ maxptforWD[5]=(tmpstring.Remove(0,namesofSetting[27].Length()+1)).Atof();
+ else if (tmpstring.Contains(namesofSetting[28]))
+ fMaxContaminationPIDMC=(tmpstring.Remove(0,namesofSetting[28].Length()+1)).Atof();
+ else if (tmpstring.Contains(namesofSetting[29]))
+ minRanges[0]=(tmpstring.Remove(0,namesofSetting[29].Length()+1)).Atof();
+ else if (tmpstring.Contains(namesofSetting[30]))
+ maxRanges[0]=(tmpstring.Remove(0,namesofSetting[30].Length()+1)).Atof();
+ else if (tmpstring.Contains(namesofSetting[31]))
+ minRanges[1]=(tmpstring.Remove(0,namesofSetting[31].Length()+1)).Atof();
+ else if (tmpstring.Contains(namesofSetting[32]))
+ maxRanges[1]=(tmpstring.Remove(0,namesofSetting[32].Length()+1)).Atof();
+ else if (tmpstring.Contains(namesofSetting[33]))
+ minRanges[2]=(tmpstring.Remove(0,namesofSetting[33].Length()+1)).Atof();
+ else if (tmpstring.Contains(namesofSetting[34]))
+ maxRanges[2]=(tmpstring.Remove(0,namesofSetting[34].Length()+1)).Atof();
+ else
+ continue;
+
+
+// Double_t minRanges[3]={0.3,0.3,0.45};
+//Double_t maxRanges[3]={1.5,1.2,2.2};
+//Double_t fMaxContaminationPIDMC=0.2;
+
+
+
+ }
+ for(int i=0;i<6;i++)
+ cout<<minptformaterial[i]<<" "<<maxptformaterial[i]<<" "<<minptforWD[i]<<" "<<maxptforWD[i]<<endl;
+ cout<<FitRange[0]<<" "<<FitRange[1]<<" "<<CutRange[0]<<CutRange[1]<<endl;
+ if(FitRange[0]>=FitRange[1])
+ {
+ cout<<"A"<<endl;
+ return false;
+ }
+ if(CutRange[0]>=CutRange[1])
+ {
+ cout<<"B"<<endl;
+ return false;
+ }
+ for(int i=0;i<6;i++)
+ {
+ if((minptformaterial[i]>maxptformaterial[i]&&minptformaterial[i]>0.0)||minptformaterial[i]<0.0||maxptformaterial[i]<0.0)
+ {
+ cout<<"C"<<endl;
+ return false;
+ }
+ if((minptforWD[i]>maxptforWD[i]&&minptforWD[i]>0.0)||minptforWD[i]<0.0||maxptforWD[i]<0.0)
+ {
+ cout<<"D"<<endl;
+ return false;
+ }
+ }
+ for(int i=0;i<3;i++)
+ if(minRanges[i]>maxRanges[i])
+ return false;
+
+
+ return true;
+}
+
+void SubHistWithFullCorr(TH1F* h1, TH1F* h2, Float_t factor1=1.0, Float_t factor2=1.0)
+{
+ if(h1->GetNbinsX()!=h2->GetNbinsX())
+ return;
+ for (int i=0;i<=h1->GetNbinsX();i++)
+ {
+ Float_t tmpvalue=factor1*h1->GetBinContent(i)-factor2*h2->GetBinContent(i);
+ Float_t tmperror=TMath::Abs(factor1*factor1*h1->GetBinError(i)*h1->GetBinError(i)-factor2*factor2*h2->GetBinError(i)*h2->GetBinError(i));
+ h1->SetBinContent(i,tmpvalue);
+ h1->SetBinError(i,TMath::Sqrt(tmperror));
+ }
+
+}
+
+void TOFMatchingForNch(TH1* h)
+{
+ if(TOFMatchingScalling[0]>0.0&&TOFMatchingScalling[1]>0.0)
+ {
+ Float_t factor=0.5*TOFMatchingScalling[0]+0.5*TOFMatchingScalling[1];
+ for(Int_t ibin=1;ibin<h->GetNbinsX();ibin++)
+ {
+ Float_t ptspectra=h->GetBinCenter(ibin);
+ if(ptspectra<tcutsdata->GetPtTOFMatching())
+ continue;
+ h->SetBinContent(ibin,(h->GetBinContent(ibin)/factor));
+ }
+
+ }
+ else
+ return;
+
+
+}
git://git.uio.no / u / mrichter / AliRoot.git / blobdiff