--- /dev/null
+#include "AliCollisionNormalization.h"
+#include "AliPhysicsSelection.h"
+#include "AliLog.h"
+#include "TFile.h"
+#include "TCanvas.h"
+#include "AliGenPythiaEventHeader.h"
+#include "AliGenDPMjetEventHeader.h"
+#include "AliGenEventHeader.h"
+#include "AliMCEvent.h"
+
+ClassImp(AliCollisionNormalization)
+
+const char * AliCollisionNormalization::fProcLabel[] = {"SD","DD","ND", "Unknown"};
+
+AliCollisionNormalization::AliCollisionNormalization() :
+ TObject(),
+ fNbinsVz(0),
+ fMinVz(0),
+ fMaxVz(0),
+ fZRange(9.99),
+ fIsMC(0),
+ fReferenceXS(0),
+ fVerbose(0),
+ fHistVzData (0),
+ fHistProcTypes (0),
+ fHistStatBin0 (0)
+{
+
+ // ctor
+
+ fNbinsVz = 30;
+ fMinVz = -15;
+ fMaxVz = +15;
+
+ for(Int_t iproc = 0; iproc < kNProcs; iproc++){
+ fHistVzMCGen[iproc] = 0;
+ fHistVzMCRec[iproc] = 0;
+ fHistVzMCTrg[iproc] = 0;
+ }
+
+
+ BookAllHistos();
+}
+
+AliCollisionNormalization::AliCollisionNormalization(Int_t nbinz, Float_t minz, Float_t maxz):
+ TObject(),
+ fNbinsVz(0),
+ fMinVz(0),
+ fMaxVz(0),
+ fZRange(9.99),
+ fIsMC(0),
+ fReferenceXS(0),
+ fVerbose(0),
+ fHistVzData (0),
+ fHistProcTypes (0),
+ fHistStatBin0 (0)
+{
+
+ // ctor, allows setting binning
+ fNbinsVz = nbinz;
+ fMinVz = minz;
+ fMaxVz = maxz;
+
+ for(Int_t iproc = 0; iproc < kNProcs; iproc++){
+ fHistVzMCGen[iproc] = 0;
+ fHistVzMCRec[iproc] = 0;
+ fHistVzMCTrg[iproc] = 0;
+ }
+
+ BookAllHistos();
+}
+
+AliCollisionNormalization::AliCollisionNormalization(const char * dataFile, const char * dataListName,
+ const char * mcFile, const char * mcListName,
+ const char * eventStatFile) :
+ TObject(),
+ fNbinsVz(0),
+ fMinVz(0),
+ fMaxVz(0),
+ fZRange(9.99),
+ fIsMC(0),
+ fReferenceXS(0),
+ fVerbose(0),
+ fHistVzData (0),
+ fHistProcTypes (0),
+ fHistStatBin0 (0)
+{
+
+ // ctor, loads histograms from file
+ for(Int_t iproc = 0; iproc < kNProcs; iproc++){
+ fHistVzMCGen[iproc] = 0;
+ fHistVzMCRec[iproc] = 0;
+ fHistVzMCTrg[iproc] = 0;
+ }
+
+
+ TFile * fdata = new TFile (dataFile);
+ TFile * fmc = new TFile (mcFile );
+
+ TList * ldata = (TList*) fdata->Get(dataListName);
+ TList * lmc = (TList*) fmc ->Get(mcListName );
+
+ AliCollisionNormalization * cndata = (AliCollisionNormalization*) ldata->FindObject("AliCollisionNormalization");
+ AliCollisionNormalization * cnmc = (AliCollisionNormalization*) lmc ->FindObject("AliCollisionNormalization");
+
+
+ // Assign or book all histos
+ for(Int_t iproc = 0; iproc < kNProcs; iproc++){
+ fHistVzMCGen[iproc]= cnmc->GetVzMCGen(iproc);
+ fHistVzMCRec[iproc]= cnmc->GetVzMCRec(iproc);
+ fHistVzMCTrg[iproc]= cnmc->GetVzMCTrg(iproc);
+ }
+ fHistVzData = cndata->GetVzData();
+ fHistProcTypes = cnmc->GetHistProcTypes();
+
+ TFile * fstat = new TFile(eventStatFile);
+ fHistStatBin0 = (TH1F*) fstat->Get("fHistStatistics_Bin0");
+
+}
+
+
+AliCollisionNormalization::~AliCollisionNormalization(){
+
+ // dtor
+ for(Int_t iproc = 0; iproc < kNProcs; iproc++){
+ if(fHistVzMCGen[iproc]) { delete fHistVzMCGen[iproc] ; fHistVzMCGen[iproc] =0;}
+ if(fHistVzMCRec[iproc]) { delete fHistVzMCRec[iproc] ; fHistVzMCRec[iproc] =0;}
+ if(fHistVzMCTrg[iproc]) { delete fHistVzMCTrg[iproc] ; fHistVzMCTrg[iproc] =0;}
+ }
+
+ if(fHistVzData ) { delete fHistVzData ; fHistVzData =0;}
+ if(fHistStatBin0 ) { delete fHistStatBin0 ; fHistStatBin0 =0;}
+ if(fHistProcTypes ) { delete fHistProcTypes ; fHistProcTypes =0;}
+
+}
+
+void AliCollisionNormalization::BookAllHistos(){
+ // books all histos
+ // Book histos of vz distributions vs multiplicity
+ // if vzOnly == kTRUE, it returns a 1d histo with vz dist only
+
+ // Do not attach histos to the current directory
+ Bool_t oldStatus = TH1::AddDirectoryStatus();
+ TH1::AddDirectory(kFALSE);
+
+ for(Int_t iproc = 0; iproc < kNProcs; iproc++){
+ fHistVzMCGen [iproc] = (TH2F*) BookVzHisto(TString("fHistVzMCGen")+ fProcLabel[iproc] ,"Vz distribution of generated events vs rec multiplicity ");
+ fHistVzMCRec [iproc] = (TH2F*) BookVzHisto(TString("fHistVzMCRec")+ fProcLabel[iproc] ,"Vz distribution of reconstructed events vs rec multiplicity");
+ fHistVzMCTrg [iproc] = (TH2F*) BookVzHisto(TString("fHistVzMCTrg")+ fProcLabel[iproc] ,"Vz distribution of triggered events vs rec multiplicity ");
+ }
+ fHistVzData = (TH2F*) BookVzHisto("fHistVzData" ,"Vz distribution of triggered events vs rec multiplicity ");
+ fHistProcTypes = new TH1F ("fHistProcTypes", "Number of events in the different process classes", kNProcs, -0.5 , kNProcs-0.5);
+
+ fHistProcTypes->GetXaxis()->SetBinLabel(kProcSD+1,"SD");
+ fHistProcTypes->GetXaxis()->SetBinLabel(kProcND+1,"ND");
+ fHistProcTypes->GetXaxis()->SetBinLabel(kProcDD+1,"DD");
+ fHistProcTypes->GetXaxis()->SetBinLabel(kProcUnknown+1,"Unknown");
+
+ TH1::AddDirectory(oldStatus);
+
+}
+
+TH1 * AliCollisionNormalization::BookVzHisto(const char * name , const char * title, Bool_t vzOnly) {
+
+ // Book histos of vz distributions vs multiplicity
+ // if vzOnly == kTRUE, it returns a 1d histo with vz dist only
+
+
+ // Do not attach histos to the current directory
+ Bool_t oldStatus = TH1::AddDirectoryStatus();
+ TH1::AddDirectory(kFALSE);
+
+ TH1 * h;
+ Double_t binLimitsVtx[] = {-30,-25,-20,-15,-10,-7,-5.5,-4,-3,-2,-1,0,1,2,3,4,5.5,7,10,15,20,25,30};
+ if (vzOnly) {
+ h = new TH1F(name,title,22,binLimitsVtx);
+ } else {
+ h = new TH2F(name,title,22,binLimitsVtx,100,Double_t(-0.5),Double_t(99.5));
+ }
+
+ h->SetXTitle("V_{z} (cm)");
+ h->SetYTitle("n_{trk}");
+ h->Sumw2();
+
+ TH1::AddDirectory(oldStatus);
+
+ return h;
+
+}
+
+TH2F * AliCollisionNormalization::GetVzMCGen (Int_t procType) {
+
+ // returns MC gen histo. If proc type < 0 sums over all proc types, reweighting the XS
+
+ if(procType>=0) return fHistVzMCGen[procType] ;
+
+ TH2F * sum = (TH2F*) fHistVzMCGen[kProcSD]->Clone();
+ sum->Reset();
+
+ for(Int_t iproc = 0; iproc < kProcUnknown; iproc++){
+ sum->Add(fHistVzMCGen[iproc],GetProcessWeight(iproc));
+ }
+
+ return sum;
+}
+TH2F * AliCollisionNormalization::GetVzMCRec (Int_t procType) {
+
+ // returns MC rec histo. If proc type < 0 sums over all proc types, reweighting the XS
+
+ if(procType>=0) return fHistVzMCRec[procType] ;
+
+ TH2F * sum = (TH2F*) fHistVzMCRec[kProcSD]->Clone();
+ sum->Reset();
+
+ for(Int_t iproc = 0; iproc < kProcUnknown; iproc++){
+ sum->Add(fHistVzMCRec[iproc],GetProcessWeight(iproc));
+ }
+
+ return sum;
+
+}
+
+
+TH2F * AliCollisionNormalization::GetVzMCTrg (Int_t procType) {
+
+ // returns MC trg histo. If proc type < 0 sums over all proc types, reweighting the XS
+
+ if(procType>=0) return fHistVzMCTrg[procType] ;
+
+ TH2F * sum = (TH2F*) fHistVzMCTrg[kProcSD]->Clone();
+ sum->Reset();
+
+ for(Int_t iproc = 0; iproc < kProcUnknown; iproc++){
+ sum->Add(fHistVzMCTrg[iproc],GetProcessWeight(iproc));
+ }
+
+ return sum;
+
+}
+
+Double_t AliCollisionNormalization::ComputeNint() {
+
+ // Compute the number of collisions applying all corrections
+ // TODO: check error propagation
+
+ TH1 * hVzData = fHistVzData->ProjectionX("_px",2,-1,"E"); // Skip zero bin
+ Int_t allEventsWithVertex = (Int_t) fHistVzData->Integral(0, fHistVzData->GetNbinsX()+1,
+ 2, fHistVzData->GetNbinsY()+1); // include under/overflow!, skip 0 bin
+
+ // Assign histos reweighted with measured XS
+ TH2F * histVzMCGen = GetVzMCGen(-1);
+ TH2F * histVzMCTrg = GetVzMCTrg(-1);
+ TH2F * histVzMCRec = GetVzMCRec(-1);
+
+ // Get or compute BG. This assumes the CINT1B suite
+ Double_t triggeredEventsWith0MultWithBG = fHistVzData->Integral(0, fHistVzData->GetNbinsX()+1,1, 1);
+ Double_t bg = 0; // This will include beam gas + accidentals
+ if (fHistStatBin0->GetNbinsY() > 4) { // FIXME: we need a better criterion to decide...
+ AliInfo("Using BG computed by Physics Selection");
+ bg = fHistStatBin0->GetBinContent(fHistStatBin0->GetNbinsX(),AliPhysicsSelection::kStatRowBG);
+ bg += fHistStatBin0->GetBinContent(fHistStatBin0->GetNbinsX(),AliPhysicsSelection::kStatRowAcc);
+ } else {
+ AliInfo("Computing BG using CINT1A/B/C/E, ignoring intensities");
+ Int_t icol = fHistStatBin0->GetNbinsX();
+ Int_t cint1B = (Int_t) fHistStatBin0->GetBinContent(icol,1);
+ Int_t cint1A = (Int_t) fHistStatBin0->GetBinContent(icol,2);
+ Int_t cint1C = (Int_t) fHistStatBin0->GetBinContent(icol,3);
+ Int_t cint1E = (Int_t) fHistStatBin0->GetBinContent(icol,4);
+ bg = cint1A + cint1C-2*cint1E ;
+ if (cint1B != triggeredEventsWith0MultWithBG) {
+ AliWarning(Form("Events in bin0 from physics selection and local counter not consistent: %d - %d", cint1B, triggeredEventsWith0MultWithBG));
+ }
+ }
+
+ Double_t triggeredEventsWith0Mult = triggeredEventsWith0MultWithBG - bg;
+ if(fVerbose > 0) AliInfo(Form("Measured events with vertex: %d",allEventsWithVertex));
+ Double_t bin0 = fHistVzData->Integral(0, fHistVzData->GetNbinsX()+1,
+ 1, 1);
+ if(fVerbose > 0) AliInfo(Form("Zero Bin, Meas: %2.2f, BG: %2.2f, Meas - BG: %2.2f", bin0, bg, triggeredEventsWith0Mult));
+
+ // This pointers are here so that I use the same names used in Jan Fiete's code (allows for faster/easier comparison)
+
+ TH2* eTrig = histVzMCTrg;
+ TH2* eTrigVtx = histVzMCRec;
+ TH1* eTrigVtx_projx = eTrigVtx->ProjectionX("eTrigVtx_projx", 2, eTrigVtx->GetNbinsX()+1);
+
+ // compute trigger and vertex efficiency
+ TH2 * effVtxTrig = (TH2*) histVzMCRec->Clone("effVtxTrig");
+ effVtxTrig->Reset();
+ effVtxTrig->Divide(histVzMCRec,histVzMCGen,1,1,"B");
+ // Apply correction to data to get corrected events
+ TH2 * correctedEvents = (TH2*) fHistVzData->Clone("correctedEvents");
+ correctedEvents->Divide(effVtxTrig);
+
+ // TH1 * correctedEvents = fHistVzData->ProjectionX("eTrigVtx_projx", 2, eTrigVtx->GetNbinsX()+1);
+
+ // loop over vertex bins
+ for (Int_t i = 1; i <= fHistVzData->GetNbinsX(); i++)
+ {
+ Double_t alpha = (Double_t) hVzData->GetBinContent(i) / allEventsWithVertex;
+ Double_t events = alpha * triggeredEventsWith0Mult;
+
+ if (histVzMCRec->GetBinContent(i,1) == 0)
+ continue;
+
+ Double_t fZ = eTrigVtx_projx->Integral(0, eTrigVtx_projx->GetNbinsX()+1) / eTrigVtx_projx->GetBinContent(i) *
+ eTrig->GetBinContent(i, 1) / eTrig->Integral(0, eTrig->GetNbinsX()+1, 1, 1);
+
+ events *= fZ;
+
+ // multiply with trigger correction
+ Double_t correction = histVzMCGen->GetBinContent(i,1)/histVzMCTrg->GetBinContent(i,1);
+ events *= correction;
+
+ if (fVerbose > 1) Printf(" Bin %d, alpha is %.2f%%, fZ is %.3f, number of events with 0 mult.: %.2f (MC comparison: %.2f)", i, alpha * 100., fZ, events,
+ histVzMCRec->GetBinContent(i,1));
+ correctedEvents->SetBinContent(i, 1, events);
+ }
+
+
+
+ // Integrate correctedEvents over full z range
+ Double_t allEvents = correctedEvents->Integral(0, correctedEvents->GetNbinsX()+1,0, correctedEvents->GetNbinsY()+1);
+ // Integrate correctedEvents over needed z range
+ Double_t allEventsZrange = correctedEvents->Integral(correctedEvents->GetXaxis()->FindBin(-fZRange), correctedEvents->GetXaxis()->FindBin(fZRange),
+ 0, correctedEvents->GetNbinsY()+1);
+
+ if(fVerbose > 1) AliInfo(Form("Results in |Vz| < %3.3f",fZRange));
+ if(fVerbose > 1) AliInfo(Form(" Events in Bin0: %2.2f, With > 1 track: %2.2f, All corrected: %2.2f",
+ correctedEvents->Integral(correctedEvents->GetXaxis()->FindBin(-fZRange), correctedEvents->GetXaxis()->FindBin(fZRange),1,1),
+ correctedEvents->Integral(correctedEvents->GetXaxis()->FindBin(-fZRange), correctedEvents->GetXaxis()->FindBin(fZRange),
+ 2,correctedEvents->GetNbinsX()+1),
+ allEventsZrange));
+ if(fVerbose > 1) {
+ Int_t nbin = histVzMCRec->GetNbinsX();
+ AliInfo(Form("Efficiency in the zero bin: %3.3f", histVzMCRec->Integral(0,nbin+1,1,1)/histVzMCGen->Integral(0,nbin+1,1,1) ));
+ }
+
+ AliInfo(Form("Number of collisions in full phase space: %2.2f", allEvents));
+
+ return allEvents;
+}
+
+Long64_t AliCollisionNormalization::Merge(TCollection* list)
+{
+ // Merge a list of AliCollisionNormalization objects with this
+ // (needed for PROOF).
+ // Returns the number of merged objects (including this).
+
+ if (!list)
+ return 0;
+
+ if (list->IsEmpty())
+ return 1;
+
+ TIterator* iter = list->MakeIterator();
+ TObject* obj;
+
+ // collections of all histograms
+ const Int_t nHists = kNProcs*3+4;
+ TList collections[nHists];
+
+ Int_t count = 0;
+ while ((obj = iter->Next())) {
+
+ AliCollisionNormalization* entry = dynamic_cast<AliCollisionNormalization*> (obj);
+ if (entry == 0)
+ continue;
+ Int_t ihist = -1;
+
+ for(Int_t iproc = 0; iproc < kNProcs; iproc++){
+ if (entry->fHistVzMCGen[iproc] ) collections[++ihist].Add(entry->fHistVzMCGen[iproc] );
+ if (entry->fHistVzMCRec[iproc] ) collections[++ihist].Add(entry->fHistVzMCRec[iproc] );
+ if (entry->fHistVzMCTrg[iproc] ) collections[++ihist].Add(entry->fHistVzMCTrg[iproc] );
+ }
+ if (entry->fHistVzData ) collections[++ihist].Add(entry->fHistVzData );
+ if (entry->fHistProcTypes ) collections[++ihist].Add(entry->fHistProcTypes );
+ if (entry->fHistStatBin0 ) collections[++ihist].Add(entry->fHistStatBin0 );
+
+ count++;
+ }
+
+ Int_t ihist = -1;
+ for(Int_t iproc = 0; iproc < kNProcs; iproc++){
+ if (fHistVzMCGen[iproc] ) fHistVzMCGen[iproc] ->Merge(&collections[++ihist]);
+ if (fHistVzMCRec[iproc] ) fHistVzMCRec[iproc] ->Merge(&collections[++ihist]);
+ if (fHistVzMCTrg[iproc] ) fHistVzMCTrg[iproc] ->Merge(&collections[++ihist]);
+
+ }
+ if (fHistVzData ) fHistVzData ->Merge(&collections[++ihist]);
+ if (fHistProcTypes ) fHistProcTypes ->Merge(&collections[++ihist]);
+ if (fHistStatBin0 ) fHistStatBin0 ->Merge(&collections[++ihist]);
+
+
+ delete iter;
+
+ return count+1;
+}
+
+void AliCollisionNormalization::FillVzMCGen(Float_t vz, Int_t ntrk, AliMCEvent * mcEvt) {
+
+ // Fill MC gen histo and the process types statistics
+
+ Int_t evtype = GetProcessType(mcEvt);
+ // When I fill gen histos, I also fill statistics of process types (used to detemine ratios afterwards).
+ fHistProcTypes->Fill(evtype);
+ fHistVzMCGen[evtype]->Fill(vz,ntrk);
+}
+
+void AliCollisionNormalization::FillVzMCRec(Float_t vz, Int_t ntrk, AliMCEvent * mcEvt){
+
+ // Fill MC rec histo
+ Int_t evtype = GetProcessType(mcEvt);
+ fHistVzMCRec[evtype]->Fill(vz,ntrk);
+
+}
+void AliCollisionNormalization::FillVzMCTrg(Float_t vz, Int_t ntrk, AliMCEvent * mcEvt) {
+
+ // Fill MC trg histo
+ Int_t evtype = GetProcessType(mcEvt);
+ fHistVzMCTrg[evtype]->Fill(vz,ntrk);
+}
+
+
+Int_t AliCollisionNormalization::GetProcessType(AliMCEvent * mcEvt) {
+
+ // Determine if the event was generated with pythia or phojet and return the process type
+
+ // Check if mcEvt is fine
+ if (!mcEvt) {
+ AliFatal("NULL mc event");
+ }
+
+ // Determine if it was a pythia or phojet header, and return the correct process type
+ AliGenPythiaEventHeader * headPy = 0;
+ AliGenDPMjetEventHeader * headPho = 0;
+ AliGenEventHeader * htmp = mcEvt->GenEventHeader();
+ if(!htmp) {
+ AliFatal("Cannot Get MC Header!!");
+ }
+ if( TString(htmp->IsA()->GetName()) == "AliGenPythiaEventHeader") {
+ headPy = (AliGenPythiaEventHeader*) htmp;
+ } else if (TString(htmp->IsA()->GetName()) == "AliGenDPMjetEventHeader") {
+ headPho = (AliGenDPMjetEventHeader*) htmp;
+ } else {
+ AliError("Unknown header");
+ }
+
+ // Determine process type
+ if(headPy) {
+ if(headPy->ProcessType() == 92 || headPy->ProcessType() == 93) {
+ // single difractive
+ return kProcSD;
+ } else if (headPy->ProcessType() == 94) {
+ // double diffractive
+ return kProcDD;
+ }
+ else if(headPy->ProcessType() != 92 && headPy->ProcessType() != 93 && headPy->ProcessType() != 94) {
+ // non difractive
+ return kProcND;
+ }
+ } else if (headPho) {
+ if(headPho->ProcessType() == 5 || headPho->ProcessType() == 6 ) {
+ // single difractive
+ return kProcSD;
+ } else if (headPho->ProcessType() == 7) {
+ // double diffractive
+ return kProcDD;
+ } else if(headPho->ProcessType() != 5 && headPho->ProcessType() != 6 && headPho->ProcessType() != 7 ) {
+ // non difractive
+ return kProcND;
+ }
+ }
+
+
+ // no process type found?
+ AliError(Form("Unknown header: %s", htmp->IsA()->GetName()));
+ return kProcUnknown;
+}
+
+Double_t AliCollisionNormalization::GetProcessWeight(Int_t proc) {
+
+ // Return a weight to be used when combining the MC histos to
+ // compute efficiency, defined as the ratio XS measured / XS in
+ // generator
+
+ Float_t ref_SD, ref_DD, ref_ND, error_SD, error_DD, error_ND;
+ GetRelativeFractions(fReferenceXS,ref_SD, ref_DD, ref_ND, error_SD, error_DD, error_ND);
+
+ static Double_t total = fHistProcTypes->Integral();
+ if (fHistProcTypes->GetBinContent(fHistProcTypes->FindBin(kProcUnknown)) > 0) {
+ AliError("There were unknown processes!!!");
+ }
+ static Double_t SD = fHistProcTypes->GetBinContent(fHistProcTypes->FindBin(kProcSD))/total;
+ static Double_t DD = fHistProcTypes->GetBinContent(fHistProcTypes->FindBin(kProcDD))/total;
+ static Double_t ND = fHistProcTypes->GetBinContent(fHistProcTypes->FindBin(kProcND))/total;
+
+ if (fVerbose > 2) {
+ AliInfo(Form("Total MC evts: %f",total));
+ AliInfo(Form(" Frac SD %4.4f", SD));
+ AliInfo(Form(" Frac DD %4.4f", DD));
+ AliInfo(Form(" Frac ND %4.4f", ND));
+ }
+
+ switch(proc) {
+ case kProcSD:
+ return ref_SD/SD;
+ break;
+ case kProcDD:
+ return ref_DD/DD;
+ break;
+ case kProcND:
+ return ref_ND/ND;
+ break;
+ default:
+ AliError("Unknown process");
+ }
+
+ return 0;
+
+}
+
+void AliCollisionNormalization::GetRelativeFractions(Int_t origin, Float_t& ref_SD, Float_t& ref_DD, Float_t& ref_ND, Float_t& error_SD, Float_t& error_DD, Float_t& error_ND)
+{
+ // Returns fraction of XS (SD, ND, DD) and corresponding error
+ // Stolen from Jan Fiete's drawSystematics macro
+
+ // origin:
+ // -1 = Pythia (test)
+ // 0 = UA5
+ // 1 = Data 1.8 TeV
+ // 2 = Tel-Aviv
+ // 3 = Durham
+ //
+
+ switch (origin)
+ {
+ case -10: // Pythia default at 7 GeV, 50% error
+ AliInfo("PYTHIA x-sections");
+ ref_SD = 0.192637; error_SD = ref_SD * 0.5;
+ ref_DD = 0.129877; error_DD = ref_DD * 0.5;
+ ref_ND = 0.677486; error_ND = 0;
+ break;
+
+ case -1: // Pythia default at 900 GeV, as test
+ AliInfo("PYTHIA x-sections");
+ ref_SD = 0.223788;
+ ref_DD = 0.123315;
+ ref_ND = 0.652897;
+ break;
+
+ case 0: // UA5
+ AliInfo("UA5 x-sections a la first paper");
+ ref_SD = 0.153; error_SD = 0.05;
+ ref_DD = 0.080; error_DD = 0.05;
+ ref_ND = 0.767; error_ND = 0;
+ break;
+
+ case 10: // UA5
+ AliInfo("UA5 x-sections hadron level definition for Pythia");
+ // Fractions in Pythia with UA5 cuts selection for SD
+ // ND: 0.688662
+ // SD: 0.188588 --> this should be 15.3
+ // DD: 0.122750
+ ref_SD = 0.224 * 0.153 / 0.189; error_SD = 0.023 * 0.224 / 0.189;
+ ref_DD = 0.095; error_DD = 0.06;
+ ref_ND = 1.0 - ref_SD - ref_DD; error_ND = 0;
+ break;
+
+ case 11: // UA5
+ AliInfo("UA5 x-sections hadron level definition for Phojet");
+ // Fractions in Phojet with UA5 cuts selection for SD
+ // ND: 0.783573
+ // SD: 0.151601 --> this should be 15.3
+ // DD: 0.064827
+ ref_SD = 0.191 * 0.153 / 0.152; error_SD = 0.023 * 0.191 / 0.152;
+ ref_DD = 0.095; error_DD = 0.06;
+ ref_ND = 1.0 - ref_SD - ref_DD; error_ND = 0;
+ break;
+
+ case 20: // E710, 1.8 TeV
+ AliInfo("E710 x-sections hadron level definition for Pythia");
+ // ND: 0.705709
+ // SD: 0.166590 --> this should be 15.9
+ // DD: 0.127701
+ ref_SD = 0.217 * 0.159 / 0.167; error_SD = 0.024 * 0.217 / 0.167;
+ ref_DD = 0.075 * 1.43; error_DD = 0.02 * 1.43;
+ ref_ND = 1.0 - ref_SD - ref_DD; error_ND = 0;
+ break;
+
+ case 21: // E710, 1.8 TeV
+ AliInfo("E710 x-sections hadron level definition for Phojet");
+ // ND: 0.817462
+ // SD: 0.125506 --> this should be 15.9
+ // DD: 0.057032
+ ref_SD = 0.161 * 0.159 / 0.126; error_SD = 0.024 * 0.161 / 0.126;
+ ref_DD = 0.075 * 1.43; error_DD = 0.02 * 1.43;
+ ref_ND = 1.0 - ref_SD - ref_DD; error_ND = 0;
+ break;
+
+ case 1: // data 1.8 TeV
+ AliInfo("??? x-sections");
+ ref_SD = 0.152;
+ ref_DD = 0.092;
+ ref_ND = 1 - ref_SD - ref_DD;
+ break;
+
+ case 2: // tel-aviv model
+ AliInfo("Tel-aviv model x-sections");
+ ref_SD = 0.171;
+ ref_DD = 0.094;
+ ref_ND = 1 - ref_SD - ref_DD;
+ break;
+
+ case 3: // durham model
+ AliInfo("Durham model x-sections");
+ ref_SD = 0.190;
+ ref_DD = 0.125;
+ ref_ND = 1 - ref_SD - ref_DD;
+ break;
+
+ default:
+ AliFatal(Form("Unknown origin %d", origin));
+ }
+}
+
git://git.uio.no / u / mrichter / AliRoot.git / blobdiff