#include <TH1F.h>
#include <TH2F.h>
#include <TH3F.h>
+#include <TProfile.h>
#include <TList.h>
#include <TLorentzVector.h>
#include <TClonesArray.h>
fLimitGenJetEta(kFALSE),
fAnalysisType(0),
fExternalWeight(1),
+ fh1Xsec(0x0),
+ fh1Trials(0x0),
fh1PtHard(0x0),
fh1PtHard_NoW(0x0),
fh1PtHard_Trials(0x0),
fLimitGenJetEta(kFALSE),
fAnalysisType(0),
fExternalWeight(1),
+ fh1Xsec(0x0),
+ fh1Trials(0x0),
fh1PtHard(0x0),
fh1PtHard_NoW(0x0),
fh1PtHard_Trials(0x0),
+Bool_t AliAnalysisTaskJetSpectrum::Notify()
+{
+ TTree *tree = AliAnalysisManager::GetAnalysisManager()->GetTree();
+ Double_t xsection = 0;
+ UInt_t ntrials = 0;
+ if(tree){
+ TFile *curfile = tree->GetCurrentFile();
+ if (!curfile) {
+ Error("Notify","No current file");
+ return kFALSE;
+ }
+ if(!fh1Xsec||!fh1Trials){
+ Printf("%s%d No Histogram fh1Xsec",(char*)__FILE__,__LINE__);
+ return kFALSE;
+ }
+
+ TString fileName(curfile->GetName());
+ if(fileName.Contains("AliESDs.root")){
+ fileName.ReplaceAll("AliESDs.root", "pyxsec.root");
+ }
+ else if(fileName.Contains("AliAOD.root")){
+ fileName.ReplaceAll("AliAOD.root", "pyxsec.root");
+ }
+ else if(fileName.Contains("galice.root")){
+ // for running with galice and kinematics alone...
+ fileName.ReplaceAll("galice.root", "pyxsec.root");
+ }
+ TFile *fxsec = TFile::Open(fileName.Data());
+ if(!fxsec){
+ Printf("%s:%d %s not found in the Input",(char*)__FILE__,__LINE__,fileName.Data());
+ // no a severe condition
+ return kTRUE;
+ }
+ TTree *xtree = (TTree*)fxsec->Get("Xsection");
+ if(!xtree){
+ Printf("%s:%d tree not found in the pyxsec.root",(char*)__FILE__,__LINE__);
+ }
+ xtree->SetBranchAddress("xsection",&xsection);
+ xtree->SetBranchAddress("ntrials",&ntrials);
+ xtree->GetEntry(0);
+ fh1Xsec->Fill("<#sigma>",xsection);
+ fh1Trials->Fill("#sum{ntrials}",ntrials);
+ }
+
+ return kTRUE;
+}
void AliAnalysisTaskJetSpectrum::UserCreateOutputObjects()
{
const Int_t nBinFrag = 25;
+ fh1Xsec = new TProfile("fh1Xsec","xsec from pyxsec.root",1,0,1);
+ fh1Xsec->GetXaxis()->SetBinLabel(1,"<#sigma>");
+
+ fh1Trials = new TH1F("fh1Trials","trials from pyxsec.root",1,0,1);
+ fh1Trials->GetXaxis()->SetBinLabel(1,"#sum{ntrials}");
+
fh1PtHard = new TH1F("fh1PtHard","PYTHIA Pt hard;p_{T,hard}",nBinPt,binLimitsPt);
fh1PtHard_NoW = new TH1F("fh1PtHard_NoW","PYTHIA Pt hard no weight;p_{T,hard}",nBinPt,binLimitsPt);
}
const Int_t saveLevel = 1; // large save level more histos
if(saveLevel>0){
+ fHistList->Add(fh1Xsec);
+ fHistList->Add(fh1Trials);
fHistList->Add(fh1PtHard);
fHistList->Add(fh1PtHard_NoW);
fHistList->Add(fh1PtHard_Trials);
if (h1){
// Printf("%s ",h1->GetName());
h1->Sumw2();
+ continue;
}
+ THnSparse *hn = dynamic_cast<THnSparse*>(fHistList->At(i));
+ if(hn)hn->Sumw2();
}
TH1::AddDirectory(oldStatus);
TFile::Open(fileNameRec);
TH2F* hist = (TH2F*) gFile->Get("unfolding/fRecSpectrum");
jetSpec->SetRecSpectrum(hist);
+ hist = (TH2F*) gFile->Get("unfolding/fGenSpectrum");
+ if(hist->GetEntries()>0)jetSpec->SetGenSpectrum(hist);
jetSpec->ApplyBayesianMethod(0.3, 20, 0, 0);
// last parameter = calculateErrors <- this method to calculate the errors takes a lot of time
jetSpec->GetRecSpectrum()->SetBinError(bine, binz, err + fhERecZRec->GetBinError(me, mz));
}
+
+ // for control again, but now from the rec file
+ // generated jets (true distribution)
+ TH2F *fhEGenZGen = (TH2F*)(tlist->FindObject("fh2EGenZGen"));
+ for (Int_t te=1; te<=fhEGenZGen->GetNbinsX(); te++)
+ for (Int_t tz=1; tz<=fhEGenZGen->GetNbinsY(); tz++)
+ {
+ Float_t ej = fhEGenZGen->GetXaxis()->GetBinCenter(te);
+ Float_t z = fhEGenZGen->GetYaxis()->GetBinCenter(tz);
+ Int_t bine = jetSpec->GetGenSpectrum()->GetXaxis()->FindBin(ej);
+ Int_t binz = jetSpec->GetGenSpectrum()->GetYaxis()->FindBin(z);
+ Float_t cont = jetSpec->GetGenSpectrum()->GetBinContent(bine,binz);
+ Float_t err = jetSpec->GetGenSpectrum()->GetBinError(bine,binz);
+ jetSpec->GetGenSpectrum()->SetBinContent(bine, binz, cont + fhEGenZGen->GetBinContent(te, tz));
+ jetSpec->GetGenSpectrum()->SetBinError(bine, binz, err + fhEGenZGen->GetBinError(te, tz));
+ }
+
+
file = TFile::Open("rec_pwg4spec.root", "RECREATE");
jetSpec->SaveHistograms();
file->Close();
}
-
-
-
void correct(){
// simple steering to correct a given distribution;
loadlibs();
- FillSpecFromFiles("pwg4spec_0000000-0010000.root","pwg4spec_15-50.root");
+ // rec and gen
+ // FillSpecFromFiles("pwg4spec_15-50_all.root","pwg4spec_allpt.root");
+ FillSpecFromFiles("pwg4spec_allpt.root","pwg4spec_allpt.root");
char name[100];
sprintf(name, "unfolded_pwg4spec.root");
//draw(name, "unfolding", 1);
}
+
+void mergeJetAnaOutput(){
+ // This is used to merge the analysis-output from different
+ // data samples/pt_hard bins
+ // in case the eventweigth was set to xsection/ntrials already this
+ // is not needed. Both methods only work in cse we do not mix different
+ // pt_hard bins, and do not have overlapping bins
+
+ const Int_t nBins = 2;
+ // LHC08q jetjet100: Mean = 1.42483e-03, RMS = 6.642e-05
+ // LHC08r jetjet50: Mean = 2.44068e-02, RMS = 1.144e-03
+ // LHC08v jetjet15-50: Mean = 2.168291 , RMS = 7.119e-02
+ const Float_t xsection[nBins] = {2.168291,2.44068e-02};
+ Float_t nTrials[nBins] = {0,0};
+ Float_t sf[nBins] = {0,0};
+
+ const char *cFile[nBins] = {"pwg4spec_15-50_all.root","pwg4spec_50_all.root"};
+
+
+ TList *lIn[nBins];
+ TFile *fIn[nBins];
+ for(int ib = 0;ib < nBins;++ib){
+ fIn[ib] = TFile::Open(cFile[ib]);
+ lIn[ib] = (TList*)fIn[ib]->Get("pwg4spec");
+ TH1* hTrials = (TH1F*)lIn[ib]->FindObject("fh1PtHard_Trials");
+ nTrials[ib] = hTrials->Integral();
+ sf[ib] = xsection[ib]/ nTrials[ib];
+ }
+
+ TFile *fOut = new TFile("pwg4spec_allpt.root","RECREATE");
+ TList *lOut = new TList();
+ lOut->SetName(lIn[0]->GetName());
+ // for the start scale all...
+ for(int ie = 0; ie < lIn[0]->GetEntries();++ie){
+ TH1 *h1Add = 0;
+ THnSparse *hnAdd = 0;
+ Printf("%d: %s",ie, lIn[0]->At(ie)->GetName());
+ for(int ib = 0;ib < nBins;++ib){
+ // dynamic cast does not work with cint
+ TObject *h = lIn[ib]->At(ie);
+ if(h->InheritsFrom("TH1")){
+ Printf("ib %d",ib);
+ TH1 *h1 = (TH1*)h;
+ if(ib==0){
+ h1Add = (TH1*)h1->Clone(h1->GetName());
+ h1Add->Scale(sf[ib]);
+ }
+ else{
+ h1Add->Add(h1,sf[ib]);
+ }
+ }
+ else if(h->InheritsFrom("THnSparse")){
+ Printf("ib %d",ib);
+ THnSparse *hn = (THnSparse*)h;
+ if(ib==0){
+ hnAdd = (THnSparse*)hn->Clone(hn->GetName());
+ hnAdd->Scale(sf[ib]);
+ }
+ else{
+ hnAdd->Add(hn,sf[ib]);
+ }
+ }
+
+
+ }// ib
+ if(h1Add)lOut->Add(h1Add);
+ else if(hnAdd)lOut->Add(hnAdd);
+ }
+ fOut->cd();
+ lOut->Write(lOut->GetName(),TObject::kSingleKey);
+ fOut->Close();
+
+
+
+}
git://git.uio.no / u / mrichter / AliRoot.git / commitdiff