From: kleinb Date: Thu, 23 Apr 2009 11:29:02 +0000 (+0000) Subject: Jet Jet Production 3 p_T hard bins X-Git-Url: http://git.uio.no/git/?p=u%2Fmrichter%2FAliRoot.git;a=commitdiff_plain;h=afea418a761beebfc7e8bd391cad69c9c7800aa3 Jet Jet Production 3 p_T hard bins --- diff --git a/prod/LHC09a1/CheckESD.C b/prod/LHC09a1/CheckESD.C new file mode 100644 index 00000000000..badd79f8c66 --- /dev/null +++ b/prod/LHC09a1/CheckESD.C @@ -0,0 +1,699 @@ +#if !defined( __CINT__) || defined(__MAKECINT__) +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#include "AliRunLoader.h" +#include "AliLoader.h" +#include "AliESDEvent.h" +#include "AliRun.h" +#include "AliStack.h" +#include "AliHeader.h" +#include "AliGenEventHeader.h" +#include "AliPID.h" +#else +// const Int_t kXiMinus = 3312; +// const Int_t kOmegaMinus = 3334; +#endif + + +TH1F* CreateHisto(const char* name, const char* title, + Int_t nBins, Double_t xMin, Double_t xMax, + const char* xLabel = NULL, const char* yLabel = NULL) +{ +// create a histogram + + TH1F* result = new TH1F(name, title, nBins, xMin, xMax); + result->SetOption("E"); + if (xLabel) result->GetXaxis()->SetTitle(xLabel); + if (yLabel) result->GetYaxis()->SetTitle(yLabel); + result->SetMarkerStyle(kFullCircle); + return result; +} + +TH1F* CreateEffHisto(TH1F* hGen, TH1F* hRec) +{ +// create an efficiency histogram + + Int_t nBins = hGen->GetNbinsX(); + TH1F* hEff = (TH1F*) hGen->Clone("hEff"); + hEff->SetTitle(""); + hEff->SetStats(kFALSE); + hEff->SetMinimum(0.); + hEff->SetMaximum(110.); + hEff->GetYaxis()->SetTitle("#epsilon [%]"); + + for (Int_t iBin = 0; iBin <= nBins; iBin++) { + Double_t nGen = hGen->GetBinContent(iBin); + Double_t nRec = hRec->GetBinContent(iBin); + if (nGen > 0) { + Double_t eff = nRec/nGen; + hEff->SetBinContent(iBin, 100. * eff); + Double_t error = sqrt(eff*(1.-eff) / nGen); + if (error == 0) error = 0.0001; + hEff->SetBinError(iBin, 100. * error); + } else { + hEff->SetBinContent(iBin, -100.); + hEff->SetBinError(iBin, 0); + } + } + + return hEff; +} + +Bool_t FitHisto(TH1* histo, Double_t& res, Double_t& resError) +{ +// fit a gaussian to a histogram + + static TF1* fitFunc = new TF1("fitFunc", "gaus"); + fitFunc->SetLineWidth(2); + fitFunc->SetFillStyle(0); + Double_t maxFitRange = 2; + + if (histo->Integral() > 50) { + Float_t mean = histo->GetMean(); + Float_t rms = histo->GetRMS(); + fitFunc->SetRange(mean - maxFitRange*rms, mean + maxFitRange*rms); + fitFunc->SetParameters(mean, rms); + histo->Fit(fitFunc, "QRI0"); + histo->GetFunction("fitFunc")->ResetBit(1<<9); + res = TMath::Abs(fitFunc->GetParameter(2)); + resError = TMath::Abs(fitFunc->GetParError(2)); + return kTRUE; + } + + return kFALSE; +} + + +Bool_t CheckESD(const char* gAliceFileName = "galice.root", + const char* esdFileName = "AliESDs.root") +{ +// check the content of the ESD + + // check values + Int_t checkNGenLow = 1; + + Double_t checkEffLow = 0.5; + Double_t checkEffSigma = 3; + Double_t checkFakeHigh = 0.5; + Double_t checkFakeSigma = 3; + + Double_t checkResPtInvHigh = 5; + Double_t checkResPtInvSigma = 3; + Double_t checkResPhiHigh = 10; + Double_t checkResPhiSigma = 3; + Double_t checkResThetaHigh = 10; + Double_t checkResThetaSigma = 3; + + Double_t checkPIDEffLow = 0.5; + Double_t checkPIDEffSigma = 3; + Double_t checkResTOFHigh = 500; + Double_t checkResTOFSigma = 3; + + Double_t checkPHOSNLow = 5; + Double_t checkPHOSEnergyLow = 0.3; + Double_t checkPHOSEnergyHigh = 1.0; + Double_t checkEMCALNLow = 50; + Double_t checkEMCALEnergyLow = 0.05; + Double_t checkEMCALEnergyHigh = 1.0; + + Double_t checkMUONNLow = 1; + Double_t checkMUONPtLow = 0.5; + Double_t checkMUONPtHigh = 10.; + + Double_t cutPtV0 = 0.3; + Double_t checkV0EffLow = 0.02; + Double_t checkV0EffSigma = 3; + Double_t cutPtCascade = 0.5; + Double_t checkCascadeEffLow = 0.01; + Double_t checkCascadeEffSigma = 3; + + // open run loader and load gAlice, kinematics and header + AliRunLoader* runLoader = AliRunLoader::Open(gAliceFileName); + if (!runLoader) { + Error("CheckESD", "getting run loader from file %s failed", + gAliceFileName); + return kFALSE; + } + runLoader->LoadgAlice(); + gAlice = runLoader->GetAliRun(); + if (!gAlice) { + Error("CheckESD", "no galice object found"); + return kFALSE; + } + runLoader->LoadKinematics(); + runLoader->LoadHeader(); + + // open the ESD file + TFile* esdFile = TFile::Open(esdFileName); + if (!esdFile || !esdFile->IsOpen()) { + Error("CheckESD", "opening ESD file %s failed", esdFileName); + return kFALSE; + } + AliESDEvent* esd = new AliESDEvent; + TTree* tree = (TTree*) esdFile->Get("esdTree"); + if (!tree) { + Error("CheckESD", "no ESD tree found"); + return kFALSE; + } + esd->ReadFromTree(tree); + + // efficiency and resolution histograms + Int_t nBinsPt = 15; + Float_t minPt = 0.1; + Float_t maxPt = 3.1; + TH1F* hGen = CreateHisto("hGen", "generated tracks", + nBinsPt, minPt, maxPt, "p_{t} [GeV/c]", "N"); + TH1F* hRec = CreateHisto("hRec", "reconstructed tracks", + nBinsPt, minPt, maxPt, "p_{t} [GeV/c]", "N"); + Int_t nGen = 0; + Int_t nRec = 0; + Int_t nFake = 0; + + TH1F* hResPtInv = CreateHisto("hResPtInv", "", 100, -10, 10, + "(p_{t,rec}^{-1}-p_{t,sim}^{-1}) / p_{t,sim}^{-1} [%]", "N"); + TH1F* hResPhi = CreateHisto("hResPhi", "", 100, -20, 20, + "#phi_{rec}-#phi_{sim} [mrad]", "N"); + TH1F* hResTheta = CreateHisto("hResTheta", "", 100, -20, 20, + "#theta_{rec}-#theta_{sim} [mrad]", "N"); + + // PID + Int_t partCode[] = + {kElectron, kMuonMinus, kPiPlus, kKPlus, kProton}; + const char* partName[] = + {"electron", "muon", "pion", "kaon", "proton", "other"}; + Double_t partFrac[] = + {0.01, 0.01, 0.85, 0.10, 0.05}; + Int_t identified[6][5]; + for (Int_t iGen = 0; iGen < 6; iGen++) { + for (Int_t iRec = 0; iRec < 5; iRec++) { + identified[iGen][iRec] = 0; + } + } + Int_t nIdentified = 0; + + // dE/dx and TOF + TH2F* hDEdxRight = new TH2F("hDEdxRight", "", 300, 0, 3, 100, 0, 400); + hDEdxRight->SetStats(kFALSE); + hDEdxRight->GetXaxis()->SetTitle("p [GeV/c]"); + hDEdxRight->GetYaxis()->SetTitle("dE/dx_{TPC}"); + hDEdxRight->SetMarkerStyle(kFullCircle); + hDEdxRight->SetMarkerSize(0.4); + TH2F* hDEdxWrong = new TH2F("hDEdxWrong", "", 300, 0, 3, 100, 0, 400); + hDEdxWrong->SetStats(kFALSE); + hDEdxWrong->GetXaxis()->SetTitle("p [GeV/c]"); + hDEdxWrong->GetYaxis()->SetTitle("dE/dx_{TPC}"); + hDEdxWrong->SetMarkerStyle(kFullCircle); + hDEdxWrong->SetMarkerSize(0.4); + hDEdxWrong->SetMarkerColor(kRed); + TH1F* hResTOFRight = CreateHisto("hResTOFRight", "", 100, -1000, 1000, + "t_{TOF}-t_{track} [ps]", "N"); + TH1F* hResTOFWrong = CreateHisto("hResTOFWrong", "", 100, -1000, 1000, + "t_{TOF}-t_{track} [ps]", "N"); + hResTOFWrong->SetLineColor(kRed); + + // calorimeters + TH1F* hEPHOS = CreateHisto("hEPHOS", "PHOS", 100, 0, 5, "E [GeV]", "N"); + TH1F* hEEMCAL = CreateHisto("hEEMCAL", "EMCAL", 100, 0, 50, "E [GeV]", "N"); + + // muons + TH1F* hPtMUON = CreateHisto("hPtMUON", "MUON", 100, 0, 20, + "p_{t} [GeV/c]", "N"); + + // V0s and cascades + TH1F* hMassK0 = CreateHisto("hMassK0", "K^{0}", 100, 0.4, 0.6, + "M(#pi^{+}#pi^{-}) [GeV/c^{2}]", "N"); + TH1F* hMassLambda = CreateHisto("hMassLambda", "#Lambda", 100, 1.0, 1.2, + "M(p#pi^{-}) [GeV/c^{2}]", "N"); + TH1F* hMassLambdaBar = CreateHisto("hMassLambdaBar", "#bar{#Lambda}", + 100, 1.0, 1.2, + "M(#bar{p}#pi^{+}) [GeV/c^{2}]", "N"); + Int_t nGenV0s = 0; + Int_t nRecV0s = 0; + TH1F* hMassXi = CreateHisto("hMassXi", "#Xi", 100, 1.2, 1.5, + "M(#Lambda#pi) [GeV/c^{2}]", "N"); + TH1F* hMassOmega = CreateHisto("hMassOmega", "#Omega", 100, 1.5, 1.8, + "M(#LambdaK) [GeV/c^{2}]", "N"); + Int_t nGenCascades = 0; + Int_t nRecCascades = 0; + + // loop over events + for (Int_t iEvent = 0; iEvent < runLoader->GetNumberOfEvents(); iEvent++) { + runLoader->GetEvent(iEvent); + + // select simulated primary particles, V0s and cascades + AliStack* stack = gAlice->Stack(); + Int_t nParticles = stack->GetNtrack(); + TArrayF vertex(3); + runLoader->GetHeader()->GenEventHeader()->PrimaryVertex(vertex); + TObjArray selParticles; + TObjArray selV0s; + TObjArray selCascades; + for (Int_t iParticle = 0; iParticle < nParticles; iParticle++) { + TParticle* particle = stack->Particle(iParticle); + if (!particle) continue; + if (particle->Pt() < 0.001) continue; + if (TMath::Abs(particle->Eta()) > 0.9) continue; + TVector3 dVertex(particle->Vx() - vertex[0], + particle->Vy() - vertex[1], + particle->Vz() - vertex[2]); + if (dVertex.Mag() > 0.0001) continue; + + switch (TMath::Abs(particle->GetPdgCode())) { + case kElectron: + case kMuonMinus: + case kPiPlus: + case kKPlus: + case kProton: { + if (particle->Pt() > minPt) { + selParticles.Add(particle); + nGen++; + hGen->Fill(particle->Pt()); + } + break; + } + case kK0Short: + case kLambda0: { + if (particle->Pt() > cutPtV0) { + nGenV0s++; + selV0s.Add(particle); + } + break; + } + case kXiMinus: + case kOmegaMinus: { + if (particle->Pt() > cutPtCascade) { + nGenCascades++; + selCascades.Add(particle); + } + break; + } + default: break; + } + } + + // get the event summary data + tree->GetEvent(iEvent); + if (!esd) { + Error("CheckESD", "no ESD object found for event %d", iEvent); + return kFALSE; + } + + // loop over tracks + for (Int_t iTrack = 0; iTrack < esd->GetNumberOfTracks(); iTrack++) { + AliESDtrack* track = esd->GetTrack(iTrack); + + // select tracks of selected particles + Int_t label = TMath::Abs(track->GetLabel()); + if (label > stack->GetNtrack()) continue; // background + TParticle* particle = stack->Particle(label); + if (!selParticles.Contains(particle)) continue; + if ((track->GetStatus() & AliESDtrack::kITSrefit) == 0) continue; + if (track->GetConstrainedChi2() > 1e9) continue; + selParticles.Remove(particle); // don't count multiple tracks + + nRec++; + hRec->Fill(particle->Pt()); + if (track->GetLabel() < 0) nFake++; + + // resolutions + Double_t p[3]; + track->GetConstrainedPxPyPz(p); + TVector3 pTrack(p); + hResPtInv->Fill(100. * (1./pTrack.Pt() - 1./particle->Pt()) * + particle->Pt()); + hResPhi->Fill(1000. * (pTrack.Phi() - particle->Phi())); + hResTheta->Fill(1000. * (pTrack.Theta() - particle->Theta())); + + // PID + if ((track->GetStatus() & AliESDtrack::kESDpid) == 0) continue; + Int_t iGen = 5; + for (Int_t i = 0; i < 5; i++) { + if (TMath::Abs(particle->GetPdgCode()) == partCode[i]) iGen = i; + } + Double_t probability[5]; + track->GetESDpid(probability); + Double_t pMax = 0; + Int_t iRec = 0; + for (Int_t i = 0; i < 5; i++) { + probability[i] *= partFrac[i]; + if (probability[i] > pMax) { + pMax = probability[i]; + iRec = i; + } + } + identified[iGen][iRec]++; + if (iGen == iRec) nIdentified++; + + // dE/dx and TOF + Double_t time[5]; + track->GetIntegratedTimes(time); + if (iGen == iRec) { + hDEdxRight->Fill(pTrack.Mag(), track->GetTPCsignal()); + if ((track->GetStatus() & AliESDtrack::kTOFpid) != 0) { + hResTOFRight->Fill(track->GetTOFsignal() - time[iRec]); + } + } else { + hDEdxWrong->Fill(pTrack.Mag(), track->GetTPCsignal()); + if ((track->GetStatus() & AliESDtrack::kTOFpid) != 0) { + hResTOFWrong->Fill(track->GetTOFsignal() - time[iRec]); + } + } + } + + // loop over muon tracks + { + for (Int_t iTrack = 0; iTrack < esd->GetNumberOfMuonTracks(); iTrack++) { + AliESDMuonTrack* muonTrack = esd->GetMuonTrack(iTrack); + Double_t ptInv = TMath::Abs(muonTrack->GetInverseBendingMomentum()); + if (ptInv > 0.001) { + hPtMUON->Fill(1./ptInv); + } + } + } + + // loop over V0s + for (Int_t iV0 = 0; iV0 < esd->GetNumberOfV0s(); iV0++) { + AliESDv0* v0 = esd->GetV0(iV0); + if (v0->GetOnFlyStatus()) continue; + v0->ChangeMassHypothesis(kK0Short); + hMassK0->Fill(v0->GetEffMass()); + v0->ChangeMassHypothesis(kLambda0); + hMassLambda->Fill(v0->GetEffMass()); + v0->ChangeMassHypothesis(kLambda0Bar); + hMassLambdaBar->Fill(v0->GetEffMass()); + + Int_t negLabel = TMath::Abs(esd->GetTrack(v0->GetNindex())->GetLabel()); + if (negLabel > stack->GetNtrack()) continue; // background + Int_t negMother = stack->Particle(negLabel)->GetMother(0); + if (negMother < 0) continue; + Int_t posLabel = TMath::Abs(esd->GetTrack(v0->GetPindex())->GetLabel()); + if (posLabel > stack->GetNtrack()) continue; // background + Int_t posMother = stack->Particle(posLabel)->GetMother(0); + if (negMother != posMother) continue; + TParticle* particle = stack->Particle(negMother); + if (!selV0s.Contains(particle)) continue; + selV0s.Remove(particle); + nRecV0s++; + } + + // loop over Cascades + for (Int_t iCascade = 0; iCascade < esd->GetNumberOfCascades(); + iCascade++) { + AliESDcascade* cascade = esd->GetCascade(iCascade); + Double_t v0q; + cascade->ChangeMassHypothesis(v0q,kXiMinus); + hMassXi->Fill(cascade->GetEffMass()); + cascade->ChangeMassHypothesis(v0q,kOmegaMinus); + hMassOmega->Fill(cascade->GetEffMass()); + + Int_t negLabel = TMath::Abs(esd->GetTrack(cascade->GetNindex()) + ->GetLabel()); + if (negLabel > stack->GetNtrack()) continue; // background + Int_t negMother = stack->Particle(negLabel)->GetMother(0); + if (negMother < 0) continue; + Int_t posLabel = TMath::Abs(esd->GetTrack(cascade->GetPindex()) + ->GetLabel()); + if (posLabel > stack->GetNtrack()) continue; // background + Int_t posMother = stack->Particle(posLabel)->GetMother(0); + if (negMother != posMother) continue; + Int_t v0Mother = stack->Particle(negMother)->GetMother(0); + if (v0Mother < 0) continue; + Int_t bacLabel = TMath::Abs(esd->GetTrack(cascade->GetBindex()) + ->GetLabel()); + if (bacLabel > stack->GetNtrack()) continue; // background + Int_t bacMother = stack->Particle(bacLabel)->GetMother(0); + if (v0Mother != bacMother) continue; + TParticle* particle = stack->Particle(v0Mother); + if (!selCascades.Contains(particle)) continue; + selCascades.Remove(particle); + nRecCascades++; + } + + // loop over the calorimeters clusters + { + Int_t nCaloCluster = esd->GetNumberOfCaloClusters(); + cout<<"CaloClusters "<GetCaloCluster(iCluster)->IsPHOS()) + hEPHOS->Fill(esd->GetCaloCluster(iCluster)->E()); + if(esd->GetCaloCluster(iCluster)->GetClusterType()==AliESDCaloCluster::kEMCALClusterv1) { + hEEMCAL->Fill(esd->GetCaloCluster(iCluster)->E()); + //cout<GetCaloCluster(iCluster)->E()< 0) { + // efficiency + Double_t eff = nRec*1./nGen; + Double_t effError = TMath::Sqrt(eff*(1.-eff) / nGen); + Double_t fake = nFake*1./nGen; + Double_t fakeError = TMath::Sqrt(fake*(1.-fake) / nGen); + Info("CheckESD", "eff = (%.1f +- %.1f) %% fake = (%.1f +- %.1f) %%", + 100.*eff, 100.*effError, 100.*fake, 100.*fakeError); + + if (eff < checkEffLow - checkEffSigma*effError) { + Warning("CheckESD", "low efficiency: (%.1f +- %.1f) %%", + 100.*eff, 100.*effError); + } + if (fake > checkFakeHigh + checkFakeSigma*fakeError) { + Warning("CheckESD", "high fake: (%.1f +- %.1f) %%", + 100.*fake, 100.*fakeError); + } + + // resolutions + Double_t res, resError; + if (FitHisto(hResPtInv, res, resError)) { + Info("CheckESD", "relative inverse pt resolution = (%.1f +- %.1f) %%", + res, resError); + if (res > checkResPtInvHigh + checkResPtInvSigma*resError) { + Warning("CheckESD", "bad pt resolution: (%.1f +- %.1f) %%", + res, resError); + } + } + + if (FitHisto(hResPhi, res, resError)) { + Info("CheckESD", "phi resolution = (%.1f +- %.1f) mrad", res, resError); + if (res > checkResPhiHigh + checkResPhiSigma*resError) { + Warning("CheckESD", "bad phi resolution: (%.1f +- %.1f) mrad", + res, resError); + } + } + + if (FitHisto(hResTheta, res, resError)) { + Info("CheckESD", "theta resolution = (%.1f +- %.1f) mrad", + res, resError); + if (res > checkResThetaHigh + checkResThetaSigma*resError) { + Warning("CheckESD", "bad theta resolution: (%.1f +- %.1f) mrad", + res, resError); + } + } + + // PID + if (nRec > 0) { + Double_t eff = nIdentified*1./nRec; + Double_t effError = TMath::Sqrt(eff*(1.-eff) / nRec); + Info("CheckESD", "PID eff = (%.1f +- %.1f) %%", + 100.*eff, 100.*effError); + if (eff < checkPIDEffLow - checkPIDEffSigma*effError) { + Warning("CheckESD", "low PID efficiency: (%.1f +- %.1f) %%", + 100.*eff, 100.*effError); + } + } + + printf("%9s:", "gen\\rec"); + for (Int_t iRec = 0; iRec < AliPID::kSPECIES; iRec++) { + printf("%9s", partName[iRec]); + } + printf("\n"); + for (Int_t iGen = 0; iGen < AliPID::kSPECIES+1; iGen++) { + printf("%9s:", partName[iGen]); + for (Int_t iRec = 0; iRec < AliPID::kSPECIES; iRec++) { + printf("%9d", identified[iGen][iRec]); + } + printf("\n"); + } + + if (FitHisto(hResTOFRight, res, resError)) { + Info("CheckESD", "TOF resolution = (%.1f +- %.1f) ps", res, resError); + if (res > checkResTOFHigh + checkResTOFSigma*resError) { + Warning("CheckESD", "bad TOF resolution: (%.1f +- %.1f) ps", + res, resError); + } + } + + // calorimeters + if (hEPHOS->Integral() < checkPHOSNLow) { + Warning("CheckESD", "low number of PHOS particles: %d", + Int_t(hEPHOS->Integral())); + } else { + Double_t mean = hEPHOS->GetMean(); + if (mean < checkPHOSEnergyLow) { + Warning("CheckESD", "low mean PHOS energy: %.1f GeV", mean); + } else if (mean > checkPHOSEnergyHigh) { + Warning("CheckESD", "high mean PHOS energy: %.1f GeV", mean); + } + } + + if (hEEMCAL->Integral() < checkEMCALNLow) { + Warning("CheckESD", "low number of EMCAL particles: %d", + Int_t(hEEMCAL->Integral())); + } else { + Double_t mean = hEEMCAL->GetMean(); + if (mean < checkEMCALEnergyLow) { + Warning("CheckESD", "low mean EMCAL energy: %.1f GeV", mean); + } else if (mean > checkEMCALEnergyHigh) { + Warning("CheckESD", "high mean EMCAL energy: %.1f GeV", mean); + } + } + + // muons + if (hPtMUON->Integral() < checkMUONNLow) { + Warning("CheckESD", "low number of MUON particles: %d", + Int_t(hPtMUON->Integral())); + } else { + Double_t mean = hPtMUON->GetMean(); + if (mean < checkMUONPtLow) { + Warning("CheckESD", "low mean MUON pt: %.1f GeV/c", mean); + } else if (mean > checkMUONPtHigh) { + Warning("CheckESD", "high mean MUON pt: %.1f GeV/c", mean); + } + } + + // V0s + if (nGenV0s > 0) { + Double_t eff = nRecV0s*1./nGenV0s; + Double_t effError = TMath::Sqrt(eff*(1.-eff) / nGenV0s); + if (effError == 0) effError = checkV0EffLow / TMath::Sqrt(1.*nGenV0s); + Info("CheckESD", "V0 eff = (%.1f +- %.1f) %%", + 100.*eff, 100.*effError); + if (eff < checkV0EffLow - checkV0EffSigma*effError) { + Warning("CheckESD", "low V0 efficiency: (%.1f +- %.1f) %%", + 100.*eff, 100.*effError); + } + } + + // Cascades + if (nGenCascades > 0) { + Double_t eff = nRecCascades*1./nGenCascades; + Double_t effError = TMath::Sqrt(eff*(1.-eff) / nGenCascades); + if (effError == 0) effError = checkV0EffLow / + TMath::Sqrt(1.*nGenCascades); + Info("CheckESD", "Cascade eff = (%.1f +- %.1f) %%", + 100.*eff, 100.*effError); + if (eff < checkCascadeEffLow - checkCascadeEffSigma*effError) { + Warning("CheckESD", "low Cascade efficiency: (%.1f +- %.1f) %%", + 100.*eff, 100.*effError); + } + } + } + + // draw the histograms if not in batch mode + if (!gROOT->IsBatch()) { + new TCanvas; + hEff->DrawCopy(); + new TCanvas; + hResPtInv->DrawCopy("E"); + new TCanvas; + hResPhi->DrawCopy("E"); + new TCanvas; + hResTheta->DrawCopy("E"); + new TCanvas; + hDEdxRight->DrawCopy(); + hDEdxWrong->DrawCopy("SAME"); + new TCanvas; + hResTOFRight->DrawCopy("E"); + hResTOFWrong->DrawCopy("SAME"); + new TCanvas; + hEPHOS->DrawCopy("E"); + new TCanvas; + hEEMCAL->DrawCopy("E"); + new TCanvas; + hPtMUON->DrawCopy("E"); + new TCanvas; + hMassK0->DrawCopy("E"); + new TCanvas; + hMassLambda->DrawCopy("E"); + new TCanvas; + hMassLambdaBar->DrawCopy("E"); + new TCanvas; + hMassXi->DrawCopy("E"); + new TCanvas; + hMassOmega->DrawCopy("E"); + } + + // write the output histograms to a file + TFile* outputFile = TFile::Open("check.root", "recreate"); + if (!outputFile || !outputFile->IsOpen()) { + Error("CheckESD", "opening output file check.root failed"); + return kFALSE; + } + hEff->Write(); + hResPtInv->Write(); + hResPhi->Write(); + hResTheta->Write(); + hDEdxRight->Write(); + hDEdxWrong->Write(); + hResTOFRight->Write(); + hResTOFWrong->Write(); + hEPHOS->Write(); + hEEMCAL->Write(); + hPtMUON->Write(); + hMassK0->Write(); + hMassLambda->Write(); + hMassLambdaBar->Write(); + hMassXi->Write(); + hMassOmega->Write(); + outputFile->Close(); + delete outputFile; + + // clean up + delete hGen; + delete hRec; + delete hEff; + delete hResPtInv; + delete hResPhi; + delete hResTheta; + delete hDEdxRight; + delete hDEdxWrong; + delete hResTOFRight; + delete hResTOFWrong; + delete hEPHOS; + delete hEEMCAL; + delete hPtMUON; + delete hMassK0; + delete hMassLambda; + delete hMassLambdaBar; + delete hMassXi; + delete hMassOmega; + + delete esd; + esdFile->Close(); + delete esdFile; + + runLoader->UnloadHeader(); + runLoader->UnloadKinematics(); + delete runLoader; + + // result of check + Info("CheckESD", "check of ESD was successfull"); + return kTRUE; +} diff --git a/prod/LHC09a1/Config.C b/prod/LHC09a1/Config.C new file mode 100644 index 00000000000..5e108b7503f --- /dev/null +++ b/prod/LHC09a1/Config.C @@ -0,0 +1,1025 @@ +// +// Configuration for the Physics Data Challenge 2007 +// +// +// Assuming inel = 70 mb (PPR v.1, p.64) +// +// 84.98% of MSEL=0 events (including diffractive) with +// QQbar switched off (these events will include injected J/psi) => kPyMbNoHvq +// +// 14.14% of MSEL=1 events with ccbar (in 4 subsamples) => kCharmpp14000wmi +// bin 1 25% (3.535%): 2.76 < pthard < 3 GeV/c +// bin 2 40% (5.656%): 3 < pthard < 4 GeV/c +// bin 3 29% (4.101%): 4 < pthard < 8 GeV/c +// bin 4 6% (0.848%): pthard > 8 GeV/c +// +// 0.73% of MSEL=1 events with bbbar (in 4 subsamples) => kBeautypp14000wmi +// bin 1 5% (0.037%): 2.76 < pthard < 4 GeV/c +// bin 2 31% (0.226%): 4 < pthard < 6 GeV/c +// bin 3 28% (0.204%): 6 < pthard < 8 GeV/c +// bin 4 36% (0.263%): pthard >8 GeV/c +// +// 0.075% of MSEL=0 events with QQbar switched off and 1 Omega- => kPyOmegaMinus +// 0.075% of MSEL=0 events with QQbar switched off and 1 OmegaBar+ => kPyOmegaPlus + +// One can use the configuration macro in compiled mode by +// root [0] gSystem->Load("libgeant321"); +// root [0] gSystem->SetIncludePath("-I$ROOTSYS/include -I$ALICE_ROOT/include\ +// -I$ALICE_ROOT -I$ALICE/geant3/TGeant3"); +// root [0] .x grun.C(1,"Config_PDC07.C++") + +class AliGenPythia; + +#if !defined(__CINT__) || defined(__MAKECINT__) +#include +#include +#include +#include +#include +#include +#include "EVGEN/AliGenCocktail.h" +#include "EVGEN/AliGenParam.h" +#include "EVGEN/AliGenMUONlib.h" +#include "STEER/AliRunLoader.h" +#include "STEER/AliRun.h" +#include "STEER/AliConfig.h" +#include "PYTHIA6/AliDecayerPythia.h" +#include "PYTHIA6/AliGenPythia.h" +#include "STEER/AliMagFMaps.h" +#include "STRUCT/AliBODY.h" +#include "STRUCT/AliMAG.h" +#include "STRUCT/AliABSOv3.h" +#include "STRUCT/AliDIPOv3.h" +#include "STRUCT/AliHALLv3.h" +#include "STRUCT/AliFRAMEv2.h" +#include "STRUCT/AliSHILv3.h" +#include "STRUCT/AliPIPEv3.h" +#include "ITS/AliITSgeom.h" +#include "ITS/AliITSv11Hybrid.h" +#include "TPC/AliTPCv2.h" +#include "TOF/AliTOFv6T0.h" +#include "HMPID/AliHMPIDv3.h" +#include "HMPID/AliHMPIDv2.h" +#include "ZDC/AliZDCv3.h" +#include "TRD/AliTRDv1.h" +#include "FMD/AliFMDv1.h" +#include "MUON/AliMUONv1.h" +#include "PHOS/AliPHOSv1.h" +#include "PMD/AliPMDv1.h" +#include "T0/AliT0v1.h" +#include "EMCAL/AliEMCALv2.h" +#include "ACORDE/AliACORDEv1.h" +#include "VZERO/AliVZEROv7.h" +#endif + + + // libraries required by geant321 + + + +enum PDC07Proc_t +{ +//--- Heavy Flavour Production --- + kCharmPbPb5500, kCharmpPb8800, kCharmpp14000, kCharmpp14000wmi, + kD0PbPb5500, kD0pPb8800, kD0pp14000, + kDPlusPbPb5500, kDPluspPb8800, kDPluspp14000, + kBeautyPbPb5500, kBeautypPb8800, kBeautypp14000, kBeautypp14000wmi, +// -- Pythia Mb + kPyMbNoHvq, kPyOmegaPlus, kPyOmegaMinus, kPyJetJet, kPyJetJetFixed, + kPyGammaJetPHOS, kPyJetJetPHOS, kPyJetJetPHOSv2, kPyGammaBremsPHOS, + kPyGammaJetEMCAL, kPyJetJetEMCAL, kPyGammaBremsEMCAL, kRunMax +}; + +const char * pprRunName[] = { + "kCharmPbPb5500", "kCharmpPb8800", "kCharmpp14000", "kCharmpp14000wmi", + "kD0PbPb5500", "kD0pPb8800", "kD0pp14000", + "kDPlusPbPb5500", "kDPluspPb8800", "kDPluspp14000", + "kBeautyPbPb5500", "kBeautypPb8800", "kBeautypp14000", "kBeautypp14000wmi", + "kPyMbNoHvq", "kPyOmegaPlus", "kPyOmegaMinus", "kPyJetJet","kPyJetJetFixed", + "kPyGammaJetPHOS", "kPyJetJetPHOS", "kPyJetJetPHOSv2", "kPyGammaBremsPHOS", + "kPyGammaJetEMCAL", "kPyJetJetEMCAL", "kPyGammaBremsEMCAL" +}; + + +//--- Decay Mode --- +enum DecayHvFl_t +{ + kNature, kHadr, kSemiEl, kSemiMu +}; +//--- Magnetic Field --- +enum Mag_t +{ + kNoField, k5kG +}; + +//--- Trigger config --- +enum TrigConf_t +{ + kDefaultPPTrig, kDefaultPbPbTrig +}; + +const char * TrigConfName[] = { + "p-p","Pb-Pb" +}; + +Float_t eCMS=14000; +PDC07Proc_t proc = kPyJetJetEMCAL; +Int_t year = 2009; + + +enum PprGeo_t + { + kHoles, kNoHoles + }; + +static PprGeo_t geo = kHoles; + +//--- Functions --- +AliGenPythia *PythiaHVQ(PDC07Proc_t proc); + +AliGenPythia* Blubb(Int_t proc) { + //*******************************************************************// + // Configuration file for hard QCD processes generation with PYTHIA // + // // + //*******************************************************************// + AliGenPythia * gener = 0x0; + + switch(proc) { + + case kPyJetJet: + comment = comment.Append(Form("pp->jet + jet at %3.0f TeV, pt hard %3.0f - %3.0f",eCMS,ptHardMin,ptHardMax)); + AliGenPythia * gener = new AliGenPythia(nEvts); + gener->SetEnergyCMS(eCMS);// Centre of mass energy + gener->SetProcess(kPyJets);// Process type + gener->SetJetEtaRange(-1.5, 1.5);// Final state kinematic cuts + gener->SetJetPhiRange(0., 360.); + gener->SetJetEtRange(10., 1000.); + gener->SetPtHard(ptHardMin, ptHardMax);// Pt transfer of the hard scattering + gener->SetStrucFunc(kCTEQ5L); + + return gener; + + case kPyGammaJetPHOS: + comment = comment.Append(" pp->jet + gamma over PHOS"); + gener = new AliGenPythia(nEvts); + gener->SetEnergyCMS(eCMS); + gener->SetProcess(kPyDirectGamma); + gener->SetStrucFunc(kCTEQ4L); + gener->SetPtHard(ptHardMin,ptHardMax); + //gener->SetYHard(-1.0,1.0); + gener->SetGammaEtaRange(-0.13,0.13); + gener->SetGammaPhiRange(218.,322.);//Over 5 modules +-2 deg + break; + case kPyJetJetPHOS: + comment = comment.Append(" pp->jet + jet over PHOS"); + gener = new AliGenPythia(nEvts); + gener->SetEnergyCMS(eCMS); + gener->SetProcess(kPyJets); + gener->SetStrucFunc(kCTEQ4L); + gener->SetPtHard(ptHardMin,ptHardMax); + //gener->SetYHard(-1.0,1.0); + gener->SetJetEtaRange(-1.,1.); + gener->SetJetPhiRange(200.,340.); + gener->SetPi0InPHOS(kTRUE); + gener->SetFragPhotonOrPi0MinPt(ptGammaPi0Min); + + printf("\n \n Event generator: Minimum pT of particle in calorimeter %f \n \n", ptGammaPi0Min); + break; + case kPyGammaBremsPHOS: + comment = comment.Append(" pp->jet + jet+bremsphoton over PHOS at 14 TeV"); + gener = new AliGenPythia(nEvts); + gener->SetEnergyCMS(eCMS); + gener->SetProcess(kPyJets); + gener->SetStrucFunc(kCTEQ4L); + gener->SetPtHard(ptHardMin,ptHardMax); + //gener->SetYHard(-1.0,1.0); + gener->SetJetEtaRange(-1.,1.); + gener->SetJetPhiRange(200.,340.); + gener->SetFragPhotonInPHOS(kTRUE); + gener->SetFragPhotonOrPi0MinPt(ptGammaPi0Min); + printf("\n \n Event generator: Minimum pT of particle in calorimeter %f \n \n", ptGammaPi0Min); + break; + case kPyJetJetPHOSv2: + comment = comment.Append(" pp->jet + jet over PHOS version2 "); + gener = new AliGenPythia(nEvts); + gener->SetEnergyCMS(eCMS); + gener->SetProcess(kPyJets); + gener->SetStrucFunc(kCTEQ4L); + gener->SetPtHard(ptHardMin,ptHardMax); + //gener->SetYHard(-1.0,1.0); + gener->SetJetEtaRange(-1.,1.); + gener->SetJetPhiRange(200.,340.); + //gener->SetPi0InPHOS(kTRUE); + gener->SetPhotonInPHOSeta(kTRUE); + gener->SetPhotonMinPt(ptGammaPi0Min); + gener->SetForceDecay(kAll); + break; + case kPyGammaJetEMCAL: + comment = comment.Append(" pp->jet + gamma over EMCAL at 14 TeV"); + gener = new AliGenPythia(nEvts); + gener->SetEnergyCMS(eCMS); + gener->SetProcess(kPyDirectGamma); + gener->SetStrucFunc(kCTEQ4L); + gener->SetPtHard(ptHardMin,ptHardMax); + //gener->SetYHard(-1.0,1.0); + gener->SetGammaEtaRange(-0.71,0.71); + gener->SetGammaPhiRange(78.,192.);//Over 6 supermodules +-2 deg + break; + case kPyJetJetEMCAL: + comment = comment.Append(" pp->jet + jet over EMCAL at 14 TeV"); + gener = new AliGenPythia(nEvts); + gener->SetEnergyCMS(eCMS); + gener->SetProcess(kPyJets); + gener->SetStrucFunc(kCTEQ4L); + gener->SetPtHard(ptHardMin,ptHardMax); + //gener->SetYHard(-1.0,1.0); + gener->SetJetEtaRange(-1,1); + gener->SetJetPhiRange(60.,210.); + gener->SetPi0InEMCAL(kTRUE); + gener->SetFragPhotonOrPi0MinPt(ptGammaPi0Min); + printf("\n \n Event generator: Minimum pT of particle in calorimeter %f \n \n", ptGammaPi0Min); + break; + case kPyGammaBremsEMCAL: + comment = comment.Append(" pp->jet + jet+bremsphoton over EMCAL at 14 TeV"); + gener = new AliGenPythia(nEvts); + gener->SetEnergyCMS(eCMS); + gener->SetProcess(kPyJets); + gener->SetStrucFunc(kCTEQ4L); + gener->SetPtHard(ptHardMin,ptHardMax); + //gener->SetYHard(-1.0,1.0); + gener->SetJetEtaRange(-1,1); + gener->SetJetPhiRange(60.,210.); //Over 2 uncovered PHOS modules + gener->SetFragPhotonInEMCAL(kTRUE); + gener->SetFragPhotonOrPi0MinPt(ptGammaPi0Min); + + printf("\n \n Event generator: Minimum pT of particle in calorimeter %f \n \n", ptGammaPi0Min); + break; + case kPyJetJetFixed: + comment = comment.Append(" pp->jet + jet+jet in Fixed bin for testing"); + gener = new AliGenPythia(nEvts); + gener->SetEnergyCMS(14000);// Centre of mass energy + gener->SetProcess(kPyJets);// Process type + gener->SetJetEtaRange(-1.1, 1.1);// Final state kinematic cuts + gener->SetJetPhiRange(0., 360.); + gener->SetJetEtRange(95., 105.); + gener->SetPtHard(80, 1000);// Pt transfer of the hard scattering + gener->SetStrucFunc(kCTEQ5L); + gener->SetGluonRadiation(1,1); + + + printf("\n \n Fixed Parameter \n"); + break; + + } + + return gener; +} + +AliGenerator *MbCocktail(); +AliGenerator *PyMbTriggered(Int_t pdg); +void ProcessEnvironmentVars(); + +// This part for configuration +static DecayHvFl_t decHvFl = kNature; +static Mag_t mag = k5kG; +static TrigConf_t trig = kDefaultPPTrig; // default pp trigger configuration +static Int_t runNumber= 0; +static Int_t eventNumber= 0; +//========================// +// Set Random Number seed // +//========================// +TDatime dt; +static UInt_t seed = dt.Get(); + +// nEvts = -1 : you get 1 QQbar pair and all the fragmentation and +// decay chain +// nEvts = N>0 : you get N charm / beauty Hadrons +Int_t nEvts = -1; +// stars = kTRUE : all heavy resonances and their decay stored +// = kFALSE: only final heavy hadrons and their decays stored +Bool_t stars = kTRUE; + +// To be used only with kCharmpp14000wmi and kBeautypp14000wmi +// To get a "reasonable" agreement with MNR results, events have to be +// generated with the minimum ptHard set to 2.76 GeV. +// To get a "perfect" agreement with MNR results, events have to be +// generated in four ptHard bins with the following relative +// normalizations: +// CHARM +// 2.76-3 GeV: 25% +// 3-4 GeV: 40% +// 4-8 GeV: 29% +// >8 GeV: 6% +// BEAUTY +// 2.76-4 GeV: 5% +// 4-6 GeV: 31% +// 6-8 GeV: 28% +// >8 GeV: 36% + +Float_t ptHardMin = 10.; +Float_t ptHardMax = 20.; +Float_t ptGammaPi0Min = 1.; +Int_t iquenching = 0; +Float_t qhat = 20.; + +// Comment line +static TString comment; + +void Config() +{ + + + // Get settings from environment variables + ProcessEnvironmentVars(); + +#if defined(__CINT__) + gSystem->Load("liblhapdf"); + gSystem->Load("libEGPythia6"); + gSystem->Load("libpythia6"); + gSystem->Load("libAliPythia6"); + gSystem->Load("libgeant321"); +#endif + + new TGeant3TGeo("C++ Interface to Geant3"); + + // Output every 100 tracks + ((TGeant3*)gMC)->SetSWIT(4,100); + + //======================================================================= + // Run loader + AliRunLoader* rl=0x0; + rl = AliRunLoader::Open("galice.root", + AliConfig::GetDefaultEventFolderName(), + "recreate"); + if (rl == 0x0) + { + gAlice->Fatal("Config.C","Can not instatiate the Run Loader"); + return; + } + rl->SetCompressionLevel(2); + rl->SetNumberOfEventsPerFile(1000); + gAlice->SetRunLoader(rl); + + // Run number + //gAlice->SetRunNumber(runNumber); + + // Set the trigger configuration + gAlice->SetTriggerDescriptor(TrigConfName[trig]); + cout<<"Trigger configuration is set to "<SetProcess("DCAY",1); + gMC->SetProcess("PAIR",1); + gMC->SetProcess("COMP",1); + gMC->SetProcess("PHOT",1); + gMC->SetProcess("PFIS",0); + gMC->SetProcess("DRAY",0); + gMC->SetProcess("ANNI",1); + gMC->SetProcess("BREM",1); + gMC->SetProcess("MUNU",1); + gMC->SetProcess("CKOV",1); + gMC->SetProcess("HADR",1); + gMC->SetProcess("LOSS",2); + gMC->SetProcess("MULS",1); + gMC->SetProcess("RAYL",1); + + Float_t cut = 1.e-3; // 1MeV cut by default + Float_t tofmax = 1.e10; + + gMC->SetCut("CUTGAM", cut); + gMC->SetCut("CUTELE", cut); + gMC->SetCut("CUTNEU", cut); + gMC->SetCut("CUTHAD", cut); + gMC->SetCut("CUTMUO", cut); + gMC->SetCut("BCUTE", cut); + gMC->SetCut("BCUTM", cut); + gMC->SetCut("DCUTE", cut); + gMC->SetCut("DCUTM", cut); + gMC->SetCut("PPCUTM", cut); + gMC->SetCut("TOFMAX", tofmax); + + //======================// + // Set External decayer // + //======================// + TVirtualMCDecayer* decayer = new AliDecayerPythia(); + // DECAYS + // + switch(decHvFl) { + case kNature: + decayer->SetForceDecay(kAll); + break; + case kHadr: + decayer->SetForceDecay(kHadronicD); + break; + case kSemiEl: + decayer->SetForceDecay(kSemiElectronic); + break; + case kSemiMu: + decayer->SetForceDecay(kSemiMuonic); + break; + } + decayer->Init(); + gMC->SetExternalDecayer(decayer); + if(proc == kPyJetJetPHOSv2) // in this case we decay the pi0 + decayer->SetForceDecay(kNeutralPion); + + //=========================// + // Generator Configuration // + //=========================// + AliGenPythia* gener = 0x0; + + if (proc <= kBeautypp14000wmi) { + AliGenPythia *pythia = PythiaHVQ(proc); + // FeedDown option + pythia->SetFeedDownHigherFamily(kFALSE); + // Stack filling option + if(!stars) pythia->SetStackFillOpt(AliGenPythia::kParentSelection); + // Set Count mode + if(nEvts>0) pythia->SetCountMode(AliGenPythia::kCountParents); + // + // DECAYS + // + switch(decHvFl) { + case kNature: + pythia->SetForceDecay(kAll); + break; + case kHadr: + pythia->SetForceDecay(kHadronicD); + break; + case kSemiEl: + pythia->SetForceDecay(kSemiElectronic); + break; + case kSemiMu: + pythia->SetForceDecay(kSemiMuonic); + break; + } + // + // GEOM & KINE CUTS + // + pythia->SetMomentumRange(0,99999999); + pythia->SetPhiRange(0., 360.); + pythia->SetThetaRange(0,180); + switch(ycut) { + case kFull: + pythia->SetYRange(-999,999); + break; + case kBarrel: + pythia->SetYRange(-2,2); + break; + case kMuonArm: + pythia->SetYRange(1,6); + break; + } + gener = pythia; + } else if (proc == kPyMbNoHvq) { + gener = MbCocktail(); + } else if (proc == kPyOmegaMinus) { + gener = PyMbTriggered(3334); + } else if (proc == kPyOmegaPlus) { + gener = PyMbTriggered(-3334); + } else if (proc <= kPyGammaBremsEMCAL) { + + AliGenPythia *pythia = Blubb((Int_t) proc ); + + // FeedDown option + pythia->SetFeedDownHigherFamily(kFALSE); + // Set Count mode + if(nEvts>0) pythia->SetCountMode(AliGenPythia::kCountParents); + + // + // GEOM & KINE CUTS + // + pythia->SetMomentumRange(0,99999999); + // pythia->SetJetEtaRange(-1.5, 1.5);// Final state kinematic cuts + // pythia->SetJetPhiRange(0., 360.); + // pythia->SetThetaRange(45,135); + + if(proc == kPyJetJetPHOSv2) + pythia->SetForceDecay(kNeutralPion); + else + pythia->SetForceDecay(kAll); + + pythia->SetPycellParameters(2., 274, 432, 0., 4., 5., 1.0); + pythia->SetPtKick(3.); // set the intrinsic kt to 5 GeV/c + pythia->SetGluonRadiation(1,1); + gener = pythia; + } + + + // PRIMARY VERTEX + + gener->SetOrigin(0., 0., 0.); // vertex position + + // Size of the interaction diamond + // Longitudinal + Float_t sigmaz = 7.55 / TMath::Sqrt(2.); // [cm] + + // Transverse + Float_t betast = 10; // beta* [m] + Float_t eps = 3.75e-6; // emittance [m] + Float_t gamma = 7000. / 0.938272; // relativistic gamma [1] + Float_t sigmaxy = TMath::Sqrt(eps * betast / gamma) / TMath::Sqrt(2.) * 100.; // [cm] + printf("\n \n Diamond size x-y: %10.3e z: %10.3e\n \n", sigmaxy, sigmaz); + + gener->SetSigma(sigmaxy, sigmaxy, sigmaz); // Sigma in (X,Y,Z) (cm) on IP position + gener->SetCutVertexZ(3.); // Truncate at 3 sigma + gener->SetVertexSmear(kPerEvent); + + gener->Init(); + + //Quenching + gener->SetQuench(iquenching); + if(iquenching == 1){ + Float_t k = 6e5*(qhat/1.7) ; //qhat=1.7, k = 6e5, default value + AliPythia::Instance()->InitQuenching(0.,0.1,k,0,0.95,6); + + } + // FIELD + + printf("\n \n Comment: %s \n \n", comment.Data()); + if (mag == kNoField) { + comment = comment.Append(" | L3 field 0.0 T"); + field = new AliMagWrapCheb("Maps","Maps", 2, 0., 10., AliMagWrapCheb::k2kG); + } else if (mag == k5kG) { + comment = comment.Append(" | L3 field 0.5 T"); + // field = new AliMagWrapCheb("Maps","Maps", 2, 1., 10., AliMagWrapCheb::k5kG); + field = new AliMagWrapCheb("Maps","Maps", 2, 1., 10., AliMagWrapCheb::k5kG,kTRUE,"$(ALICE_ROOT)/data/maps/mfchebKGI_sym.root"); + } + + printf("\n \n Comment: %s \n \n", comment.Data()); + rl->CdGAFile(); + gAlice->SetField(field); + + + + Int_t iABSO = 1; + Int_t iACORDE = 0; + Int_t iDIPO = 1; + Int_t iEMCAL = 1; + Int_t iFMD = 1; + Int_t iFRAME = 1; + Int_t iHALL = 1; + Int_t iITS = 1; + Int_t iMAG = 1; + Int_t iMUON = 1; + Int_t iPHOS = 1; + Int_t iPIPE = 1; + Int_t iPMD = 0 ; + Int_t iHMPID = 1; + Int_t iSHIL = 1; + Int_t iT0 = 1; + Int_t iTOF = 1; + Int_t iTPC = 1; + Int_t iTRD = 1; + Int_t iVZERO = 1; + Int_t iZDC = 1; + + + //=================== Alice BODY parameters ============================= + AliBODY *BODY = new AliBODY("BODY", "Alice envelop"); + + + if (iMAG) + { + //=================== MAG parameters ============================ + // --- Start with Magnet since detector layouts may be depending --- + // --- on the selected Magnet dimensions --- + AliMAG *MAG = new AliMAG("MAG", "Magnet"); + } + + + if (iABSO) + { + //=================== ABSO parameters ============================ + AliABSO *ABSO = new AliABSOv3("ABSO", "Muon Absorber"); + } + + if (iDIPO) + { + //=================== DIPO parameters ============================ + + AliDIPO *DIPO = new AliDIPOv3("DIPO", "Dipole version 3"); + } + + if (iHALL) + { + //=================== HALL parameters ============================ + + AliHALL *HALL = new AliHALLv3("HALL", "Alice Hall"); + } + + + if (iFRAME) + { + //=================== FRAME parameters ============================ + + AliFRAMEv2 *FRAME = new AliFRAMEv2("FRAME", "Space Frame"); + if (geo == kHoles) FRAME->SetHoles(1); + else FRAME->SetHoles(0); + + } + + if (iSHIL) + { + //=================== SHIL parameters ============================ + + AliSHIL *SHIL = new AliSHILv3("SHIL", "Shielding Version 3"); + } + + + if (iPIPE) + { + //=================== PIPE parameters ============================ + + AliPIPE *PIPE = new AliPIPEv3("PIPE", "Beam Pipe"); + } + + if(iITS) { + + //=================== ITS parameters ============================ + // + AliITS *ITS = new AliITSv11Hybrid("ITS","ITS v11Hybrid"); + // AliITSvPPRasymmFMD *ITS = new AliITSvPPRasymmFMD("ITS","New ITS PPR detailed version with asymmetric services"); + } + + if (iTPC) + { + //============================ TPC parameters ===================== + AliTPC *TPC = new AliTPCv2("TPC", "Default"); + } + + + if (iTOF) { + //=================== TOF parameters ============================ + AliTOF *TOF = new AliTOFv6T0("TOF", "normal TOF"); + } + + + if (iHMPID) + { + //=================== HMPID parameters =========================== + // AliHMPID *HMPID = new AliHMPIDv2("HMPID", "normal HMPID"); + AliHMPID *HMPID = new AliHMPIDv3("HMPID", "normal HMPID"); + + } + + + if (iZDC) + { + //=================== ZDC parameters ============================ + + AliZDC *ZDC = new AliZDCv3("ZDC", "normal ZDC"); + } + + if (iTRD) + { + //=================== TRD parameters ============================ + + AliTRD *TRD = new AliTRDv1("TRD", "TRD slow simulator"); + AliTRDgeometry *geoTRD = TRD->GetGeometry(); + // Partial geometry: modules at 0,8,9,17 + // Partial geometry: modules at 1,7,10,16 expected for 2009 + // starting at 3h in positive direction + + + geoTRD->SetSMstatus(2,0); + geoTRD->SetSMstatus(3,0); + geoTRD->SetSMstatus(4,0); + geoTRD->SetSMstatus(5,0); + geoTRD->SetSMstatus(6,0); + geoTRD->SetSMstatus(11,0); + geoTRD->SetSMstatus(12,0); + geoTRD->SetSMstatus(13,0); + geoTRD->SetSMstatus(14,0); + geoTRD->SetSMstatus(15,0); + + + + + /* + // Partial geometry: modules at 2,3,4,6,11,12,14,15 + // starting at 6h in positive direction + //Hole on SM 13 and 14 for PHOS + geoTRD->SetSMstatus(13,0); + geoTRD->SetSMstatus(14,0); + geoTRD->SetSMstatus(15,0); + */ + + } + + if (iFMD) + { + //=================== FMD parameters ============================ + AliFMD *FMD = new AliFMDv1("FMD", "normal FMD"); + } + + if (iMUON) + { + //=================== MUON parameters =========================== + // New MUONv1 version (geometry defined via builders) + AliMUON *MUON = new AliMUONv1("MUON", "default"); + } + //=================== PHOS parameters =========================== + + if (iPHOS) + { + AliPHOS *PHOS = new AliPHOSv1("PHOS", "IHEP"); + // Default for cold phos... + } + + + if (iPMD) + { + //=================== PMD parameters ============================ + AliPMD *PMD = new AliPMDv1("PMD", "normal PMD"); + } + + if (iT0) + { + //=================== T0 parameters ============================ + AliT0 *T0 = new AliT0v1("T0", "T0 Detector"); + } + + if (iEMCAL) + { + //=================== EMCAL parameters ============================ + AliEMCAL *EMCAL = new AliEMCALv2("EMCAL", "EMCAL_COMPLETE"); + } + + if (iACORDE) + { + //=================== CRT parameters ============================ + AliACORDE *ACORDE = new AliACORDEv1("CRT", "normal ACORDE"); + } + + if (iVZERO) + { + //=================== CRT parameters ============================ + AliVZERO *VZERO = new AliVZEROv7("VZERO", "normal VZERO"); + } +} + +// PYTHIA + +AliGenPythia *PythiaHVQ(PDC07Proc_t proc) { +//*******************************************************************// +// Configuration file for charm / beauty generation with PYTHIA // +// // +// The parameters have been tuned in order to reproduce the inclusive// +// heavy quark pt distribution given by the NLO pQCD calculation by // +// Mangano, Nason and Ridolfi. // +// // +// For details and for the NORMALIZATION of the yields see: // +// N.Carrer and A.Dainese, // +// "Charm and beauty production at the LHC", // +// ALICE-INT-2003-019, [arXiv:hep-ph/0311225]; // +// PPR Chapter 6.6, CERN/LHCC 2005-030 (2005). // +//*******************************************************************// + AliGenPythia * gener = 0x0; + + switch(proc) { + case kCharmPbPb5500: + comment = comment.Append(" Charm in Pb-Pb at 5.5 TeV"); + gener = new AliGenPythia(nEvts); + gener->SetProcess(kPyCharmPbPbMNR); + gener->SetStrucFunc(kCTEQ4L); + gener->SetPtHard(2.1,-1.0); + gener->SetEnergyCMS(5500.); + gener->SetNuclei(208,208); + break; + case kCharmpPb8800: + comment = comment.Append(" Charm in p-Pb at 8.8 TeV"); + gener = new AliGenPythia(nEvts); + gener->SetProcess(kPyCharmpPbMNR); + gener->SetStrucFunc(kCTEQ4L); + gener->SetPtHard(2.1,-1.0); + gener->SetEnergyCMS(8800.); + gener->SetProjectile("P",1,1); + gener->SetTarget("Pb",208,82); + break; + case kCharmpp14000: + comment = comment.Append(" Charm in pp at 14 TeV"); + gener = new AliGenPythia(nEvts); + gener->SetProcess(kPyCharmppMNR); + gener->SetStrucFunc(kCTEQ4L); + gener->SetPtHard(2.1,-1.0); + gener->SetEnergyCMS(14000.); + break; + case kCharmpp14000wmi: + comment = comment.Append(" Charm in pp at 14 TeV with mult. interactions"); + gener = new AliGenPythia(-1); + gener->SetProcess(kPyCharmppMNRwmi); + gener->SetStrucFunc(kCTEQ5L); + gener->SetPtHard(ptHardMin,ptHardMax); + gener->SetEnergyCMS(14000.); + break; + case kD0PbPb5500: + comment = comment.Append(" D0 in Pb-Pb at 5.5 TeV"); + gener = new AliGenPythia(nEvts); + gener->SetProcess(kPyD0PbPbMNR); + gener->SetStrucFunc(kCTEQ4L); + gener->SetPtHard(2.1,-1.0); + gener->SetEnergyCMS(5500.); + gener->SetNuclei(208,208); + break; + case kD0pPb8800: + comment = comment.Append(" D0 in p-Pb at 8.8 TeV"); + gener = new AliGenPythia(nEvts); + gener->SetProcess(kPyD0pPbMNR); + gener->SetStrucFunc(kCTEQ4L); + gener->SetPtHard(2.1,-1.0); + gener->SetEnergyCMS(8800.); + gener->SetProjectile("P",1,1); + gener->SetTarget("Pb",208,82); + break; + case kD0pp14000: + comment = comment.Append(" D0 in pp at 14 TeV"); + gener = new AliGenPythia(nEvts); + gener->SetProcess(kPyD0ppMNR); + gener->SetStrucFunc(kCTEQ4L); + gener->SetPtHard(2.1,-1.0); + gener->SetEnergyCMS(14000.); + break; + case kDPlusPbPb5500: + comment = comment.Append(" DPlus in Pb-Pb at 5.5 TeV"); + gener = new AliGenPythia(nEvts); + gener->SetProcess(kPyDPlusPbPbMNR); + gener->SetStrucFunc(kCTEQ4L); + gener->SetPtHard(2.1,-1.0); + gener->SetEnergyCMS(5500.); + gener->SetNuclei(208,208); + break; + case kDPluspPb8800: + comment = comment.Append(" DPlus in p-Pb at 8.8 TeV"); + gener = new AliGenPythia(nEvts); + gener->SetProcess(kPyDPluspPbMNR); + gener->SetStrucFunc(kCTEQ4L); + gener->SetPtHard(2.1,-1.0); + gener->SetEnergyCMS(8800.); + gener->SetProjectile("P",1,1); + gener->SetTarget("Pb",208,82); + break; + case kDPluspp14000: + comment = comment.Append(" DPlus in pp at 14 TeV"); + gener = new AliGenPythia(nEvts); + gener->SetProcess(kPyDPlusppMNR); + gener->SetStrucFunc(kCTEQ4L); + gener->SetPtHard(2.1,-1.0); + gener->SetEnergyCMS(14000.); + break; + case kBeautyPbPb5500: + comment = comment.Append(" Beauty in Pb-Pb at 5.5 TeV"); + gener = new AliGenPythia(nEvts); + gener->SetProcess(kPyBeautyPbPbMNR); + gener->SetStrucFunc(kCTEQ4L); + gener->SetPtHard(2.75,-1.0); + gener->SetEnergyCMS(5500.); + gener->SetNuclei(208,208); + break; + case kBeautypPb8800: + comment = comment.Append(" Beauty in p-Pb at 8.8 TeV"); + gener = new AliGenPythia(nEvts); + gener->SetProcess(kPyBeautypPbMNR); + gener->SetStrucFunc(kCTEQ4L); + gener->SetPtHard(2.75,-1.0); + gener->SetEnergyCMS(8800.); + gener->SetProjectile("P",1,1); + gener->SetTarget("Pb",208,82); + break; + case kBeautypp14000: + comment = comment.Append(" Beauty in pp at 14 TeV"); + gener = new AliGenPythia(nEvts); + gener->SetProcess(kPyBeautyppMNR); + gener->SetStrucFunc(kCTEQ4L); + gener->SetPtHard(2.75,-1.0); + gener->SetEnergyCMS(14000.); + break; + case kBeautypp14000wmi: + comment = comment.Append(" Beauty in pp at 14 TeV with mult. interactions"); + gener = new AliGenPythia(-1); + gener->SetProcess(kPyBeautyppMNRwmi); + gener->SetStrucFunc(kCTEQ5L); + gener->SetPtHard(ptHardMin,ptHardMax); + gener->SetEnergyCMS(14000.); + break; + } + + return gener; +} + +AliGenerator* MbCocktail() +{ + comment = comment.Append(" pp at 14 TeV: Pythia low-pt, no heavy quarks + J/Psi from parameterisation"); + AliGenCocktail * gener = new AliGenCocktail(); + gener->UsePerEventRates(); + +// Pythia + + AliGenPythia* pythia = new AliGenPythia(-1); + pythia->SetMomentumRange(0, 999999.); + pythia->SetThetaRange(0., 180.); + pythia->SetYRange(-12.,12.); + pythia->SetPtRange(0,1000.); + pythia->SetProcess(kPyMb); + pythia->SetEnergyCMS(14000.); + pythia->SwitchHFOff(); + +// J/Psi parameterisation + + AliGenParam* jpsi = new AliGenParam(1, AliGenMUONlib::kJpsi, "CDF scaled", "Jpsi"); + jpsi->SetPtRange(0.,100.); + jpsi->SetYRange(-8., 8.); + jpsi->SetPhiRange(0., 360.); + jpsi->SetForceDecay(kAll); + + gener->AddGenerator(pythia, "Pythia", 1.); + gener->AddGenerator(jpsi, "J/Psi", 8.e-4); + + return gener; +} + +AliGenerator* PyMbTriggered(Int_t pdg) +{ + AliGenPythia* pythia = new AliGenPythia(-1); + pythia->SetMomentumRange(0, 999999.); + pythia->SetThetaRange(0., 180.); + pythia->SetYRange(-12.,12.); + pythia->SetPtRange(0,1000.); + pythia->SetProcess(kPyMb); + pythia->SetEnergyCMS(14000.); + pythia->SetTriggerParticle(pdg, 0.9); + return pythia; +} + +void ProcessEnvironmentVars() +{ + cout << "######################################" << endl; + cout << "## Processing environment variables ##" << endl; + cout << "######################################" << endl; + + // Run Number + if (gSystem->Getenv("DC_RUN")) { + runNumber = atoi(gSystem->Getenv("DC_RUN")); + } + //cout<<"Run number "<Getenv("DC_EVENT")) { + eventNumber = atoi(gSystem->Getenv("DC_EVENT")); + } + //cout<<"Event number "<Getenv("CONFIG_SEED")) { + seed = atoi(gSystem->Getenv("CONFIG_SEED")); + } + else if(gSystem->Getenv("DC_EVENT") && gSystem->Getenv("DC_RUN")){ + seed = runNumber * 100000 + eventNumber; + } + + gRandom->SetSeed(seed); + + cout<<"////////////////////////////////////////////////////////////////////////////////////"<GetSeed()<Getenv("DC_RUN_TYPE")) { + cout<<"run type "<Getenv("DC_RUN_TYPE")<Getenv("DC_RUN_TYPE"), pprRunName[iRun])==0) { + proc = (PDC07Proc_t)iRun; + } + } + } + else + cout << "Environment variable DC_RUN_TYPE is not defined" << endl; + + if (gSystem->Getenv("ECMS")) + eCMS = atof(gSystem->Getenv("ECMS")); + if (gSystem->Getenv("PTHARDMIN")) + ptHardMin = atof(gSystem->Getenv("PTHARDMIN")); + if (gSystem->Getenv("PTHARDMAX")) + ptHardMax = atof(gSystem->Getenv("PTHARDMAX")); + if (gSystem->Getenv("PTGAMMAPI0MIN")) + ptGammaPi0Min = atof(gSystem->Getenv("PTGAMMAPI0MIN")); + if (gSystem->Getenv("QUENCHING")) + iquenching = atof(gSystem->Getenv("QUENCHING")); + if (gSystem->Getenv("QHAT")) + qhat = atof(gSystem->Getenv("QHAT")); + + cout<<">> Run type set to "<> Collision energy set to "<> ptHard limits: "<> pt gamma/pi0 threshold "<< ptGammaPi0Min<<" GeV "<> quenching on? "<< iquenching<<" qhat "<SetComment("AOD Reader"); + jrh->SetPtCut(0.); + jrh->SetTestFilterMask(1<<1); // hard cuts... + // Define reader and set its header + AliJetAODReader *er = new AliJetAODReader(); + er->SetReaderHeader(jrh); + + // Define jet header + AliUA1JetHeaderV1 *jh=new AliUA1JetHeaderV1(); + jh->SetComment("UA1 jet code with default parameters"); + jh->BackgMode(0); + jh->SetRadius(0.4); + jh->SetEtSeed(2.); + jh->SetLegoNbinPhi(420.); + jh->SetLegoNbinEta(120.); + jh->SetLegoEtaMin(-1.9); + jh->SetLegoEtaMax(+1.9); + jh->SetMinJetEt(5.); + jh->SetJetEtaMin(-1.5); + jh->SetJetEtaMax(1.5); + + + + // Define jet finder. Set its header and reader + jetFinder = new AliUA1JetFinderV1(); + jetFinder->SetJetHeader(jh); + jetFinder->SetJetReader(er); + /* + jetFinder->SetPlotMode(kTRUE); + jetFinder->SetOutputFile("jets.root"); + */ + // + return jetFinder; +} diff --git a/prod/LHC09a1/ConfigJetAnalysisMC.C b/prod/LHC09a1/ConfigJetAnalysisMC.C new file mode 100644 index 00000000000..585d021e91d --- /dev/null +++ b/prod/LHC09a1/ConfigJetAnalysisMC.C @@ -0,0 +1,39 @@ +AliJetFinder* ConfigJetAnalysis() +{ + // + // Configuration goes here + // + printf("ConfigJetAnalysis() \n"); + Printf("For MC events"); + AliJetKineReaderHeader *jrh = new AliJetKineReaderHeader(); + jrh->SetComment("MC full Kinematics"); + jrh->SetFastSimTPC(kFALSE); + jrh->SetFastSimEMCAL(kFALSE); + jrh->SetFiducialEta(-2.1,2.1); + jrh->SetPtCut(0.); + + // Define reader and set its header + AliJetKineReader *er = new AliJetKineReader(); + er->SetReaderHeader(jrh); + + + // Define jet header + AliUA1JetHeaderV1 *jh=new AliUA1JetHeaderV1(); + jh->SetComment("UA1 jet code with MC parameters"); + jh->BackgMode(0); + jh->SetRadius(1.0); + jh->SetEtSeed(4.); + jh->SetLegoNbinPhi(432.); + jh->SetLegoNbinEta(274.); + jh->SetLegoEtaMin(-2); + jh->SetLegoEtaMax(+2); + jh->SetJetEtaMin(-1.5); + jh->SetJetEtaMax(1.5); + jh->SetMinJetEt(5.); + + // Define jet finder. Set its header and reader + jetFinder = new AliUA1JetFinderV1(); + jetFinder->SetJetHeader(jh); + jetFinder->SetJetReader(er); + return jetFinder; +} diff --git a/prod/LHC09a1/JDL b/prod/LHC09a1/JDL new file mode 100644 index 00000000000..4b4be8b56c1 --- /dev/null +++ b/prod/LHC09a1/JDL @@ -0,0 +1,38 @@ +Executable = "aliroot_new"; + +Jobtag={"comment:pp, jet-jet events 15 < pt hard < 50 GeV, ID #76"}; + +# sequence is really essential here!! +Packages={"VO_ALICE@AliRoot::v4-16-Rev-09","VO_ALICE@GEANT3::v1-9-8","VO_ALICE@ROOT::v5-23-02","VO_ALICE@APISCONFIG::V2.4"}; +TTL = "72000"; +PRICE = 1; + +ValidationCommand="/alice/cern.ch/user/a/aliprod/LHC09a1/validation.sh"; + +InputFile= {"LF:/alice/cern.ch/user/a/aliprod/LHC09a1/Config.C", + "LF:/alice/cern.ch/user/a/aliprod/LHC09a1/simrun.C", + "LF:/alice/cern.ch/user/a/aliprod/LHC09a1/sim.C", + "LF:/alice/cern.ch/user/a/aliprod/LHC09a1/rec.C", + "LF:/alice/cern.ch/user/a/aliprod/LHC09a1/rec2.C", + "LF:/alice/cern.ch/user/a/aliprod/LHC09a1/tag.C", + "LF:/alice/cern.ch/user/a/aliprod/LHC09a1/ConfigJetAnalysis.C", + "LF:/alice/cern.ch/user/a/aliprod/LHC09a1/ConfigJetAnalysisMC.C", + "LF:/alice/cern.ch/user/a/aliprod/LHC09a1/runAODFilterJets.C", + "LF:/alice/cern.ch/user/a/aliprod/LHC09a1/CheckESD.C"}; + +OutputArchive={"log_archive:stdout,stderr,*.log@ALICE::Grenoble::DPM","root_archive.zip:galice.root,geometry.root,Kinematics.root,TrackRefs.root,AliESDs*.root,AliESDfriends*.root,AliAOD*.root,*tag.root,*QA*.root,pyxsec.root@ALICE::FZK::SE,ALICE::Torino::DPM,ALICE::IPNO::DPM"}; + +OutputFile={"Run*.root@ALICE::FZK::SE,ALICE::Torino::DPM,ALICE::IPNO::DPM"}; + +OutputDir="/alice/sim/PDC_08b/LHC09a1/$1/#alien_counter_03i#"; + +JDLVariables={"Packages", "OutputDir"}; +GUIDFILE="guid.txt"; + + +# pthardbin 0: 15- 50 1: 50-100 2: 100 - inf defined in simrun.C +# change for diferent procution cycles!! +splitarguments="simrun.C --run $1 --event 150 --aliencounter #alien_counter# --process kPyJetJet --pthardbin 0 --quench 0 --qhat 0" ; +split="production:1-1000"; + +Workdirectorysize={"4000MB"}; diff --git a/prod/LHC09a1/rec.C b/prod/LHC09a1/rec.C new file mode 100644 index 00000000000..094905c7181 --- /dev/null +++ b/prod/LHC09a1/rec.C @@ -0,0 +1,37 @@ +void rec() { + + AliReconstruction reco; + reco.SetWriteESDfriend(); + reco.SetWriteAlignmentData(); + + reco.SetRecoParam("ITS",AliITSRecoParam::GetLowFluxParam()); + reco.SetRecoParam("TPC",AliTPCRecoParam::GetLowFluxParam()); + reco.SetRecoParam("TRD",AliTRDrecoParam::GetLowFluxParam()); + reco.SetRecoParam("PHOS",AliPHOSRecoParam::GetDefaultParameters()); + reco.SetRecoParam("MUON",AliMUONRecoParam::GetLowFluxParam()); + reco.SetRecoParam("EMCAL",AliEMCALRecParam::GetLowFluxParam()); + + reco.SetDefaultStorage("alien://Folder=/alice/simulation/2008/v4-15-Release/Residual/"); + + + /* + reco.SetSpecificStorage("GRP/GRP/Data", + "alien://Folder=/alice/simulation/2008/v4-15-Release/Ideal/"); + */ + // No write access to the OCDB => local specific storage + reco.SetSpecificStorage("GRP/GRP/Data", + Form("local://%s/../",gSystem->pwd())); + + + // add TRD standalone tracks + reco.SetOption("TRD", "sl_tr_1"); + + reco.SetRunQA("ALL:ALL"); + + + TStopwatch timer; + timer.Start(); + reco.Run(); + timer.Stop(); + timer.Print(); +} diff --git a/prod/LHC09a1/rec2.C b/prod/LHC09a1/rec2.C new file mode 100644 index 00000000000..388f1ec4728 --- /dev/null +++ b/prod/LHC09a1/rec2.C @@ -0,0 +1,36 @@ +void rec2() { + + + AliReconstruction reco; + reco.SetWriteESDfriend(); + reco.SetWriteAlignmentData(); + + reco.SetRecoParam("ITS",AliITSRecoParam::GetLowFluxParam()); + reco.SetRecoParam("TPC",AliTPCRecoParam::GetLowFluxParam()); + reco.SetRecoParam("TRD",AliTRDrecoParam::GetLowFluxParam()); + reco.SetRecoParam("PHOS",AliPHOSRecoParam::GetDefaultParameters()); + reco.SetRecoParam("MUON",AliMUONRecoParam::GetLowFluxParam()); + reco.SetRecoParam("EMCAL",AliEMCALRecParam::GetLowFluxParam()); + + reco.SetDefaultStorage("alien://Folder=/alice/simulation/2008/v4-15-Release/Residual/"); + + /* + reco.SetSpecificStorage("GRP/GRP/Data", + "alien://Folder=/alice/simulation/2008/v4-15-Release/Ideal/"); + */ + // No write access to the OCDB => local specific storage + reco.SetSpecificStorage("GRP/GRP/Data", + Form("local://%s/../",gSystem->pwd())); + + // TRD standalone tracks and reco + reco.SetOption("TRD", "sl_tr_1"); + reco.SetRunTracking("TRD"); + + reco.SetRunQA("ALL:ALL"); + + TStopwatch timer; + timer.Start(); + reco.Run(); + timer.Stop(); + timer.Print(); +} diff --git a/prod/LHC09a1/runAODFilterJets.C b/prod/LHC09a1/runAODFilterJets.C new file mode 100644 index 00000000000..a3208290be5 --- /dev/null +++ b/prod/LHC09a1/runAODFilterJets.C @@ -0,0 +1,143 @@ +void runAODFilterJets() +{ + + bool bKineFilter = true; + bool bJets = true; + + gSystem->Load("libTree.so"); + gSystem->Load("libPhysics.so"); + gSystem->Load("libGeom.so"); + gSystem->Load("libVMC.so"); + + gSystem->Load("libSTEERBase.so"); + gSystem->Load("libESD.so"); + gSystem->Load("libAOD.so"); + gSystem->Load("libANALYSIS.so"); + gSystem->Load("libANALYSISalice.so"); + if(bJets)gSystem->Load("libJETAN.so"); + + // Create the chain + // + TChain *chain = new TChain("esdTree"); + chain->Add("AliESDs.root"); + + //////////////////////////////////////////////////////////////////////// + // Create the analysis manager + // + // Input + AliESDInputHandler* inpHandler = new AliESDInputHandler(); + // Output + AliAODHandler* aodHandler = new AliAODHandler(); + aodHandler->SetOutputFileName("AliAOD.root"); + // MC Truth + AliMCEventHandler* mcHandler = new AliMCEventHandler(); + AliAnalysisManager *mgr = new AliAnalysisManager("Filter Manager", "Filter Manager"); + if(bKineFilter){ + mgr->SetMCtruthEventHandler(mcHandler); + } + + mgr->SetInputEventHandler (inpHandler); + mgr->SetOutputEventHandler (aodHandler); + + mgr->SetDebugLevel(10); + + // Filtering of MC particles (decays conversions etc) + // this task is also needed to set the MCEventHandler + // to the AODHandler, this will not be needed when + // AODHandler goes to ANALYSISalice + AliAnalysisTaskMCParticleFilter *kinefilter = new AliAnalysisTaskMCParticleFilter("Particle Filter"); + if(bKineFilter)mgr->AddTask(kinefilter); + + + // + // soft + AliESDtrackCuts* esdTrackCutsL = new AliESDtrackCuts("AliESDtrackCuts", "Loose"); + esdTrackCutsL->SetMinNClustersTPC(50); + esdTrackCutsL->SetMaxChi2PerClusterTPC(3.5); + esdTrackCutsL->SetMaxCovDiagonalElements(2,2,0.5,0.5,2); + esdTrackCutsL->SetRequireTPCRefit(kTRUE); + esdTrackCutsL->SetMaxDCAToVertexZ(3.0); + esdTrackCutsL->SetMaxDCAToVertexXY(3.0); + esdTrackCutsL->SetRequireSigmaToVertex(kFALSE); + esdTrackCutsL->SetAcceptKingDaughters(kFALSE); + + // + // hard + AliESDtrackCuts* esdTrackCutsH = new AliESDtrackCuts("AliESDtrackCuts", "Har +d"); + esdTrackCutsH->SetMinNClustersTPC(100); + esdTrackCutsH->SetMaxChi2PerClusterTPC(2.0); + esdTrackCutsH->SetMaxCovDiagonalElements(2,2,0.5,0.5,2); + esdTrackCutsH->SetRequireTPCRefit(kTRUE); + esdTrackCutsL->SetMaxDCAToVertexZ(3.0); + esdTrackCutsL->SetMaxDCAToVertexXY(3.0); + esdTrackCutsH->SetRequireSigmaToVertex(kFALSE); + esdTrackCutsH->SetAcceptKingDaughters(kFALSE); + // + + AliAnalysisTaskJets *jetana = 0; + AliAnalysisTaskJets *jetanaMC = 0; + + if(bJets){ + jetana = new AliAnalysisTaskJets("JetAnalysis"); + jetana->SetConfigFile("../ConfigJetAnalysis.C"); + jetana->SetDebugLevel(10); + + jetanaMC = new AliAnalysisTaskJets("JetAnalysisMC"); + jetanaMC->SetDebugLevel(10); + jetanaMC->SetConfigFile("../ConfigJetAnalysisMC.C"); + jetanaMC->SetNonStdBranch("jetsMC"); + mgr->AddTask(jetanaMC); + mgr->AddTask(jetana); + } + + AliAnalysisFilter* trackFilter = new AliAnalysisFilter("trackFilter"); + trackFilter->AddCuts(esdTrackCutsL); + trackFilter->AddCuts(esdTrackCutsH); + + AliAnalysisTaskESDfilter *esdfilter = new AliAnalysisTaskESDfilter("ESD Filter"); + esdfilter->SetTrackFilter(trackFilter); + + mgr->AddTask(esdfilter); + + + // + // Create containers for input/output + AliAnalysisDataContainer *cinput1 = mgr->CreateContainer("cchain",TChain::Class(), + AliAnalysisManager::kInputContainer); + + AliAnalysisDataContainer *coutput1 = mgr->CreateContainer("tree", TTree::Class(), + AliAnalysisManager::kOutputContainer, "default"); + + coutput1->SetSpecialOutput(); + + if(bKineFilter){ + mgr->ConnectInput (kinefilter, 0, cinput1 ); + mgr->ConnectOutput (kinefilter, 0, coutput1 ); + } + + mgr->ConnectInput (esdfilter, 0, cinput1 ); + mgr->ConnectOutput (esdfilter, 0, coutput1 ); + + if(bJets){ + AliAnalysisDataContainer *coutput2 = mgr->CreateContainer("histos", TList::Class(), + AliAnalysisManager::kOutputContainer, "histosJets.root"); + + AliAnalysisDataContainer *coutputMC2 = mgr->CreateContainer("histosMC", TList::Class(), AliAnalysisManager::kOutputContainer, "histosJetsMC.root"); + + mgr->ConnectInput (jetana, 0, cinput1 ); + mgr->ConnectOutput (jetana, 0, coutput1 ); + mgr->ConnectOutput (jetana, 1, coutput2 ); + mgr->ConnectInput (jetanaMC, 0, cinput1 ); + mgr->ConnectOutput (jetanaMC, 0, coutput1 ); + mgr->ConnectOutput (jetanaMC, 1, coutputMC2 ); + } + + // + // Run the analysis + // + mgr->InitAnalysis(); + mgr->PrintStatus(); + mgr->StartAnalysis("local",chain); + +} diff --git a/prod/LHC09a1/sim.C b/prod/LHC09a1/sim.C new file mode 100644 index 00000000000..cbd9396563c --- /dev/null +++ b/prod/LHC09a1/sim.C @@ -0,0 +1,43 @@ +void sim(Int_t nev=150) { + AliSimulation simu; + simu.SetMakeSDigits("TRD TOF PHOS HMPID EMCAL MUON FMD ZDC T0 VZERO"); + simu.SetMakeDigitsFromHits("ITS TPC"); + + simu.SetDefaultStorage("alien://Folder=/alice/simulation/2008/v4-15-Release/Ideal/"); + + // No write access to the OCDB => local specific storage + simu.SetSpecificStorage("GRP/GRP/Data", + Form("local://%s/",gSystem->pwd())); + + Printf("%s:%d %d events",(char*)__FILE__,__LINE__,nev); + + TStopwatch timer; + timer.Start(); + simu.Run(nev); + WriteXsection(); + timer.Stop(); + timer.Print(); +} + +WriteXsection() +{ + TPythia6 *pythia = TPythia6::Instance(); + pythia->Pystat(1); + Double_t xsection = pythia->GetPARI(1); + UInt_t ntrials = pythia->GetMSTI(5); + + TFile *file = new TFile("pyxsec.root","recreate"); + TTree *tree = new TTree("Xsection","Pythia cross section"); + TBranch *branch = tree->Branch("xsection", &xsection, "X/D"); + TBranch *branch = tree->Branch("ntrials" , &ntrials , "X/i"); + tree->Fill(); + + + + tree->Write(); + file->Close(); + + cout << "Pythia cross section: " << xsection + << ", number of trials: " << ntrials << endl; +} + diff --git a/prod/LHC09a1/simrun.C b/prod/LHC09a1/simrun.C new file mode 100644 index 00000000000..9edc7b34899 --- /dev/null +++ b/prod/LHC09a1/simrun.C @@ -0,0 +1,234 @@ +//#define VERBOSEARGS +{ + // set job and simulation variables as : alienoutner needed for seed + // root run.C --run --event --aliencounter --process --minhard --maxhard --minpt + + + char* program = "aliroot"; + // char* program = "alienaliroot"; + + int nrun = 0; + int nevent = 0; + int seed = 0; + int naliencounter = 0; + + char sseed[1024]; + char srun[1024]; + char sevent[1024]; + char saliencounter[1024]; + char sprocess[1024]; + char sminpthard[1024]; + char smaxpthard[1024]; + char sminptgammapi0[1024]; + char squench[1024]; + char sqhat[1024]; + + Int_t npthardbin = 0; + const Int_t nPtHardBins = 3; + const Int_t ptHardLo[nPtHardBins] = {15,50,100}; + const Int_t ptHardHi[nPtHardBins] = {50,100,-1}; + + + + + sprintf(srun,""); + sprintf(sevent,""); + sprintf(saliencounter,""); + sprintf(sprocess,""); + sprintf(sminpthard,""); + sprintf(smaxpthard,""); + sprintf(sminptgammapi0,""); + sprintf(squench,""); + sprintf(sqhat,""); + + for (int i=0; i< gApplication->Argc();i++){ +#ifdef VERBOSEARGS + printf("Arg %d: %s\n",i,gApplication->Argv(i)); +#endif + if (!(strcmp(gApplication->Argv(i),"--run"))) + nrun = atoi(gApplication->Argv(i+1)); + sprintf(srun,"%d",nrun); + if (!(strcmp(gApplication->Argv(i),"--event"))) + nevent = atoi(gApplication->Argv(i+1)); + sprintf(sevent,"%d",nevent); + if (!(strcmp(gApplication->Argv(i),"--aliencounter"))) + naliencounter = atoi(gApplication->Argv(i+1)); + sprintf(saliencounter,"%d",naliencounter); + + if (!(strcmp(gApplication->Argv(i),"--process"))) + sprintf(sprocess, gApplication->Argv(i+1)); + + if (!(strcmp(gApplication->Argv(i),"--pthardbin"))) + npthardbin = atoi(gApplication->Argv(i+1)); + + + if (!(strcmp(gApplication->Argv(i),"--minhard"))) + sprintf(sminpthard,gApplication->Argv(i+1)); + + if (!(strcmp(gApplication->Argv(i),"--maxhard"))) + sprintf(smaxpthard,gApplication->Argv(i+1)); + + + if (!(strcmp(gApplication->Argv(i),"--minpt"))) + sprintf(sminptgammapi0,gApplication->Argv(i+1)); + + if (!(strcmp(gApplication->Argv(i),"--quench"))) + sprintf(squench,gApplication->Argv(i+1)); + + if (!(strcmp(gApplication->Argv(i),"--qhat"))) + sprintf(sqhat,gApplication->Argv(i+1)); + + } + + if(!(strcmp(sminpthard,""))){ + sprintf(sminpthard,"%d",ptHardLo[npthardbin%nPtHardBins]); + Printf(">>> MinPtHard %s %d",sminpthard,npthardbin%nPtHardBins); + } + if(!(strcmp(smaxpthard,""))){ + sprintf(smaxpthard,"%d",ptHardHi[npthardbin%nPtHardBins]); + Printf(">>> MaxPtHard %s %d",smaxpthard,npthardbin%nPtHardBins); + } + + + + /* + char cmd[200] ; + sprintf(cmd, ".! tar zxvf DBpp.tgz") ; + gROOT->ProcessLine(cmd) ; + sprintf(cmd, ".! tar zxvf TPCCalib_v4-16-Rev-01.tgz") ; + gROOT->ProcessLine(cmd) ; + */ + gSystem->Exec("echo $PWD "); + gSystem->Exec("ls -l "); + + seed = nrun * 100000 + naliencounter * 1000 + nevent; + sprintf(sseed,"%d",seed); + + if (seed==0) { + fprintf(stderr,"!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!\n"); + fprintf(stderr,"!!!! WARNING! Seeding variable for MC is 0 !!!!\n"); + fprintf(stderr,"!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!\n"); + } else { + fprintf(stdout,"!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!\n"); + fprintf(stdout,"!!! MC Seed is %d \n",seed); + fprintf(stdout,"!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!\n"); + } + + // set the seed environment variable + gSystem->Setenv("CONFIG_SEED",sseed); + gSystem->Setenv("DC_RUN",srun); + gSystem->Setenv("DC_EVENT",sevent); + gSystem->Setenv("DC_RUN_TYPE",sprocess);//"kPyGammaJetPHOS"); + gSystem->Setenv("PTHARDMIN",sminpthard);//"20"); + gSystem->Setenv("PTHARDMAX",smaxpthard);//"30"); +// gSystem->Setenv("DC_RUN_TYPE","kPyJetJet"); +// gSystem->Setenv("PTHARDMIN","40"); +// gSystem->Setenv("PTHARDMAX","-1"); + + gSystem->Setenv("PTGAMMAPI0MIN",sminptgammapi0);//"1"); + gSystem->Setenv("QUENCHING",squench); + gSystem->Setenv("QHAT",sqhat); +// gSystem->Setenv("QUENCHING","0"); +// gSystem->Setenv("QHAT","0"); + + gSystem->Setenv("ECMS","10000"); +// gSystem->Setenv("ECMS","14000"); + gSystem->Setenv("ALIMDC_RAWDB1","./mdc1"); + gSystem->Setenv("ALIMDC_RAWDB2","./mdc2"); + gSystem->Setenv("ALIMDC_TAGDB","./mdc1/tag"); + gSystem->Setenv("ALIMDC_RUNDB","./mdc1/meta"); + cout<< "SIMRUN:: Run " << gSystem->Getenv("DC_RUN") << " Event " << gSystem->Getenv("DC_EVENT") + << " Process " << gSystem->Getenv("DC_RUN_TYPE") + << " minpthard " << gSystem->Getenv("PTHARDMIN") + << " maxpthard " << gSystem->Getenv("PTHARDMAX") + << " minpt " << gSystem->Getenv("PTGAMMAPI0MIN") + << endl; + //gSystem->Exec("cp $ROOTSYS/etc/system.rootrc .rootrc"); + cout<<">>>>> SIMULATION <<<<<"<Exec("echo ALICE_ROOT $ALICE_ROOT"); + gSystem->Exec("echo ROOTSYS $ROOTSYS"); + gSystem->Exec("echo LD_LIB $LD_LIBRARY_PATH"); + gSystem->Exec("echo PATH $PATH"); + gSystem->Exec("ls -l $ALICE_ROOT/bin/"); + gSystem->Exec("ls -l $ALICE_ROOT/bin/*"); + + gSystem->Exec(Form("%s -b -q \"sim.C(%d)\" > sim.log 2>&1",program,nevent)); + + + // residual + + + gSystem->MakeDirectory("residual"); + gSystem->Exec("rm residual/*.root"); + gSystem->Exec("ln -s ${PWD}/*.root residual/"); + + // residual_trd + gSystem->MakeDirectory("residual_trd"); + gSystem->Exec("rm residual_trd/*.root"); + gSystem->Exec("ln -s ${PWD}/*.root residual_trd/"); + + gSystem->ChangeDirectory("residual"); + cout<<">>>>> RECONSTRUCTION <<<<<"<Exec(Form("%s -b -q ../rec.C > rec.log 2>&1",program)); + Printf("%d",iStatus); + cout<<">>>>> TAG <<<<<"<Exec("aliroot -b -q ../tag.C 2>&1 tag.log"); + iStatus = gSystem->Exec(Form("%s -b -q ../tag.C > tag.log 2>&1",program)); + Printf("%d",iStatus); + cout<<">>>>> CHECK ESD <<<<<"<Exec(Form("%s -b -q ../CheckESD.C > check.log 2>&1",program)); + cout<<">>>>> CREATE AOD <<<<<"<Exec(Form("%s -b -q ../runAODFilterJets.C > aod.log 2>&1",program)); + Printf("%d",iStatus); + + gSystem->Exec(Form("mv Run*tag.root ../")); + gSystem->Exec(Form("mv *RecPoints.root ../")); // for debugging + TString addString(""); + // logs + gSystem->Exec(Form("mv rec.log ../rec%s.log",addString.Data())); + gSystem->Exec(Form("mv tag.log ../tag%s.log",addString.Data())); + gSystem->Exec(Form("mv check.log ../check%s.log",addString.Data())); + gSystem->Exec(Form("mv aod.log ../aod%s.log",addString.Data())); + // roots + gSystem->Exec(Form("mv AliESDs.root ../AliESDs%s.root",addString.Data())); + gSystem->Exec(Form("mv AliESDfriends.root ../AliESDfriends%s.root",addString.Data())); + gSystem->Exec(Form("mv AliAOD.root ../AliAOD%s.root",addString.Data())); + gSystem->ChangeDirectory("../"); + + + + + // residual_trd + + gSystem->ChangeDirectory("residual_trd"); + /* + sprintf(cmd, ".! tar zxvf ../DBpp.tgz") ; + gROOT->ProcessLine(cmd) ; + sprintf(cmd, ".! tar zxvf ../TPCCalib_v4-16-Rev-01.tgz") ; + gROOT->ProcessLine(cmd) ; + */ + gSystem->Exec("ls -l"); + cout<<">>>>> REC <<<<<"<Exec(Form("%s -b -q ../rec2.C > rec.log 2>&1",program)); + + Printf("%d",iStatus); + cout<<">>>>> CHECK ESD <<<<<"<Exec(Form("%s -b -q ../CheckESD.C > check.log 2>&1",program)); + Printf("%d",iStatus); + cout<<">>>>> CREATE AOD <<<<<"<Exec(Form("%s -b -q ../runAODFilterJets.C > aod.log 2>&1",program)); + Printf("%d",iStatus); + TString addString("TRD"); + // logs + gSystem->Exec(Form("mv rec.log ../rec%s.log",addString.Data())); + gSystem->Exec(Form("mv check.log ../check%s.log",addString.Data())); + gSystem->Exec(Form("mv aod.log ../aod%s.log",addString.Data())); + // roots + gSystem->Exec(Form("mv AliESDs.root ../AliESDs%s.root",addString.Data())); + gSystem->Exec(Form("mv AliESDfriends.root ../AliESDfriends%s.root",addString.Data())); + gSystem->Exec(Form("mv AliAOD.root ../AliAOD%s.root",addString.Data())); + gSystem->ChangeDirectory("../"); + gSystem->Exec("ls -l"); +} diff --git a/prod/LHC09a1/tag.C b/prod/LHC09a1/tag.C new file mode 100644 index 00000000000..3cd2e7f2da5 --- /dev/null +++ b/prod/LHC09a1/tag.C @@ -0,0 +1,109 @@ +void tag() { + const char* turl = gSystem->Getenv("ALIEN_JDL_OUTPUTDIR"); + TString fESDFileName = "alien://"; + fESDFileName += turl; + fESDFileName += "/AliESDs.root"; + + TString fGUID = 0; + GetGUID(fGUID); + + TString fAliroot, fRoot, fGeant; + GetVersions(fAliroot,fRoot,fGeant); + + UpdateTag(fAliroot,fRoot,fGeant,fESDFileName,fGUID); +} + +//_____________________________________// +GetVersions(TString &fAliroot, TString &froot, TString &fgeant) { + const char* fver = gSystem->Getenv("ALIEN_JDL_PACKAGES"); + TString fS = fver; + Int_t fFirst = fS.First("#"); + + while(fFirst != -1) { + Int_t fTotalLength = fS.Length(); + TString tmp = fS; + TString fS1 = fS(0,fFirst); + tmp = fS(fFirst+2,fTotalLength); + fS = tmp; + + if(fS1.Contains("Root")) fAliroot = fS1; + if(fS1.Contains("ROOT")) froot = fS1; + if(fS1.Contains("GEANT")) fgeant = fS1; + + if(tmp.Contains("Root")) fAliroot = tmp; + if(tmp.Contains("ROOT")) froot = tmp; + if(tmp.Contains("GEANT")) fgeant = tmp; + + fFirst = tmp.First("#"); + } +} + +//_____________________________________// +GetGUID(TString &guid) { + ofstream myfile ("guid.txt"); + if (myfile.is_open()) { + TFile *f = TFile::Open("AliESDs.root","read"); + if(f->IsOpen()) { + guid = f->GetUUID().AsString(); + myfile << "AliESDs.root \t"<GetUUID().AsString(); + cout<OpenDirectory(gSystem->pwd()); + const char * name = 0x0; + // Add all files matching *pattern* to the chain + while((name = gSystem->GetDirEntry(dirp))) { + if (strstr(name,tagPattern)) { + TFile *f = TFile::Open(name,"read") ; + + AliRunTag *tag = new AliRunTag; + AliEventTag *evTag = new AliEventTag; + TTree *fTree = (TTree *)f->Get("T"); + fTree->SetBranchAddress("AliTAG",&tag); + + //Defining new tag objects + AliRunTag *newTag = new AliRunTag(); + TTree ttag("T","A Tree with event tags"); + TBranch * btag = ttag.Branch("AliTAG", &newTag); + btag->SetCompressionLevel(9); + for(Int_t iTagFiles = 0; iTagFiles < fTree->GetEntries(); iTagFiles++) { + fTree->GetEntry(iTagFiles); + newTag->SetRunId(tag->GetRunId()); + newTag->SetAlirootVersion(faliroot); + newTag->SetRootVersion(froot); + newTag->SetGeant3Version(fgeant); + const TClonesArray *tagList = tag->GetEventTags(); + for(Int_t j = 0; j < tagList->GetEntries(); j++) { + evTag = (AliEventTag *) tagList->At(j); + evTag->SetTURL(turl); + evTag->SetGUID(guid); + newTag->AddEventTag(*evTag); + } + ttag.Fill(); + newTag->Clear(); + }//tag file loop + + TFile* ftag = TFile::Open(name, "recreate"); + ftag->cd(); + ttag.Write(); + ftag->Close(); + + delete tag; + delete newTag; + }//pattern check + }//directory loop + return kTRUE; +} + diff --git a/prod/LHC09a1/validation.sh b/prod/LHC09a1/validation.sh new file mode 100644 index 00000000000..20b7f14f428 --- /dev/null +++ b/prod/LHC09a1/validation.sh @@ -0,0 +1,77 @@ +#!/bin/sh +################################################## +validateout=`dirname $0` +validatetime=`date` +validated="0"; +error=1 +if [ -z $validateout ] +then + validateout="." +fi + +cd $validateout; +validateworkdir=`pwd`; + +echo "*******************************************************" >> stdout; +echo "* AliRoot Validation Script V1.0 *" >> stdout; +echo "* Time: $validatetime " >> stdout; +echo "* Dir: $validateout" >> stdout; +echo "* Workdir: $validateworkdir" >> stdout; +echo "* ----------------------------------------------------*" >> stdout; +ls -la ./ >> stdout; +echo "* ----------------------------------------------------*" >> stdout; + +################################################## +if [ -f rec.log ] && [ -f sim.log ] && [ -f check.log ] && [ -f tag.log ] +then +sv=`grep -i "Segmentation violation" *.log` +if [ "$sv" = "" ] + then + sf=`grep -i "Segmentation fault" *.log` + if [ "$sf" = "" ] + then + be=`grep -i "Bus error" *.log` + if [ "$be" = "" ] + then + ab=`grep -i "Abort" *.log` + if [ "$ab" = "" ] + then + fp=`grep -i "Floating point exception" *.log` + if [ "$fp" = "" ] + then + kl=`grep -i "Killed" *.log` + if [ "$kl" = "" ] + then + bf=`grep -i "busy flag cleared" *.log` + if [ "$bf" = "" ] + then + ch=`grep -i "check of ESD was successfull" check.log` + if [ "$ch" = "" ] + then + echo "* # The ESD was not successfully checked *" >>stdout; + else + echo "* ---------------- Job Validated ------------------*" >> stdout; + error="0"; + fi + else + echo "* # Check Macro failed ! #" >> stdout; + fi; + fi; + fi + fi + fi + fi +fi +else + echo "* ########## Job not validated - no rec.log or sim.log or check.log ###" >> stdout; +fi +if [ "$error" = "1" ] + then + echo "* ################ Job not validated ################" >> stdout; +fi +echo "* ----------------------------------------------------*" >> stdout; +echo "*******************************************************" >> stdout; +sleep 15; +cd - +exit $error +