From: morsch Date: Wed, 25 Nov 2009 09:23:25 +0000 (+0000) Subject: Recent first physics productions. X-Git-Url: http://git.uio.no/git/?a=commitdiff_plain;h=b0e91ec7c0a8ad3bf09b36f66c03aea1841f039d;p=u%2Fmrichter%2FAliRoot.git Recent first physics productions. --- diff --git a/prod/LHC09d2/CheckESD.C b/prod/LHC09d2/CheckESD.C new file mode 100644 index 00000000000..1f700d32b40 --- /dev/null +++ b/prod/LHC09d2/CheckESD.C @@ -0,0 +1,693 @@ +#if !defined( __CINT__) || defined(__MAKECINT__) +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#include "AliRunLoader.h" +#include "AliLoader.h" +#include "AliESDEvent.h" +#include "AliESDv0.h" +#include "AliESDcascade.h" +#include "AliESDMuonTrack.h" +#include "AliESDCaloCluster.h" +#include "AliRun.h" +#include "AliStack.h" +#include "AliHeader.h" +#include "AliGenEventHeader.h" +#include "AliPID.h" +#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[AliPID::kSPECIES] = + {kElectron, kMuonMinus, kPiPlus, kKPlus, kProton}; + const char* partName[AliPID::kSPECIES+1] = + {"electron", "muon", "pion", "kaon", "proton", "other"}; + Double_t partFrac[AliPID::kSPECIES] = + {0.01, 0.01, 0.85, 0.10, 0.05}; + Int_t identified[AliPID::kSPECIES+1][AliPID::kSPECIES]; + for (Int_t iGen = 0; iGen < AliPID::kSPECIES+1; iGen++) { + for (Int_t iRec = 0; iRec < AliPID::kSPECIES; 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, 50, "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 = runLoader->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 + hResPtInv->Fill(100. * (TMath::Abs(track->GetSigned1Pt()) - 1./particle->Pt()) * + particle->Pt()); + hResPhi->Fill(1000. * (track->Phi() - particle->Phi())); + hResTheta->Fill(1000. * (track->Theta() - particle->Theta())); + + // PID + if ((track->GetStatus() & AliESDtrack::kESDpid) == 0) continue; + Int_t iGen = 5; + for (Int_t i = 0; i < AliPID::kSPECIES; i++) { + if (TMath::Abs(particle->GetPdgCode()) == partCode[i]) iGen = i; + } + Double_t probability[AliPID::kSPECIES]; + track->GetESDpid(probability); + Double_t pMax = 0; + Int_t iRec = 0; + for (Int_t i = 0; i < AliPID::kSPECIES; 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[AliPID::kSPECIES]; + track->GetIntegratedTimes(time); + if (iGen == iRec) { + hDEdxRight->Fill(particle->P(), track->GetTPCsignal()); + if ((track->GetStatus() & AliESDtrack::kTOFpid) != 0) { + hResTOFRight->Fill(track->GetTOFsignal() - time[iRec]); + } + } else { + hDEdxWrong->Fill(particle->P(), 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->GetEffMassXi()); + cascade->ChangeMassHypothesis(v0q,kOmegaMinus); + hMassOmega->Fill(cascade->GetEffMassXi()); + + 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 clusters + { + for (Int_t iCluster=0; iClusterGetNumberOfCaloClusters(); iCluster++) { + AliESDCaloCluster * clust = esd->GetCaloCluster(iCluster); + if (clust->IsPHOS()) hEPHOS->Fill(clust->E()); + if (clust->IsEMCAL()) hEEMCAL->Fill(clust->E()); + } + } + + } + + // perform checks + if (nGen < checkNGenLow) { + Warning("CheckESD", "low number of generated particles: %d", Int_t(nGen)); + } + + TH1F* hEff = CreateEffHisto(hGen, hRec); + + Info("CheckESD", "%d out of %d tracks reconstructed including %d " + "fake tracks", nRec, nGen, nFake); + if (nGen > 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/LHC09d2/Config.C b/prod/LHC09d2/Config.C new file mode 100644 index 00000000000..20921a46ac1 --- /dev/null +++ b/prod/LHC09d2/Config.C @@ -0,0 +1,557 @@ +// +// Configuration for the first physics production 2008 +// + +// 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.C++") + +#if !defined(__CINT__) || defined(__MAKECINT__) +#include +#include +#include +#include +#include +#include +#include "STEER/AliRunLoader.h" +#include "STEER/AliRun.h" +#include "STEER/AliConfig.h" +#include "PYTHIA6/AliDecayerPythia.h" +#include "PYTHIA6/AliGenPythia.h" +#include "TDPMjet/AliGenDPMjet.h" +#include "STEER/AliMagFCheb.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/AliITSv11Hybrid.h" +#include "TPC/AliTPCv2.h" +#include "TOF/AliTOFv6T0.h" +#include "HMPID/AliHMPIDv3.h" +#include "ZDC/AliZDCv3.h" +#include "TRD/AliTRDv1.h" +#include "TRD/AliTRDgeometry.h" +#include "FMD/AliFMDv1.h" +#include "MUON/AliMUONv1.h" +#include "PHOS/AliPHOSv1.h" +#include "PHOS/AliPHOSSimParam.h" +#include "PMD/AliPMDv1.h" +#include "T0/AliT0v1.h" +#include "EMCAL/AliEMCALv2.h" +#include "ACORDE/AliACORDEv1.h" +#include "VZERO/AliVZEROv7.h" +#endif + + +enum PDC06Proc_t +{ + kPythia6, kPythia6D6T, kPhojet, kRunMax +}; + +const char * pprRunName[] = { + "kPythia6", "kPythia6D6T", "kPhojet" +}; + +enum Mag_t +{ + kNoField, k5kG, kFieldMax +}; + +const char * pprField[] = { + "kNoField", "k5kG" +}; + +//--- Functions --- +class AliGenPythia; +AliGenerator *MbPythia(); +AliGenerator *MbPythiaTuneD6T(); +AliGenerator *MbPhojet(); +void ProcessEnvironmentVars(); + +// Geterator, field, beam energy +static PDC06Proc_t proc = kPhojet; +static Mag_t mag = k5kG; +static Float_t energy = 10000; // energy in CMS +//========================// +// Set Random Number seed // +//========================// +TDatime dt; +static UInt_t seed = dt.Get(); + +// Comment line +static TString comment; + +void Config() +{ + + + // Get settings from environment variables + ProcessEnvironmentVars(); + + gRandom->SetSeed(seed); + cerr<<"Seed for random number generation= "<Load("liblhapdf"); // Parton density functions + gSystem->Load("libEGPythia6"); // TGenerator interface + if (proc != kPythia6D6T) { + gSystem->Load("libpythia6"); // Pythia 6.2 + } else { + gSystem->Load("libqpythia"); // Pythia 6.4 + } + gSystem->Load("libAliPythia6"); // ALICE specific implementations + gSystem->Load("libgeant321"); +#endif + + new TGeant3TGeo("C++ Interface to Geant3"); + + //======================================================================= + // Create the output file + + + AliRunLoader* rl=0x0; + + cout<<"Config.C: Creating Run Loader ..."<Fatal("Config.C","Can not instatiate the Run Loader"); + return; + } + rl->SetCompressionLevel(2); + rl->SetNumberOfEventsPerFile(1000); + gAlice->SetRunLoader(rl); + // gAlice->SetGeometryFromFile("geometry.root"); + // gAlice->SetGeometryFromCDB(); + + // Set the trigger configuration: proton-proton + gAlice->SetTriggerDescriptor("p-p"); + + // + //======================================================================= + // ************* STEERING parameters FOR ALICE SIMULATION ************** + // --- Specify event type to be tracked through the ALICE setup + // --- All positions are in cm, angles in degrees, and P and E in GeV + + + gMC->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(); + decayer->SetForceDecay(kAll); + decayer->Init(); + gMC->SetExternalDecayer(decayer); + + //=========================// + // Generator Configuration // + //=========================// + AliGenerator* gener = 0x0; + + if (proc == kPythia6) { + gener = MbPythia(); + } else if (proc == kPythia6D6T) { + gener = MbPythiaTuneD6T(); + } else if (proc == kPhojet) { + gener = MbPhojet(); + } + + + + // PRIMARY VERTEX + // + gener->SetOrigin(0., 0., 0.); // vertex position + // + // + // Size of the interaction diamond + // Longitudinal + Float_t sigmaz = 5.4 / TMath::Sqrt(2.); // [cm] + if (energy == 900) + sigmaz = 10.5 / TMath::Sqrt(2.); // [cm] + + if (energy == 7000) + sigmaz = 6.3 / TMath::Sqrt(2.); // [cm] + + // + // Transverse + Float_t betast = 10; // beta* [m] + Float_t eps = 3.75e-6; // emittance [m] + Float_t gamma = energy / 2.0 / 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(); + if (proc == kPythia6D6T) { + // F I X + AliPythia::Instance()->Pytune(109); + AliPythia::Instance()->Initialize("CMS","p","p", energy); + // F I X + } + // FIELD + // + AliMagF* field = 0x0; + if (mag == kNoField) { + comment = comment.Append(" | L3 field 0.0 T"); + field = new AliMagF("Maps","Maps", 0., 0., AliMagF::k5kGUniform,AliMagF::kBeamTypepp, energy/2.0); + } else if (mag == k5kG) { + comment = comment.Append(" | L3 field 0.5 T"); + field = new AliMagF("Maps","Maps", -1., -1., AliMagF::k5kG, AliMagF::kBeamTypepp, energy/2.0); + } + + printf("\n \n Comment: %s \n \n", comment.Data()); + + TGeoGlobalMagField::Instance()->SetField(field); + + rl->CdGAFile(); + + 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 = 1; + 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"); + FRAME->SetHoles(1); + } + + 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"); + } + + 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 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,1,7,8,9,16,17 + // 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); + geoTRD->SetSMstatus(16,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"); + } + + if (iPHOS) + { + //=================== PHOS parameters =========================== + + AliPHOS *PHOS = new AliPHOSv1("PHOS", "IHEP"); + //Set simulation parameters different from the default ones. + AliPHOSSimParam* simEmc = AliPHOSSimParam::GetInstance() ; + + // APD noise of warm (+20C) PHOS: + // a2 = a1*(Y1/Y2)*(M1/M2), where a1 = 0.012 is APD noise at -25C, + // Y1 = 4.3 photo-electrons/MeV, Y2 = 1.7 p.e/MeV - light yields at -25C and +20C, + // M1 = 50, M2 = 50 - APD gain factors chosen for t1 = -25C and t2 = +20C, + // Y = MeanLightYield*APDEfficiency. + + Float_t apdNoise = 0.012*2.5; + simEmc->SetAPDNoise(apdNoise); + + //Raw Light Yield at +20C + simEmc->SetMeanLightYield(18800); + + //ADC channel width at +18C. + simEmc->SetADCchannelW(0.0125); + } + + + 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) + { + //=================== ACORDE parameters ============================ + + AliACORDE *ACORDE = new AliACORDEv1("ACORDE", "normal ACORDE"); + } + + if (iVZERO) + { + //=================== ACORDE parameters ============================ + + AliVZERO *VZERO = new AliVZEROv7("VZERO", "normal VZERO"); + } +} +// +// PYTHIA +// + +AliGenerator* MbPythia() +{ + comment = comment.Append(" pp at 14 TeV: Pythia low-pt"); +// +// 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(energy); + + return pythia; +} + +AliGenerator* MbPythiaTuneD6T() +{ + comment = comment.Append(" pp at 14 TeV: Pythia low-pt"); +// +// 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(energy); +// Tune +// 109 D6T : Rick Field's CDF Tune D6T (NB: needs CTEQ6L pdfs externally) +// pythia->SetTune(109); // F I X + pythia->SetStrucFunc(kCTEQ6l); +// + return pythia; +} + +AliGenerator* MbPhojet() +{ + comment = comment.Append(" pp at 14 TeV: Pythia low-pt"); +// +// DPMJET +#if defined(__CINT__) + gSystem->Load("libdpmjet"); // Parton density functions + gSystem->Load("libTDPMjet"); // Parton density functions +#endif + AliGenDPMjet* dpmjet = new AliGenDPMjet(-1); + dpmjet->SetMomentumRange(0, 999999.); + dpmjet->SetThetaRange(0., 180.); + dpmjet->SetYRange(-12.,12.); + dpmjet->SetPtRange(0,1000.); + dpmjet->SetProcess(kDpmMb); + dpmjet->SetEnergyCMS(energy); + + return dpmjet; +} + +void ProcessEnvironmentVars() +{ + // Run type + if (gSystem->Getenv("CONFIG_RUN_TYPE")) { + for (Int_t iRun = 0; iRun < kRunMax; iRun++) { + if (strcmp(gSystem->Getenv("CONFIG_RUN_TYPE"), pprRunName[iRun])==0) { + proc = (PDC06Proc_t)iRun; + cout<<"Run type set to "<Getenv("CONFIG_FIELD")) { + for (Int_t iField = 0; iField < kFieldMax; iField++) { + if (strcmp(gSystem->Getenv("CONFIG_FIELD"), pprField[iField])==0) { + mag = (Mag_t)iField; + cout<<"Field set to "<Getenv("CONFIG_ENERGY")) { + energy = atoi(gSystem->Getenv("CONFIG_ENERGY")); + cout<<"Energy set to "<Getenv("CONFIG_SEED")) { + seed = atoi(gSystem->Getenv("CONFIG_SEED")); + } +} diff --git a/prod/LHC09d2/CreateAODfromESD.C b/prod/LHC09d2/CreateAODfromESD.C new file mode 100644 index 00000000000..9c4148736af --- /dev/null +++ b/prod/LHC09d2/CreateAODfromESD.C @@ -0,0 +1,137 @@ +#if !defined(__CINT__) || defined(__MAKECINT__) +#include +#include +#include "AliAnalysisManager.h" +#include "AliESDInputHandler.h" +#include "AliAODHandler.h" +#include "AliAnalysisTaskESDfilter.h" +#include "AliAnalysisDataContainer.h" +#endif + +void CreateAODfromESD(const char *inFileName = "AliESDs.root", + const char *outFileName = "AliAODs.root", + Bool_t bKineFilter = kTRUE) +{ + + gSystem->Load("libTree"); + gSystem->Load("libGeom"); + gSystem->Load("libPhysics"); + gSystem->Load("libVMC"); + gSystem->Load("libSTEERBase"); + gSystem->Load("libESD"); + gSystem->Load("libAOD"); + + gSystem->Load("libANALYSIS"); + gSystem->Load("libANALYSISalice"); + gSystem->Load("libCORRFW"); + gSystem->Load("libPWG3muon"); + + TChain *chain = new TChain("esdTree"); + // Steering input chain + chain->Add(inFileName); + AliAnalysisManager *mgr = new AliAnalysisManager("ESD to AOD", "Analysis Manager"); + + // Input + AliESDInputHandler* inpHandler = new AliESDInputHandler(); + inpHandler->SetReadTags(); + mgr->SetInputEventHandler (inpHandler); + // Output + AliAODHandler* aodHandler = new AliAODHandler(); + aodHandler->SetOutputFileName(outFileName); + mgr->SetOutputEventHandler(aodHandler); + + // MC Truth + if(bKineFilter){ + AliMCEventHandler* mcHandler = new AliMCEventHandler(); + mgr->SetMCtruthEventHandler(mcHandler); + } + + + // Tasks + // 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); + + // Barrel Tracks + AliAnalysisTaskESDfilter *filter = new AliAnalysisTaskESDfilter("Filter"); + mgr->AddTask(filter); + // Muons + AliAnalysisTaskESDMuonFilter *esdmuonfilter = new AliAnalysisTaskESDMuonFilter("ESD Muon Filter"); + mgr->AddTask(esdmuonfilter); + + // Cuts on primary tracks + AliESDtrackCuts* esdTrackCutsL = new AliESDtrackCuts("AliESDtrackCuts", "Standard"); + esdTrackCutsL->SetMinNClustersTPC(50); + esdTrackCutsL->SetMaxChi2PerClusterTPC(3.5); + esdTrackCutsL->SetMaxCovDiagonalElements(2, 2, 0.5, 0.5, 2); + esdTrackCutsL->SetRequireTPCRefit(kTRUE); + esdTrackCutsL->SetMaxDCAToVertexXY(3.0); + esdTrackCutsL->SetMaxDCAToVertexZ(3.0); + esdTrackCutsL->SetDCAToVertex2D(kTRUE); + esdTrackCutsL->SetRequireSigmaToVertex(kFALSE); + esdTrackCutsL->SetAcceptKinkDaughters(kFALSE); + // ITS stand-alone tracks + AliESDtrackCuts* esdTrackCutsITSsa = new AliESDtrackCuts("AliESDtrackCuts", "ITS stand-alone"); + esdTrackCutsITSsa->SetRequireITSStandAlone(kTRUE); + + AliAnalysisFilter* trackFilter = new AliAnalysisFilter("trackFilter"); + trackFilter->AddCuts(esdTrackCutsL); + trackFilter->AddCuts(esdTrackCutsITSsa); + + // Cuts on V0s + AliESDv0Cuts* esdV0Cuts = new AliESDv0Cuts("AliESDv0Cuts", "Standard pp"); + esdV0Cuts->SetMinRadius(0.2); + esdV0Cuts->SetMaxRadius(200); + esdV0Cuts->SetMinDcaPosToVertex(0.05); + esdV0Cuts->SetMinDcaNegToVertex(0.05); + esdV0Cuts->SetMaxDcaV0Daughters(1.0); + esdV0Cuts->SetMinCosinePointingAngle(0.99); + AliAnalysisFilter* v0Filter = new AliAnalysisFilter("v0Filter"); + v0Filter->AddCuts(esdV0Cuts); + + +// + filter->SetTrackFilter(trackFilter); + filter->SetV0Filter(v0Filter); + + +// Create AOD Tags + AliAnalysisTaskTagCreator* tagTask = new AliAnalysisTaskTagCreator("AOD Tag Creator"); + mgr->AddTask(tagTask); + + // Pipelining + AliAnalysisDataContainer *cinput1 = mgr->GetCommonInputContainer(); + AliAnalysisDataContainer *coutput1 = mgr->GetCommonOutputContainer(); + + + AliAnalysisDataContainer *coutputT + = mgr->CreateContainer("cTag", TTree::Class(), AliAnalysisManager::kOutputContainer, "AOD.tag.root"); + + coutput1->SetSpecialOutput(); + coutputT->SetSpecialOutput(); + + if(bKineFilter) { + mgr->ConnectInput (kinefilter, 0, cinput1 ); + mgr->ConnectOutput (kinefilter, 0, coutput1 ); + } + + mgr->ConnectInput (filter, 0, cinput1 ); + mgr->ConnectOutput(filter, 0, coutput1); + + mgr->ConnectInput (esdmuonfilter, 0, cinput1 ); +// mgr->ConnectOutput(esdmuonfilter, 0, coutput1); + + mgr->ConnectInput (tagTask, 0, cinput1); + mgr->ConnectOutput(tagTask, 1, coutputT); + + // + // Run the analysis + // + mgr->InitAnalysis(); + mgr->PrintStatus(); + mgr->StartAnalysis("local", chain); +} + diff --git a/prod/LHC09d2/JDL b/prod/LHC09d2/JDL new file mode 100644 index 00000000000..9e11552b567 --- /dev/null +++ b/prod/LHC09d2/JDL @@ -0,0 +1,31 @@ +Executable = "aliroot_new"; +Jobtag={"comment:Early physics (2009 - stage 3) pp, Pythia6 Tune D6T, 0T, 900GeV, ID #129"}; + +Packages={"VO_ALICE@AliRoot::v4-17-Rev-15","VO_ALICE@GEANT3::v1-11-4","VO_ALICE@ROOT::v20091104-1","VO_ALICE@APISCONFIG::V1.0x"}; + +TTL = "72000"; + +Requirements = member(other.GridPartitions,"fastmc"); + +Validationcommand ="/alice/cern.ch/user/a/aliprod/prod2007/configs_pbpb_hijing/validation.sh"; + +InputFile= {"LF:/alice/cern.ch/user/a/aliprod/LHC09d2/CheckESD.C", + "LF:/alice/cern.ch/user/a/aliprod/LHC09d2/Config.C", + "LF:/alice/cern.ch/user/a/aliprod/LHC09d2/rec.C", + "LF:/alice/cern.ch/user/a/aliprod/LHC09d2/sim.C", + "LF:/alice/cern.ch/user/a/aliprod/LHC09d2/simrun.C", + "LF:/alice/cern.ch/user/a/aliprod/LHC09d2/tag.C", + "LF:/alice/cern.ch/user/a/aliprod/LHC09d2/CreateAODfromESD.C"}; + + +OutputArchive={"log_archive:*.log,stdout,stderr@ALICE::JINR::SE","root_archive.zip:galice.root,Kinematics.root,TrackRefs.root,AliESDs.root,AliESDfriends.root,AliAODs.root,*QA*.root,ITS.RecPoints.root,AliITSPlaneEffSPDtracklet.root,AliITSPlaneEffSPDtrackletBkg.root,TrackletsMCpred.root,TrackleterSPDHistos.root,Run*.root@ALICE::CERN::ALICEDISK,ALICE::FZK::SE,ALICE::CNAF::SE"}; + +OutputDir="/alice/sim/LHC09d2/$1/#alien_counter_03i#"; + +JDLVariables={"Packages", "OutputDir"}; +GUIDFILE="guid.txt"; + +splitarguments="simrun.C --run $1 --event #alien_counter# --process kPythia6D6T --field kNoField --energy 900" ; +split="production:1-1000"; + +Workdirectorysize={"10000MB"}; diff --git a/prod/LHC09d2/rec.C b/prod/LHC09d2/rec.C new file mode 100644 index 00000000000..f1810f3ac19 --- /dev/null +++ b/prod/LHC09d2/rec.C @@ -0,0 +1,24 @@ +void rec() { + + AliReconstruction reco; + + reco.SetWriteESDfriend(); + reco.SetWriteAlignmentData(); + reco.SetRunPlaneEff(kTRUE); + AliITSRecoParam *itspar = AliITSRecoParam::GetPlaneEffParam(-1); + itspar->SetOptTrackletsPlaneEff(kTRUE); + reco.SetRecoParam("ITS",itspar); + + reco.SetDefaultStorage("alien://Folder=/alice/simulation/2008/v4-15-Release/Residual/"); +// reco.SetSpecificStorage("GRP/GRP/Data", +// Form("local://%s",gSystem->pwd())); + // We store the object in AliEn during the simulation + reco.SetSpecificStorage("GRP/GRP/Data", + "alien://Folder=/alice/simulation/2008/v4-15-Release/Ideal/"); + + TStopwatch timer; + timer.Start(); + reco.Run(); + timer.Stop(); + timer.Print(); +} diff --git a/prod/LHC09d2/sim.C b/prod/LHC09d2/sim.C new file mode 100644 index 00000000000..46720580c86 --- /dev/null +++ b/prod/LHC09d2/sim.C @@ -0,0 +1,20 @@ +void sim(Int_t nev=200) { + + AliSimulation simulator; + simulator.SetMakeSDigits("TRD TOF PHOS HMPID EMCAL MUON FMD ZDC PMD T0 VZERO"); + simulator.SetMakeDigitsFromHits("ITS TPC"); + + // The raw data are not written due to the huge increase of the + // virtual memory in HLT + // simulator.SetWriteRawData("ALL","raw.root",kTRUE); + + simulator.SetDefaultStorage("alien://Folder=/alice/simulation/2008/v4-15-Release/Ideal/"); +// simulator.SetSpecificStorage("GRP/GRP/Data", +// Form("local://%s",gSystem->pwd())); + + TStopwatch timer; + timer.Start(); + simulator.Run(nev); + timer.Stop(); + timer.Print(); +} diff --git a/prod/LHC09d2/simrun.C b/prod/LHC09d2/simrun.C new file mode 100644 index 00000000000..068cb6b07ee Binary files /dev/null and b/prod/LHC09d2/simrun.C differ diff --git a/prod/LHC09d2/tag.C b/prod/LHC09d2/tag.C new file mode 100644 index 00000000000..668d95ae835 --- /dev/null +++ b/prod/LHC09d2/tag.C @@ -0,0 +1,108 @@ +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/LHC09d3/CheckESD.C b/prod/LHC09d3/CheckESD.C new file mode 100644 index 00000000000..1f700d32b40 --- /dev/null +++ b/prod/LHC09d3/CheckESD.C @@ -0,0 +1,693 @@ +#if !defined( __CINT__) || defined(__MAKECINT__) +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#include "AliRunLoader.h" +#include "AliLoader.h" +#include "AliESDEvent.h" +#include "AliESDv0.h" +#include "AliESDcascade.h" +#include "AliESDMuonTrack.h" +#include "AliESDCaloCluster.h" +#include "AliRun.h" +#include "AliStack.h" +#include "AliHeader.h" +#include "AliGenEventHeader.h" +#include "AliPID.h" +#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[AliPID::kSPECIES] = + {kElectron, kMuonMinus, kPiPlus, kKPlus, kProton}; + const char* partName[AliPID::kSPECIES+1] = + {"electron", "muon", "pion", "kaon", "proton", "other"}; + Double_t partFrac[AliPID::kSPECIES] = + {0.01, 0.01, 0.85, 0.10, 0.05}; + Int_t identified[AliPID::kSPECIES+1][AliPID::kSPECIES]; + for (Int_t iGen = 0; iGen < AliPID::kSPECIES+1; iGen++) { + for (Int_t iRec = 0; iRec < AliPID::kSPECIES; 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, 50, "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 = runLoader->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 + hResPtInv->Fill(100. * (TMath::Abs(track->GetSigned1Pt()) - 1./particle->Pt()) * + particle->Pt()); + hResPhi->Fill(1000. * (track->Phi() - particle->Phi())); + hResTheta->Fill(1000. * (track->Theta() - particle->Theta())); + + // PID + if ((track->GetStatus() & AliESDtrack::kESDpid) == 0) continue; + Int_t iGen = 5; + for (Int_t i = 0; i < AliPID::kSPECIES; i++) { + if (TMath::Abs(particle->GetPdgCode()) == partCode[i]) iGen = i; + } + Double_t probability[AliPID::kSPECIES]; + track->GetESDpid(probability); + Double_t pMax = 0; + Int_t iRec = 0; + for (Int_t i = 0; i < AliPID::kSPECIES; 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[AliPID::kSPECIES]; + track->GetIntegratedTimes(time); + if (iGen == iRec) { + hDEdxRight->Fill(particle->P(), track->GetTPCsignal()); + if ((track->GetStatus() & AliESDtrack::kTOFpid) != 0) { + hResTOFRight->Fill(track->GetTOFsignal() - time[iRec]); + } + } else { + hDEdxWrong->Fill(particle->P(), 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->GetEffMassXi()); + cascade->ChangeMassHypothesis(v0q,kOmegaMinus); + hMassOmega->Fill(cascade->GetEffMassXi()); + + 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 clusters + { + for (Int_t iCluster=0; iClusterGetNumberOfCaloClusters(); iCluster++) { + AliESDCaloCluster * clust = esd->GetCaloCluster(iCluster); + if (clust->IsPHOS()) hEPHOS->Fill(clust->E()); + if (clust->IsEMCAL()) hEEMCAL->Fill(clust->E()); + } + } + + } + + // perform checks + if (nGen < checkNGenLow) { + Warning("CheckESD", "low number of generated particles: %d", Int_t(nGen)); + } + + TH1F* hEff = CreateEffHisto(hGen, hRec); + + Info("CheckESD", "%d out of %d tracks reconstructed including %d " + "fake tracks", nRec, nGen, nFake); + if (nGen > 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/LHC09d3/Config.C b/prod/LHC09d3/Config.C new file mode 100644 index 00000000000..20921a46ac1 --- /dev/null +++ b/prod/LHC09d3/Config.C @@ -0,0 +1,557 @@ +// +// Configuration for the first physics production 2008 +// + +// 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.C++") + +#if !defined(__CINT__) || defined(__MAKECINT__) +#include +#include +#include +#include +#include +#include +#include "STEER/AliRunLoader.h" +#include "STEER/AliRun.h" +#include "STEER/AliConfig.h" +#include "PYTHIA6/AliDecayerPythia.h" +#include "PYTHIA6/AliGenPythia.h" +#include "TDPMjet/AliGenDPMjet.h" +#include "STEER/AliMagFCheb.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/AliITSv11Hybrid.h" +#include "TPC/AliTPCv2.h" +#include "TOF/AliTOFv6T0.h" +#include "HMPID/AliHMPIDv3.h" +#include "ZDC/AliZDCv3.h" +#include "TRD/AliTRDv1.h" +#include "TRD/AliTRDgeometry.h" +#include "FMD/AliFMDv1.h" +#include "MUON/AliMUONv1.h" +#include "PHOS/AliPHOSv1.h" +#include "PHOS/AliPHOSSimParam.h" +#include "PMD/AliPMDv1.h" +#include "T0/AliT0v1.h" +#include "EMCAL/AliEMCALv2.h" +#include "ACORDE/AliACORDEv1.h" +#include "VZERO/AliVZEROv7.h" +#endif + + +enum PDC06Proc_t +{ + kPythia6, kPythia6D6T, kPhojet, kRunMax +}; + +const char * pprRunName[] = { + "kPythia6", "kPythia6D6T", "kPhojet" +}; + +enum Mag_t +{ + kNoField, k5kG, kFieldMax +}; + +const char * pprField[] = { + "kNoField", "k5kG" +}; + +//--- Functions --- +class AliGenPythia; +AliGenerator *MbPythia(); +AliGenerator *MbPythiaTuneD6T(); +AliGenerator *MbPhojet(); +void ProcessEnvironmentVars(); + +// Geterator, field, beam energy +static PDC06Proc_t proc = kPhojet; +static Mag_t mag = k5kG; +static Float_t energy = 10000; // energy in CMS +//========================// +// Set Random Number seed // +//========================// +TDatime dt; +static UInt_t seed = dt.Get(); + +// Comment line +static TString comment; + +void Config() +{ + + + // Get settings from environment variables + ProcessEnvironmentVars(); + + gRandom->SetSeed(seed); + cerr<<"Seed for random number generation= "<Load("liblhapdf"); // Parton density functions + gSystem->Load("libEGPythia6"); // TGenerator interface + if (proc != kPythia6D6T) { + gSystem->Load("libpythia6"); // Pythia 6.2 + } else { + gSystem->Load("libqpythia"); // Pythia 6.4 + } + gSystem->Load("libAliPythia6"); // ALICE specific implementations + gSystem->Load("libgeant321"); +#endif + + new TGeant3TGeo("C++ Interface to Geant3"); + + //======================================================================= + // Create the output file + + + AliRunLoader* rl=0x0; + + cout<<"Config.C: Creating Run Loader ..."<Fatal("Config.C","Can not instatiate the Run Loader"); + return; + } + rl->SetCompressionLevel(2); + rl->SetNumberOfEventsPerFile(1000); + gAlice->SetRunLoader(rl); + // gAlice->SetGeometryFromFile("geometry.root"); + // gAlice->SetGeometryFromCDB(); + + // Set the trigger configuration: proton-proton + gAlice->SetTriggerDescriptor("p-p"); + + // + //======================================================================= + // ************* STEERING parameters FOR ALICE SIMULATION ************** + // --- Specify event type to be tracked through the ALICE setup + // --- All positions are in cm, angles in degrees, and P and E in GeV + + + gMC->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(); + decayer->SetForceDecay(kAll); + decayer->Init(); + gMC->SetExternalDecayer(decayer); + + //=========================// + // Generator Configuration // + //=========================// + AliGenerator* gener = 0x0; + + if (proc == kPythia6) { + gener = MbPythia(); + } else if (proc == kPythia6D6T) { + gener = MbPythiaTuneD6T(); + } else if (proc == kPhojet) { + gener = MbPhojet(); + } + + + + // PRIMARY VERTEX + // + gener->SetOrigin(0., 0., 0.); // vertex position + // + // + // Size of the interaction diamond + // Longitudinal + Float_t sigmaz = 5.4 / TMath::Sqrt(2.); // [cm] + if (energy == 900) + sigmaz = 10.5 / TMath::Sqrt(2.); // [cm] + + if (energy == 7000) + sigmaz = 6.3 / TMath::Sqrt(2.); // [cm] + + // + // Transverse + Float_t betast = 10; // beta* [m] + Float_t eps = 3.75e-6; // emittance [m] + Float_t gamma = energy / 2.0 / 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(); + if (proc == kPythia6D6T) { + // F I X + AliPythia::Instance()->Pytune(109); + AliPythia::Instance()->Initialize("CMS","p","p", energy); + // F I X + } + // FIELD + // + AliMagF* field = 0x0; + if (mag == kNoField) { + comment = comment.Append(" | L3 field 0.0 T"); + field = new AliMagF("Maps","Maps", 0., 0., AliMagF::k5kGUniform,AliMagF::kBeamTypepp, energy/2.0); + } else if (mag == k5kG) { + comment = comment.Append(" | L3 field 0.5 T"); + field = new AliMagF("Maps","Maps", -1., -1., AliMagF::k5kG, AliMagF::kBeamTypepp, energy/2.0); + } + + printf("\n \n Comment: %s \n \n", comment.Data()); + + TGeoGlobalMagField::Instance()->SetField(field); + + rl->CdGAFile(); + + 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 = 1; + 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"); + FRAME->SetHoles(1); + } + + 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"); + } + + 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 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,1,7,8,9,16,17 + // 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); + geoTRD->SetSMstatus(16,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"); + } + + if (iPHOS) + { + //=================== PHOS parameters =========================== + + AliPHOS *PHOS = new AliPHOSv1("PHOS", "IHEP"); + //Set simulation parameters different from the default ones. + AliPHOSSimParam* simEmc = AliPHOSSimParam::GetInstance() ; + + // APD noise of warm (+20C) PHOS: + // a2 = a1*(Y1/Y2)*(M1/M2), where a1 = 0.012 is APD noise at -25C, + // Y1 = 4.3 photo-electrons/MeV, Y2 = 1.7 p.e/MeV - light yields at -25C and +20C, + // M1 = 50, M2 = 50 - APD gain factors chosen for t1 = -25C and t2 = +20C, + // Y = MeanLightYield*APDEfficiency. + + Float_t apdNoise = 0.012*2.5; + simEmc->SetAPDNoise(apdNoise); + + //Raw Light Yield at +20C + simEmc->SetMeanLightYield(18800); + + //ADC channel width at +18C. + simEmc->SetADCchannelW(0.0125); + } + + + 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) + { + //=================== ACORDE parameters ============================ + + AliACORDE *ACORDE = new AliACORDEv1("ACORDE", "normal ACORDE"); + } + + if (iVZERO) + { + //=================== ACORDE parameters ============================ + + AliVZERO *VZERO = new AliVZEROv7("VZERO", "normal VZERO"); + } +} +// +// PYTHIA +// + +AliGenerator* MbPythia() +{ + comment = comment.Append(" pp at 14 TeV: Pythia low-pt"); +// +// 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(energy); + + return pythia; +} + +AliGenerator* MbPythiaTuneD6T() +{ + comment = comment.Append(" pp at 14 TeV: Pythia low-pt"); +// +// 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(energy); +// Tune +// 109 D6T : Rick Field's CDF Tune D6T (NB: needs CTEQ6L pdfs externally) +// pythia->SetTune(109); // F I X + pythia->SetStrucFunc(kCTEQ6l); +// + return pythia; +} + +AliGenerator* MbPhojet() +{ + comment = comment.Append(" pp at 14 TeV: Pythia low-pt"); +// +// DPMJET +#if defined(__CINT__) + gSystem->Load("libdpmjet"); // Parton density functions + gSystem->Load("libTDPMjet"); // Parton density functions +#endif + AliGenDPMjet* dpmjet = new AliGenDPMjet(-1); + dpmjet->SetMomentumRange(0, 999999.); + dpmjet->SetThetaRange(0., 180.); + dpmjet->SetYRange(-12.,12.); + dpmjet->SetPtRange(0,1000.); + dpmjet->SetProcess(kDpmMb); + dpmjet->SetEnergyCMS(energy); + + return dpmjet; +} + +void ProcessEnvironmentVars() +{ + // Run type + if (gSystem->Getenv("CONFIG_RUN_TYPE")) { + for (Int_t iRun = 0; iRun < kRunMax; iRun++) { + if (strcmp(gSystem->Getenv("CONFIG_RUN_TYPE"), pprRunName[iRun])==0) { + proc = (PDC06Proc_t)iRun; + cout<<"Run type set to "<Getenv("CONFIG_FIELD")) { + for (Int_t iField = 0; iField < kFieldMax; iField++) { + if (strcmp(gSystem->Getenv("CONFIG_FIELD"), pprField[iField])==0) { + mag = (Mag_t)iField; + cout<<"Field set to "<Getenv("CONFIG_ENERGY")) { + energy = atoi(gSystem->Getenv("CONFIG_ENERGY")); + cout<<"Energy set to "<Getenv("CONFIG_SEED")) { + seed = atoi(gSystem->Getenv("CONFIG_SEED")); + } +} diff --git a/prod/LHC09d3/CreateAODfromESD.C b/prod/LHC09d3/CreateAODfromESD.C new file mode 100644 index 00000000000..9c4148736af --- /dev/null +++ b/prod/LHC09d3/CreateAODfromESD.C @@ -0,0 +1,137 @@ +#if !defined(__CINT__) || defined(__MAKECINT__) +#include +#include +#include "AliAnalysisManager.h" +#include "AliESDInputHandler.h" +#include "AliAODHandler.h" +#include "AliAnalysisTaskESDfilter.h" +#include "AliAnalysisDataContainer.h" +#endif + +void CreateAODfromESD(const char *inFileName = "AliESDs.root", + const char *outFileName = "AliAODs.root", + Bool_t bKineFilter = kTRUE) +{ + + gSystem->Load("libTree"); + gSystem->Load("libGeom"); + gSystem->Load("libPhysics"); + gSystem->Load("libVMC"); + gSystem->Load("libSTEERBase"); + gSystem->Load("libESD"); + gSystem->Load("libAOD"); + + gSystem->Load("libANALYSIS"); + gSystem->Load("libANALYSISalice"); + gSystem->Load("libCORRFW"); + gSystem->Load("libPWG3muon"); + + TChain *chain = new TChain("esdTree"); + // Steering input chain + chain->Add(inFileName); + AliAnalysisManager *mgr = new AliAnalysisManager("ESD to AOD", "Analysis Manager"); + + // Input + AliESDInputHandler* inpHandler = new AliESDInputHandler(); + inpHandler->SetReadTags(); + mgr->SetInputEventHandler (inpHandler); + // Output + AliAODHandler* aodHandler = new AliAODHandler(); + aodHandler->SetOutputFileName(outFileName); + mgr->SetOutputEventHandler(aodHandler); + + // MC Truth + if(bKineFilter){ + AliMCEventHandler* mcHandler = new AliMCEventHandler(); + mgr->SetMCtruthEventHandler(mcHandler); + } + + + // Tasks + // 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); + + // Barrel Tracks + AliAnalysisTaskESDfilter *filter = new AliAnalysisTaskESDfilter("Filter"); + mgr->AddTask(filter); + // Muons + AliAnalysisTaskESDMuonFilter *esdmuonfilter = new AliAnalysisTaskESDMuonFilter("ESD Muon Filter"); + mgr->AddTask(esdmuonfilter); + + // Cuts on primary tracks + AliESDtrackCuts* esdTrackCutsL = new AliESDtrackCuts("AliESDtrackCuts", "Standard"); + esdTrackCutsL->SetMinNClustersTPC(50); + esdTrackCutsL->SetMaxChi2PerClusterTPC(3.5); + esdTrackCutsL->SetMaxCovDiagonalElements(2, 2, 0.5, 0.5, 2); + esdTrackCutsL->SetRequireTPCRefit(kTRUE); + esdTrackCutsL->SetMaxDCAToVertexXY(3.0); + esdTrackCutsL->SetMaxDCAToVertexZ(3.0); + esdTrackCutsL->SetDCAToVertex2D(kTRUE); + esdTrackCutsL->SetRequireSigmaToVertex(kFALSE); + esdTrackCutsL->SetAcceptKinkDaughters(kFALSE); + // ITS stand-alone tracks + AliESDtrackCuts* esdTrackCutsITSsa = new AliESDtrackCuts("AliESDtrackCuts", "ITS stand-alone"); + esdTrackCutsITSsa->SetRequireITSStandAlone(kTRUE); + + AliAnalysisFilter* trackFilter = new AliAnalysisFilter("trackFilter"); + trackFilter->AddCuts(esdTrackCutsL); + trackFilter->AddCuts(esdTrackCutsITSsa); + + // Cuts on V0s + AliESDv0Cuts* esdV0Cuts = new AliESDv0Cuts("AliESDv0Cuts", "Standard pp"); + esdV0Cuts->SetMinRadius(0.2); + esdV0Cuts->SetMaxRadius(200); + esdV0Cuts->SetMinDcaPosToVertex(0.05); + esdV0Cuts->SetMinDcaNegToVertex(0.05); + esdV0Cuts->SetMaxDcaV0Daughters(1.0); + esdV0Cuts->SetMinCosinePointingAngle(0.99); + AliAnalysisFilter* v0Filter = new AliAnalysisFilter("v0Filter"); + v0Filter->AddCuts(esdV0Cuts); + + +// + filter->SetTrackFilter(trackFilter); + filter->SetV0Filter(v0Filter); + + +// Create AOD Tags + AliAnalysisTaskTagCreator* tagTask = new AliAnalysisTaskTagCreator("AOD Tag Creator"); + mgr->AddTask(tagTask); + + // Pipelining + AliAnalysisDataContainer *cinput1 = mgr->GetCommonInputContainer(); + AliAnalysisDataContainer *coutput1 = mgr->GetCommonOutputContainer(); + + + AliAnalysisDataContainer *coutputT + = mgr->CreateContainer("cTag", TTree::Class(), AliAnalysisManager::kOutputContainer, "AOD.tag.root"); + + coutput1->SetSpecialOutput(); + coutputT->SetSpecialOutput(); + + if(bKineFilter) { + mgr->ConnectInput (kinefilter, 0, cinput1 ); + mgr->ConnectOutput (kinefilter, 0, coutput1 ); + } + + mgr->ConnectInput (filter, 0, cinput1 ); + mgr->ConnectOutput(filter, 0, coutput1); + + mgr->ConnectInput (esdmuonfilter, 0, cinput1 ); +// mgr->ConnectOutput(esdmuonfilter, 0, coutput1); + + mgr->ConnectInput (tagTask, 0, cinput1); + mgr->ConnectOutput(tagTask, 1, coutputT); + + // + // Run the analysis + // + mgr->InitAnalysis(); + mgr->PrintStatus(); + mgr->StartAnalysis("local", chain); +} + diff --git a/prod/LHC09d3/JDL b/prod/LHC09d3/JDL new file mode 100644 index 00000000000..a2308348270 --- /dev/null +++ b/prod/LHC09d3/JDL @@ -0,0 +1,30 @@ +Executable = "aliroot_new"; +Jobtag={"comment:Early physics (2009 - stage 3) pp, Phojet, 0.5 T, 900 GeV, ID #130"}; + +Packages={"VO_ALICE@AliRoot::v4-17-Rev-15","VO_ALICE@GEANT3::v1-11-4","VO_ALICE@ROOT::v20091104-1","VO_ALICE@APISCONFIG::V1.0x"}; + +TTL = "72000"; + +Requirements = member(other.GridPartitions,"fastmc"); + +Validationcommand ="/alice/cern.ch/user/a/aliprod/prod2007/configs_pbpb_hijing/validation.sh"; + +InputFile= {"LF:/alice/cern.ch/user/a/aliprod/LHC09d3/CheckESD.C", + "LF:/alice/cern.ch/user/a/aliprod/LHC09d3/Config.C", + "LF:/alice/cern.ch/user/a/aliprod/LHC09d3/rec.C", + "LF:/alice/cern.ch/user/a/aliprod/LHC09d3/sim.C", + "LF:/alice/cern.ch/user/a/aliprod/LHC09d3/simrun.C", + "LF:/alice/cern.ch/user/a/aliprod/LHC09d3/tag.C", + "LF:/alice/cern.ch/user/a/aliprod/LHC09d3/CreateAODfromESD.C"}; + +OutputArchive={"log_archive:*.log,stdout,stderr@ALICE::JINR::SE","root_archive.zip:galice.root,Kinematics.root,TrackRefs.root,AliESDs.root,AliESDfriends.root,AliAODs.root,*QA*.root,ITS.RecPoints.root,AliITSPlaneEffSPDtracklet.root,AliITSPlaneEffSPDtrackletBkg.root,TrackletsMCpred.root,TrackleterSPDHistos.root,Run*.root@ALICE::CERN::ALICEDISK,ALICE::FZK::SE,ALICE::CNAF::SE"}; + +OutputDir="/alice/sim/LHC09d3/$1/#alien_counter_03i#"; + +JDLVariables={"Packages", "OutputDir"}; +GUIDFILE="guid.txt"; + +splitarguments="simrun.C --run $1 --event #alien_counter# --process kPhojet --field k5kG --energy 900" ; +split="production:1-1000"; + +Workdirectorysize={"10000MB"}; diff --git a/prod/LHC09d3/rec.C b/prod/LHC09d3/rec.C new file mode 100644 index 00000000000..f1810f3ac19 --- /dev/null +++ b/prod/LHC09d3/rec.C @@ -0,0 +1,24 @@ +void rec() { + + AliReconstruction reco; + + reco.SetWriteESDfriend(); + reco.SetWriteAlignmentData(); + reco.SetRunPlaneEff(kTRUE); + AliITSRecoParam *itspar = AliITSRecoParam::GetPlaneEffParam(-1); + itspar->SetOptTrackletsPlaneEff(kTRUE); + reco.SetRecoParam("ITS",itspar); + + reco.SetDefaultStorage("alien://Folder=/alice/simulation/2008/v4-15-Release/Residual/"); +// reco.SetSpecificStorage("GRP/GRP/Data", +// Form("local://%s",gSystem->pwd())); + // We store the object in AliEn during the simulation + reco.SetSpecificStorage("GRP/GRP/Data", + "alien://Folder=/alice/simulation/2008/v4-15-Release/Ideal/"); + + TStopwatch timer; + timer.Start(); + reco.Run(); + timer.Stop(); + timer.Print(); +} diff --git a/prod/LHC09d3/sim.C b/prod/LHC09d3/sim.C new file mode 100644 index 00000000000..46720580c86 --- /dev/null +++ b/prod/LHC09d3/sim.C @@ -0,0 +1,20 @@ +void sim(Int_t nev=200) { + + AliSimulation simulator; + simulator.SetMakeSDigits("TRD TOF PHOS HMPID EMCAL MUON FMD ZDC PMD T0 VZERO"); + simulator.SetMakeDigitsFromHits("ITS TPC"); + + // The raw data are not written due to the huge increase of the + // virtual memory in HLT + // simulator.SetWriteRawData("ALL","raw.root",kTRUE); + + simulator.SetDefaultStorage("alien://Folder=/alice/simulation/2008/v4-15-Release/Ideal/"); +// simulator.SetSpecificStorage("GRP/GRP/Data", +// Form("local://%s",gSystem->pwd())); + + TStopwatch timer; + timer.Start(); + simulator.Run(nev); + timer.Stop(); + timer.Print(); +} diff --git a/prod/LHC09d3/simrun.C b/prod/LHC09d3/simrun.C new file mode 100644 index 00000000000..068cb6b07ee Binary files /dev/null and b/prod/LHC09d3/simrun.C differ diff --git a/prod/LHC09d3/tag.C b/prod/LHC09d3/tag.C new file mode 100644 index 00000000000..668d95ae835 --- /dev/null +++ b/prod/LHC09d3/tag.C @@ -0,0 +1,108 @@ +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/LHC09d4/CheckESD.C b/prod/LHC09d4/CheckESD.C new file mode 100644 index 00000000000..1f700d32b40 --- /dev/null +++ b/prod/LHC09d4/CheckESD.C @@ -0,0 +1,693 @@ +#if !defined( __CINT__) || defined(__MAKECINT__) +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#include "AliRunLoader.h" +#include "AliLoader.h" +#include "AliESDEvent.h" +#include "AliESDv0.h" +#include "AliESDcascade.h" +#include "AliESDMuonTrack.h" +#include "AliESDCaloCluster.h" +#include "AliRun.h" +#include "AliStack.h" +#include "AliHeader.h" +#include "AliGenEventHeader.h" +#include "AliPID.h" +#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[AliPID::kSPECIES] = + {kElectron, kMuonMinus, kPiPlus, kKPlus, kProton}; + const char* partName[AliPID::kSPECIES+1] = + {"electron", "muon", "pion", "kaon", "proton", "other"}; + Double_t partFrac[AliPID::kSPECIES] = + {0.01, 0.01, 0.85, 0.10, 0.05}; + Int_t identified[AliPID::kSPECIES+1][AliPID::kSPECIES]; + for (Int_t iGen = 0; iGen < AliPID::kSPECIES+1; iGen++) { + for (Int_t iRec = 0; iRec < AliPID::kSPECIES; 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, 50, "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 = runLoader->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 + hResPtInv->Fill(100. * (TMath::Abs(track->GetSigned1Pt()) - 1./particle->Pt()) * + particle->Pt()); + hResPhi->Fill(1000. * (track->Phi() - particle->Phi())); + hResTheta->Fill(1000. * (track->Theta() - particle->Theta())); + + // PID + if ((track->GetStatus() & AliESDtrack::kESDpid) == 0) continue; + Int_t iGen = 5; + for (Int_t i = 0; i < AliPID::kSPECIES; i++) { + if (TMath::Abs(particle->GetPdgCode()) == partCode[i]) iGen = i; + } + Double_t probability[AliPID::kSPECIES]; + track->GetESDpid(probability); + Double_t pMax = 0; + Int_t iRec = 0; + for (Int_t i = 0; i < AliPID::kSPECIES; 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[AliPID::kSPECIES]; + track->GetIntegratedTimes(time); + if (iGen == iRec) { + hDEdxRight->Fill(particle->P(), track->GetTPCsignal()); + if ((track->GetStatus() & AliESDtrack::kTOFpid) != 0) { + hResTOFRight->Fill(track->GetTOFsignal() - time[iRec]); + } + } else { + hDEdxWrong->Fill(particle->P(), 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->GetEffMassXi()); + cascade->ChangeMassHypothesis(v0q,kOmegaMinus); + hMassOmega->Fill(cascade->GetEffMassXi()); + + 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 clusters + { + for (Int_t iCluster=0; iClusterGetNumberOfCaloClusters(); iCluster++) { + AliESDCaloCluster * clust = esd->GetCaloCluster(iCluster); + if (clust->IsPHOS()) hEPHOS->Fill(clust->E()); + if (clust->IsEMCAL()) hEEMCAL->Fill(clust->E()); + } + } + + } + + // perform checks + if (nGen < checkNGenLow) { + Warning("CheckESD", "low number of generated particles: %d", Int_t(nGen)); + } + + TH1F* hEff = CreateEffHisto(hGen, hRec); + + Info("CheckESD", "%d out of %d tracks reconstructed including %d " + "fake tracks", nRec, nGen, nFake); + if (nGen > 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/LHC09d4/Config.C b/prod/LHC09d4/Config.C new file mode 100644 index 00000000000..20921a46ac1 --- /dev/null +++ b/prod/LHC09d4/Config.C @@ -0,0 +1,557 @@ +// +// Configuration for the first physics production 2008 +// + +// 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.C++") + +#if !defined(__CINT__) || defined(__MAKECINT__) +#include +#include +#include +#include +#include +#include +#include "STEER/AliRunLoader.h" +#include "STEER/AliRun.h" +#include "STEER/AliConfig.h" +#include "PYTHIA6/AliDecayerPythia.h" +#include "PYTHIA6/AliGenPythia.h" +#include "TDPMjet/AliGenDPMjet.h" +#include "STEER/AliMagFCheb.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/AliITSv11Hybrid.h" +#include "TPC/AliTPCv2.h" +#include "TOF/AliTOFv6T0.h" +#include "HMPID/AliHMPIDv3.h" +#include "ZDC/AliZDCv3.h" +#include "TRD/AliTRDv1.h" +#include "TRD/AliTRDgeometry.h" +#include "FMD/AliFMDv1.h" +#include "MUON/AliMUONv1.h" +#include "PHOS/AliPHOSv1.h" +#include "PHOS/AliPHOSSimParam.h" +#include "PMD/AliPMDv1.h" +#include "T0/AliT0v1.h" +#include "EMCAL/AliEMCALv2.h" +#include "ACORDE/AliACORDEv1.h" +#include "VZERO/AliVZEROv7.h" +#endif + + +enum PDC06Proc_t +{ + kPythia6, kPythia6D6T, kPhojet, kRunMax +}; + +const char * pprRunName[] = { + "kPythia6", "kPythia6D6T", "kPhojet" +}; + +enum Mag_t +{ + kNoField, k5kG, kFieldMax +}; + +const char * pprField[] = { + "kNoField", "k5kG" +}; + +//--- Functions --- +class AliGenPythia; +AliGenerator *MbPythia(); +AliGenerator *MbPythiaTuneD6T(); +AliGenerator *MbPhojet(); +void ProcessEnvironmentVars(); + +// Geterator, field, beam energy +static PDC06Proc_t proc = kPhojet; +static Mag_t mag = k5kG; +static Float_t energy = 10000; // energy in CMS +//========================// +// Set Random Number seed // +//========================// +TDatime dt; +static UInt_t seed = dt.Get(); + +// Comment line +static TString comment; + +void Config() +{ + + + // Get settings from environment variables + ProcessEnvironmentVars(); + + gRandom->SetSeed(seed); + cerr<<"Seed for random number generation= "<Load("liblhapdf"); // Parton density functions + gSystem->Load("libEGPythia6"); // TGenerator interface + if (proc != kPythia6D6T) { + gSystem->Load("libpythia6"); // Pythia 6.2 + } else { + gSystem->Load("libqpythia"); // Pythia 6.4 + } + gSystem->Load("libAliPythia6"); // ALICE specific implementations + gSystem->Load("libgeant321"); +#endif + + new TGeant3TGeo("C++ Interface to Geant3"); + + //======================================================================= + // Create the output file + + + AliRunLoader* rl=0x0; + + cout<<"Config.C: Creating Run Loader ..."<Fatal("Config.C","Can not instatiate the Run Loader"); + return; + } + rl->SetCompressionLevel(2); + rl->SetNumberOfEventsPerFile(1000); + gAlice->SetRunLoader(rl); + // gAlice->SetGeometryFromFile("geometry.root"); + // gAlice->SetGeometryFromCDB(); + + // Set the trigger configuration: proton-proton + gAlice->SetTriggerDescriptor("p-p"); + + // + //======================================================================= + // ************* STEERING parameters FOR ALICE SIMULATION ************** + // --- Specify event type to be tracked through the ALICE setup + // --- All positions are in cm, angles in degrees, and P and E in GeV + + + gMC->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(); + decayer->SetForceDecay(kAll); + decayer->Init(); + gMC->SetExternalDecayer(decayer); + + //=========================// + // Generator Configuration // + //=========================// + AliGenerator* gener = 0x0; + + if (proc == kPythia6) { + gener = MbPythia(); + } else if (proc == kPythia6D6T) { + gener = MbPythiaTuneD6T(); + } else if (proc == kPhojet) { + gener = MbPhojet(); + } + + + + // PRIMARY VERTEX + // + gener->SetOrigin(0., 0., 0.); // vertex position + // + // + // Size of the interaction diamond + // Longitudinal + Float_t sigmaz = 5.4 / TMath::Sqrt(2.); // [cm] + if (energy == 900) + sigmaz = 10.5 / TMath::Sqrt(2.); // [cm] + + if (energy == 7000) + sigmaz = 6.3 / TMath::Sqrt(2.); // [cm] + + // + // Transverse + Float_t betast = 10; // beta* [m] + Float_t eps = 3.75e-6; // emittance [m] + Float_t gamma = energy / 2.0 / 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(); + if (proc == kPythia6D6T) { + // F I X + AliPythia::Instance()->Pytune(109); + AliPythia::Instance()->Initialize("CMS","p","p", energy); + // F I X + } + // FIELD + // + AliMagF* field = 0x0; + if (mag == kNoField) { + comment = comment.Append(" | L3 field 0.0 T"); + field = new AliMagF("Maps","Maps", 0., 0., AliMagF::k5kGUniform,AliMagF::kBeamTypepp, energy/2.0); + } else if (mag == k5kG) { + comment = comment.Append(" | L3 field 0.5 T"); + field = new AliMagF("Maps","Maps", -1., -1., AliMagF::k5kG, AliMagF::kBeamTypepp, energy/2.0); + } + + printf("\n \n Comment: %s \n \n", comment.Data()); + + TGeoGlobalMagField::Instance()->SetField(field); + + rl->CdGAFile(); + + 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 = 1; + 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"); + FRAME->SetHoles(1); + } + + 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"); + } + + 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 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,1,7,8,9,16,17 + // 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); + geoTRD->SetSMstatus(16,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"); + } + + if (iPHOS) + { + //=================== PHOS parameters =========================== + + AliPHOS *PHOS = new AliPHOSv1("PHOS", "IHEP"); + //Set simulation parameters different from the default ones. + AliPHOSSimParam* simEmc = AliPHOSSimParam::GetInstance() ; + + // APD noise of warm (+20C) PHOS: + // a2 = a1*(Y1/Y2)*(M1/M2), where a1 = 0.012 is APD noise at -25C, + // Y1 = 4.3 photo-electrons/MeV, Y2 = 1.7 p.e/MeV - light yields at -25C and +20C, + // M1 = 50, M2 = 50 - APD gain factors chosen for t1 = -25C and t2 = +20C, + // Y = MeanLightYield*APDEfficiency. + + Float_t apdNoise = 0.012*2.5; + simEmc->SetAPDNoise(apdNoise); + + //Raw Light Yield at +20C + simEmc->SetMeanLightYield(18800); + + //ADC channel width at +18C. + simEmc->SetADCchannelW(0.0125); + } + + + 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) + { + //=================== ACORDE parameters ============================ + + AliACORDE *ACORDE = new AliACORDEv1("ACORDE", "normal ACORDE"); + } + + if (iVZERO) + { + //=================== ACORDE parameters ============================ + + AliVZERO *VZERO = new AliVZEROv7("VZERO", "normal VZERO"); + } +} +// +// PYTHIA +// + +AliGenerator* MbPythia() +{ + comment = comment.Append(" pp at 14 TeV: Pythia low-pt"); +// +// 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(energy); + + return pythia; +} + +AliGenerator* MbPythiaTuneD6T() +{ + comment = comment.Append(" pp at 14 TeV: Pythia low-pt"); +// +// 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(energy); +// Tune +// 109 D6T : Rick Field's CDF Tune D6T (NB: needs CTEQ6L pdfs externally) +// pythia->SetTune(109); // F I X + pythia->SetStrucFunc(kCTEQ6l); +// + return pythia; +} + +AliGenerator* MbPhojet() +{ + comment = comment.Append(" pp at 14 TeV: Pythia low-pt"); +// +// DPMJET +#if defined(__CINT__) + gSystem->Load("libdpmjet"); // Parton density functions + gSystem->Load("libTDPMjet"); // Parton density functions +#endif + AliGenDPMjet* dpmjet = new AliGenDPMjet(-1); + dpmjet->SetMomentumRange(0, 999999.); + dpmjet->SetThetaRange(0., 180.); + dpmjet->SetYRange(-12.,12.); + dpmjet->SetPtRange(0,1000.); + dpmjet->SetProcess(kDpmMb); + dpmjet->SetEnergyCMS(energy); + + return dpmjet; +} + +void ProcessEnvironmentVars() +{ + // Run type + if (gSystem->Getenv("CONFIG_RUN_TYPE")) { + for (Int_t iRun = 0; iRun < kRunMax; iRun++) { + if (strcmp(gSystem->Getenv("CONFIG_RUN_TYPE"), pprRunName[iRun])==0) { + proc = (PDC06Proc_t)iRun; + cout<<"Run type set to "<Getenv("CONFIG_FIELD")) { + for (Int_t iField = 0; iField < kFieldMax; iField++) { + if (strcmp(gSystem->Getenv("CONFIG_FIELD"), pprField[iField])==0) { + mag = (Mag_t)iField; + cout<<"Field set to "<Getenv("CONFIG_ENERGY")) { + energy = atoi(gSystem->Getenv("CONFIG_ENERGY")); + cout<<"Energy set to "<Getenv("CONFIG_SEED")) { + seed = atoi(gSystem->Getenv("CONFIG_SEED")); + } +} diff --git a/prod/LHC09d4/CreateAODfromESD.C b/prod/LHC09d4/CreateAODfromESD.C new file mode 100644 index 00000000000..9c4148736af --- /dev/null +++ b/prod/LHC09d4/CreateAODfromESD.C @@ -0,0 +1,137 @@ +#if !defined(__CINT__) || defined(__MAKECINT__) +#include +#include +#include "AliAnalysisManager.h" +#include "AliESDInputHandler.h" +#include "AliAODHandler.h" +#include "AliAnalysisTaskESDfilter.h" +#include "AliAnalysisDataContainer.h" +#endif + +void CreateAODfromESD(const char *inFileName = "AliESDs.root", + const char *outFileName = "AliAODs.root", + Bool_t bKineFilter = kTRUE) +{ + + gSystem->Load("libTree"); + gSystem->Load("libGeom"); + gSystem->Load("libPhysics"); + gSystem->Load("libVMC"); + gSystem->Load("libSTEERBase"); + gSystem->Load("libESD"); + gSystem->Load("libAOD"); + + gSystem->Load("libANALYSIS"); + gSystem->Load("libANALYSISalice"); + gSystem->Load("libCORRFW"); + gSystem->Load("libPWG3muon"); + + TChain *chain = new TChain("esdTree"); + // Steering input chain + chain->Add(inFileName); + AliAnalysisManager *mgr = new AliAnalysisManager("ESD to AOD", "Analysis Manager"); + + // Input + AliESDInputHandler* inpHandler = new AliESDInputHandler(); + inpHandler->SetReadTags(); + mgr->SetInputEventHandler (inpHandler); + // Output + AliAODHandler* aodHandler = new AliAODHandler(); + aodHandler->SetOutputFileName(outFileName); + mgr->SetOutputEventHandler(aodHandler); + + // MC Truth + if(bKineFilter){ + AliMCEventHandler* mcHandler = new AliMCEventHandler(); + mgr->SetMCtruthEventHandler(mcHandler); + } + + + // Tasks + // 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); + + // Barrel Tracks + AliAnalysisTaskESDfilter *filter = new AliAnalysisTaskESDfilter("Filter"); + mgr->AddTask(filter); + // Muons + AliAnalysisTaskESDMuonFilter *esdmuonfilter = new AliAnalysisTaskESDMuonFilter("ESD Muon Filter"); + mgr->AddTask(esdmuonfilter); + + // Cuts on primary tracks + AliESDtrackCuts* esdTrackCutsL = new AliESDtrackCuts("AliESDtrackCuts", "Standard"); + esdTrackCutsL->SetMinNClustersTPC(50); + esdTrackCutsL->SetMaxChi2PerClusterTPC(3.5); + esdTrackCutsL->SetMaxCovDiagonalElements(2, 2, 0.5, 0.5, 2); + esdTrackCutsL->SetRequireTPCRefit(kTRUE); + esdTrackCutsL->SetMaxDCAToVertexXY(3.0); + esdTrackCutsL->SetMaxDCAToVertexZ(3.0); + esdTrackCutsL->SetDCAToVertex2D(kTRUE); + esdTrackCutsL->SetRequireSigmaToVertex(kFALSE); + esdTrackCutsL->SetAcceptKinkDaughters(kFALSE); + // ITS stand-alone tracks + AliESDtrackCuts* esdTrackCutsITSsa = new AliESDtrackCuts("AliESDtrackCuts", "ITS stand-alone"); + esdTrackCutsITSsa->SetRequireITSStandAlone(kTRUE); + + AliAnalysisFilter* trackFilter = new AliAnalysisFilter("trackFilter"); + trackFilter->AddCuts(esdTrackCutsL); + trackFilter->AddCuts(esdTrackCutsITSsa); + + // Cuts on V0s + AliESDv0Cuts* esdV0Cuts = new AliESDv0Cuts("AliESDv0Cuts", "Standard pp"); + esdV0Cuts->SetMinRadius(0.2); + esdV0Cuts->SetMaxRadius(200); + esdV0Cuts->SetMinDcaPosToVertex(0.05); + esdV0Cuts->SetMinDcaNegToVertex(0.05); + esdV0Cuts->SetMaxDcaV0Daughters(1.0); + esdV0Cuts->SetMinCosinePointingAngle(0.99); + AliAnalysisFilter* v0Filter = new AliAnalysisFilter("v0Filter"); + v0Filter->AddCuts(esdV0Cuts); + + +// + filter->SetTrackFilter(trackFilter); + filter->SetV0Filter(v0Filter); + + +// Create AOD Tags + AliAnalysisTaskTagCreator* tagTask = new AliAnalysisTaskTagCreator("AOD Tag Creator"); + mgr->AddTask(tagTask); + + // Pipelining + AliAnalysisDataContainer *cinput1 = mgr->GetCommonInputContainer(); + AliAnalysisDataContainer *coutput1 = mgr->GetCommonOutputContainer(); + + + AliAnalysisDataContainer *coutputT + = mgr->CreateContainer("cTag", TTree::Class(), AliAnalysisManager::kOutputContainer, "AOD.tag.root"); + + coutput1->SetSpecialOutput(); + coutputT->SetSpecialOutput(); + + if(bKineFilter) { + mgr->ConnectInput (kinefilter, 0, cinput1 ); + mgr->ConnectOutput (kinefilter, 0, coutput1 ); + } + + mgr->ConnectInput (filter, 0, cinput1 ); + mgr->ConnectOutput(filter, 0, coutput1); + + mgr->ConnectInput (esdmuonfilter, 0, cinput1 ); +// mgr->ConnectOutput(esdmuonfilter, 0, coutput1); + + mgr->ConnectInput (tagTask, 0, cinput1); + mgr->ConnectOutput(tagTask, 1, coutputT); + + // + // Run the analysis + // + mgr->InitAnalysis(); + mgr->PrintStatus(); + mgr->StartAnalysis("local", chain); +} + diff --git a/prod/LHC09d4/JDL b/prod/LHC09d4/JDL new file mode 100644 index 00000000000..32ddb57d274 --- /dev/null +++ b/prod/LHC09d4/JDL @@ -0,0 +1,30 @@ +Executable = "aliroot_new"; +Jobtag={"comment:Early physics (2009 - stage 3) pp, Phojet, 0.0T, 900GeV, ID #131"}; + +Packages={"VO_ALICE@AliRoot::v4-17-Rev-15","VO_ALICE@GEANT3::v1-11-4","VO_ALICE@ROOT::v20091104-1","VO_ALICE@APISCONFIG::V1.0x"}; + +TTL = "72000"; + +Requirements = member(other.GridPartitions,"fastmc"); + +Validationcommand ="/alice/cern.ch/user/a/aliprod/prod2007/configs_pbpb_hijing/validation.sh"; + +InputFile= {"LF:/alice/cern.ch/user/a/aliprod/LHC09d4/CheckESD.C", + "LF:/alice/cern.ch/user/a/aliprod/LHC09d4/Config.C", + "LF:/alice/cern.ch/user/a/aliprod/LHC09d4/rec.C", + "LF:/alice/cern.ch/user/a/aliprod/LHC09d4/sim.C", + "LF:/alice/cern.ch/user/a/aliprod/LHC09d4/simrun.C", + "LF:/alice/cern.ch/user/a/aliprod/LHC09d4/tag.C", + "LF:/alice/cern.ch/user/a/aliprod/LHC09d4/CreateAODfromESD.C"}; + +OutputArchive={"log_archive:*.log,stdout,stderr@ALICE::JINR::SE","root_archive.zip:galice.root,Kinematics.root,TrackRefs.root,AliESDs.root,AliESDfriends.root,AliAODs.root,*QA*.root,ITS.RecPoints.root,AliITSPlaneEffSPDtracklet.root,AliITSPlaneEffSPDtrackletBkg.root,TrackletsMCpred.root,TrackleterSPDHistos.root,Run*.root@ALICE::CERN::ALICEDISK,ALICE::FZK::SE,ALICE::CNAF::SE"}; + +OutputDir="/alice/sim/LHC09d4/$1/#alien_counter_03i#"; + +JDLVariables={"Packages", "OutputDir"}; +GUIDFILE="guid.txt"; + +splitarguments="simrun.C --run $1 --event #alien_counter# --process kPhojet --field kNoField --energy 900" ; +split="production:1-1000"; + +Workdirectorysize={"10000MB"}; diff --git a/prod/LHC09d4/rec.C b/prod/LHC09d4/rec.C new file mode 100644 index 00000000000..f1810f3ac19 --- /dev/null +++ b/prod/LHC09d4/rec.C @@ -0,0 +1,24 @@ +void rec() { + + AliReconstruction reco; + + reco.SetWriteESDfriend(); + reco.SetWriteAlignmentData(); + reco.SetRunPlaneEff(kTRUE); + AliITSRecoParam *itspar = AliITSRecoParam::GetPlaneEffParam(-1); + itspar->SetOptTrackletsPlaneEff(kTRUE); + reco.SetRecoParam("ITS",itspar); + + reco.SetDefaultStorage("alien://Folder=/alice/simulation/2008/v4-15-Release/Residual/"); +// reco.SetSpecificStorage("GRP/GRP/Data", +// Form("local://%s",gSystem->pwd())); + // We store the object in AliEn during the simulation + reco.SetSpecificStorage("GRP/GRP/Data", + "alien://Folder=/alice/simulation/2008/v4-15-Release/Ideal/"); + + TStopwatch timer; + timer.Start(); + reco.Run(); + timer.Stop(); + timer.Print(); +} diff --git a/prod/LHC09d4/sim.C b/prod/LHC09d4/sim.C new file mode 100644 index 00000000000..46720580c86 --- /dev/null +++ b/prod/LHC09d4/sim.C @@ -0,0 +1,20 @@ +void sim(Int_t nev=200) { + + AliSimulation simulator; + simulator.SetMakeSDigits("TRD TOF PHOS HMPID EMCAL MUON FMD ZDC PMD T0 VZERO"); + simulator.SetMakeDigitsFromHits("ITS TPC"); + + // The raw data are not written due to the huge increase of the + // virtual memory in HLT + // simulator.SetWriteRawData("ALL","raw.root",kTRUE); + + simulator.SetDefaultStorage("alien://Folder=/alice/simulation/2008/v4-15-Release/Ideal/"); +// simulator.SetSpecificStorage("GRP/GRP/Data", +// Form("local://%s",gSystem->pwd())); + + TStopwatch timer; + timer.Start(); + simulator.Run(nev); + timer.Stop(); + timer.Print(); +} diff --git a/prod/LHC09d4/simrun.C b/prod/LHC09d4/simrun.C new file mode 100644 index 00000000000..068cb6b07ee Binary files /dev/null and b/prod/LHC09d4/simrun.C differ diff --git a/prod/LHC09d4/tag.C b/prod/LHC09d4/tag.C new file mode 100644 index 00000000000..668d95ae835 --- /dev/null +++ b/prod/LHC09d4/tag.C @@ -0,0 +1,108 @@ +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/LHC09d5/CheckESD.C b/prod/LHC09d5/CheckESD.C new file mode 100644 index 00000000000..1f700d32b40 --- /dev/null +++ b/prod/LHC09d5/CheckESD.C @@ -0,0 +1,693 @@ +#if !defined( __CINT__) || defined(__MAKECINT__) +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#include "AliRunLoader.h" +#include "AliLoader.h" +#include "AliESDEvent.h" +#include "AliESDv0.h" +#include "AliESDcascade.h" +#include "AliESDMuonTrack.h" +#include "AliESDCaloCluster.h" +#include "AliRun.h" +#include "AliStack.h" +#include "AliHeader.h" +#include "AliGenEventHeader.h" +#include "AliPID.h" +#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[AliPID::kSPECIES] = + {kElectron, kMuonMinus, kPiPlus, kKPlus, kProton}; + const char* partName[AliPID::kSPECIES+1] = + {"electron", "muon", "pion", "kaon", "proton", "other"}; + Double_t partFrac[AliPID::kSPECIES] = + {0.01, 0.01, 0.85, 0.10, 0.05}; + Int_t identified[AliPID::kSPECIES+1][AliPID::kSPECIES]; + for (Int_t iGen = 0; iGen < AliPID::kSPECIES+1; iGen++) { + for (Int_t iRec = 0; iRec < AliPID::kSPECIES; 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, 50, "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 = runLoader->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 + hResPtInv->Fill(100. * (TMath::Abs(track->GetSigned1Pt()) - 1./particle->Pt()) * + particle->Pt()); + hResPhi->Fill(1000. * (track->Phi() - particle->Phi())); + hResTheta->Fill(1000. * (track->Theta() - particle->Theta())); + + // PID + if ((track->GetStatus() & AliESDtrack::kESDpid) == 0) continue; + Int_t iGen = 5; + for (Int_t i = 0; i < AliPID::kSPECIES; i++) { + if (TMath::Abs(particle->GetPdgCode()) == partCode[i]) iGen = i; + } + Double_t probability[AliPID::kSPECIES]; + track->GetESDpid(probability); + Double_t pMax = 0; + Int_t iRec = 0; + for (Int_t i = 0; i < AliPID::kSPECIES; 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[AliPID::kSPECIES]; + track->GetIntegratedTimes(time); + if (iGen == iRec) { + hDEdxRight->Fill(particle->P(), track->GetTPCsignal()); + if ((track->GetStatus() & AliESDtrack::kTOFpid) != 0) { + hResTOFRight->Fill(track->GetTOFsignal() - time[iRec]); + } + } else { + hDEdxWrong->Fill(particle->P(), 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->GetEffMassXi()); + cascade->ChangeMassHypothesis(v0q,kOmegaMinus); + hMassOmega->Fill(cascade->GetEffMassXi()); + + 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 clusters + { + for (Int_t iCluster=0; iClusterGetNumberOfCaloClusters(); iCluster++) { + AliESDCaloCluster * clust = esd->GetCaloCluster(iCluster); + if (clust->IsPHOS()) hEPHOS->Fill(clust->E()); + if (clust->IsEMCAL()) hEEMCAL->Fill(clust->E()); + } + } + + } + + // perform checks + if (nGen < checkNGenLow) { + Warning("CheckESD", "low number of generated particles: %d", Int_t(nGen)); + } + + TH1F* hEff = CreateEffHisto(hGen, hRec); + + Info("CheckESD", "%d out of %d tracks reconstructed including %d " + "fake tracks", nRec, nGen, nFake); + if (nGen > 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/LHC09d5/Config.C b/prod/LHC09d5/Config.C new file mode 100644 index 00000000000..b9cf8e85a9b --- /dev/null +++ b/prod/LHC09d5/Config.C @@ -0,0 +1,557 @@ +// +// Configuration for the first physics production 2008 +// + +// 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.C++") + +#if !defined(__CINT__) || defined(__MAKECINT__) +#include +#include +#include +#include +#include +#include +#include "STEER/AliRunLoader.h" +#include "STEER/AliRun.h" +#include "STEER/AliConfig.h" +#include "PYTHIA6/AliDecayerPythia.h" +#include "PYTHIA6/AliGenPythia.h" +#include "TDPMjet/AliGenDPMjet.h" +#include "STEER/AliMagFCheb.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/AliITSv11Hybrid.h" +#include "TPC/AliTPCv2.h" +#include "TOF/AliTOFv6T0.h" +#include "HMPID/AliHMPIDv3.h" +#include "ZDC/AliZDCv3.h" +#include "TRD/AliTRDv1.h" +#include "TRD/AliTRDgeometry.h" +#include "FMD/AliFMDv1.h" +#include "MUON/AliMUONv1.h" +#include "PHOS/AliPHOSv1.h" +#include "PHOS/AliPHOSSimParam.h" +#include "PMD/AliPMDv1.h" +#include "T0/AliT0v1.h" +#include "EMCAL/AliEMCALv2.h" +#include "ACORDE/AliACORDEv1.h" +#include "VZERO/AliVZEROv7.h" +#endif + + +enum PDC06Proc_t +{ + kPythia6, kPythia6D6T, kPhojet, kRunMax +}; + +const char * pprRunName[] = { + "kPythia6", "kPythia6D6T", "kPhojet" +}; + +enum Mag_t +{ + kNoField, k5kG, kFieldMax +}; + +const char * pprField[] = { + "kNoField", "k5kG" +}; + +//--- Functions --- +class AliGenPythia; +AliGenerator *MbPythia(); +AliGenerator *MbPythiaTuneD6T(); +AliGenerator *MbPhojet(); +void ProcessEnvironmentVars(); + +// Geterator, field, beam energy +static PDC06Proc_t proc = kPhojet; +static Mag_t mag = k5kG; +static Float_t energy = 10000; // energy in CMS +//========================// +// Set Random Number seed // +//========================// +TDatime dt; +static UInt_t seed = dt.Get(); + +// Comment line +static TString comment; + +void Config() +{ + + + // Get settings from environment variables + ProcessEnvironmentVars(); + + gRandom->SetSeed(seed); + cerr<<"Seed for random number generation= "<Load("liblhapdf"); // Parton density functions + gSystem->Load("libEGPythia6"); // TGenerator interface + if (proc != kPythia6D6T) { + gSystem->Load("libpythia6"); // Pythia 6.2 + } else { + gSystem->Load("libqpythia"); // Pythia 6.4 + } + gSystem->Load("libAliPythia6"); // ALICE specific implementations + gSystem->Load("libgeant321"); +#endif + + new TGeant3TGeo("C++ Interface to Geant3"); + + //======================================================================= + // Create the output file + + + AliRunLoader* rl=0x0; + + cout<<"Config.C: Creating Run Loader ..."<Fatal("Config.C","Can not instatiate the Run Loader"); + return; + } + rl->SetCompressionLevel(2); + rl->SetNumberOfEventsPerFile(1000); + gAlice->SetRunLoader(rl); + // gAlice->SetGeometryFromFile("geometry.root"); + // gAlice->SetGeometryFromCDB(); + + // Set the trigger configuration: proton-proton + gAlice->SetTriggerDescriptor("p-p"); + + // + //======================================================================= + // ************* STEERING parameters FOR ALICE SIMULATION ************** + // --- Specify event type to be tracked through the ALICE setup + // --- All positions are in cm, angles in degrees, and P and E in GeV + + + gMC->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(); + decayer->SetForceDecay(kAll); + decayer->Init(); + gMC->SetExternalDecayer(decayer); + + //=========================// + // Generator Configuration // + //=========================// + AliGenerator* gener = 0x0; + + if (proc == kPythia6) { + gener = MbPythia(); + } else if (proc == kPythia6D6T) { + gener = MbPythiaTuneD6T(); + } else if (proc == kPhojet) { + gener = MbPhojet(); + } + + + + // PRIMARY VERTEX + // + gener->SetOrigin(-0.07, 0.25, -1.5); // vertex position + // + // + // Size of the interaction diamond + // Longitudinal + Float_t sigmaz = 5.4 / TMath::Sqrt(2.); // [cm] + if (energy == 900) + sigmaz = 10.5 / TMath::Sqrt(2.); // [cm] + + if (energy == 7000) + sigmaz = 6.3 / TMath::Sqrt(2.); // [cm] + + // + // Transverse + Float_t betast = 10; // beta* [m] + Float_t eps = 3.75e-6; // emittance [m] + Float_t gamma = energy / 2.0 / 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(); + if (proc == kPythia6D6T) { + // F I X + AliPythia::Instance()->Pytune(109); + AliPythia::Instance()->Initialize("CMS","p","p", energy); + // F I X + } + // FIELD + // + AliMagF* field = 0x0; + if (mag == kNoField) { + comment = comment.Append(" | L3 field 0.0 T"); + field = new AliMagF("Maps","Maps", 0., 0., AliMagF::k5kGUniform,AliMagF::kBeamTypepp, energy/2.0); + } else if (mag == k5kG) { + comment = comment.Append(" | L3 field 0.5 T"); + field = new AliMagF("Maps","Maps", -1., -1., AliMagF::k5kG, AliMagF::kBeamTypepp, energy/2.0); + } + + printf("\n \n Comment: %s \n \n", comment.Data()); + + TGeoGlobalMagField::Instance()->SetField(field); + + rl->CdGAFile(); + + 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 = 1; + 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"); + FRAME->SetHoles(1); + } + + 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"); + } + + 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 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,1,7,8,9,16,17 + // 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(10,0); + geoTRD->SetSMstatus(11,0); + geoTRD->SetSMstatus(12,0); + 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"); + } + + if (iPHOS) + { + //=================== PHOS parameters =========================== + + AliPHOS *PHOS = new AliPHOSv1("PHOS", "IHEP"); + //Set simulation parameters different from the default ones. + AliPHOSSimParam* simEmc = AliPHOSSimParam::GetInstance() ; + + // APD noise of warm (+20C) PHOS: + // a2 = a1*(Y1/Y2)*(M1/M2), where a1 = 0.012 is APD noise at -25C, + // Y1 = 4.3 photo-electrons/MeV, Y2 = 1.7 p.e/MeV - light yields at -25C and +20C, + // M1 = 50, M2 = 50 - APD gain factors chosen for t1 = -25C and t2 = +20C, + // Y = MeanLightYield*APDEfficiency. + + Float_t apdNoise = 0.012*2.5; + simEmc->SetAPDNoise(apdNoise); + + //Raw Light Yield at +20C + simEmc->SetMeanLightYield(18800); + + //ADC channel width at +18C. + simEmc->SetADCchannelW(0.0125); + } + + + 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) + { + //=================== ACORDE parameters ============================ + + AliACORDE *ACORDE = new AliACORDEv1("ACORDE", "normal ACORDE"); + } + + if (iVZERO) + { + //=================== ACORDE parameters ============================ + + AliVZERO *VZERO = new AliVZEROv7("VZERO", "normal VZERO"); + } +} +// +// PYTHIA +// + +AliGenerator* MbPythia() +{ + comment = comment.Append(" pp at 14 TeV: Pythia low-pt"); +// +// 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(energy); + + return pythia; +} + +AliGenerator* MbPythiaTuneD6T() +{ + comment = comment.Append(" pp at 14 TeV: Pythia low-pt"); +// +// 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(energy); +// Tune +// 109 D6T : Rick Field's CDF Tune D6T (NB: needs CTEQ6L pdfs externally) +// pythia->SetTune(109); // F I X + pythia->SetStrucFunc(kCTEQ6l); +// + return pythia; +} + +AliGenerator* MbPhojet() +{ + comment = comment.Append(" pp at 14 TeV: Pythia low-pt"); +// +// DPMJET +#if defined(__CINT__) + gSystem->Load("libdpmjet"); // Parton density functions + gSystem->Load("libTDPMjet"); // Parton density functions +#endif + AliGenDPMjet* dpmjet = new AliGenDPMjet(-1); + dpmjet->SetMomentumRange(0, 999999.); + dpmjet->SetThetaRange(0., 180.); + dpmjet->SetYRange(-12.,12.); + dpmjet->SetPtRange(0,1000.); + dpmjet->SetProcess(kDpmMb); + dpmjet->SetEnergyCMS(energy); + + return dpmjet; +} + +void ProcessEnvironmentVars() +{ + // Run type + if (gSystem->Getenv("CONFIG_RUN_TYPE")) { + for (Int_t iRun = 0; iRun < kRunMax; iRun++) { + if (strcmp(gSystem->Getenv("CONFIG_RUN_TYPE"), pprRunName[iRun])==0) { + proc = (PDC06Proc_t)iRun; + cout<<"Run type set to "<Getenv("CONFIG_FIELD")) { + for (Int_t iField = 0; iField < kFieldMax; iField++) { + if (strcmp(gSystem->Getenv("CONFIG_FIELD"), pprField[iField])==0) { + mag = (Mag_t)iField; + cout<<"Field set to "<Getenv("CONFIG_ENERGY")) { + energy = atoi(gSystem->Getenv("CONFIG_ENERGY")); + cout<<"Energy set to "<Getenv("CONFIG_SEED")) { + seed = atoi(gSystem->Getenv("CONFIG_SEED")); + } +} diff --git a/prod/LHC09d5/CreateAODfromESD.C b/prod/LHC09d5/CreateAODfromESD.C new file mode 100644 index 00000000000..9c4148736af --- /dev/null +++ b/prod/LHC09d5/CreateAODfromESD.C @@ -0,0 +1,137 @@ +#if !defined(__CINT__) || defined(__MAKECINT__) +#include +#include +#include "AliAnalysisManager.h" +#include "AliESDInputHandler.h" +#include "AliAODHandler.h" +#include "AliAnalysisTaskESDfilter.h" +#include "AliAnalysisDataContainer.h" +#endif + +void CreateAODfromESD(const char *inFileName = "AliESDs.root", + const char *outFileName = "AliAODs.root", + Bool_t bKineFilter = kTRUE) +{ + + gSystem->Load("libTree"); + gSystem->Load("libGeom"); + gSystem->Load("libPhysics"); + gSystem->Load("libVMC"); + gSystem->Load("libSTEERBase"); + gSystem->Load("libESD"); + gSystem->Load("libAOD"); + + gSystem->Load("libANALYSIS"); + gSystem->Load("libANALYSISalice"); + gSystem->Load("libCORRFW"); + gSystem->Load("libPWG3muon"); + + TChain *chain = new TChain("esdTree"); + // Steering input chain + chain->Add(inFileName); + AliAnalysisManager *mgr = new AliAnalysisManager("ESD to AOD", "Analysis Manager"); + + // Input + AliESDInputHandler* inpHandler = new AliESDInputHandler(); + inpHandler->SetReadTags(); + mgr->SetInputEventHandler (inpHandler); + // Output + AliAODHandler* aodHandler = new AliAODHandler(); + aodHandler->SetOutputFileName(outFileName); + mgr->SetOutputEventHandler(aodHandler); + + // MC Truth + if(bKineFilter){ + AliMCEventHandler* mcHandler = new AliMCEventHandler(); + mgr->SetMCtruthEventHandler(mcHandler); + } + + + // Tasks + // 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); + + // Barrel Tracks + AliAnalysisTaskESDfilter *filter = new AliAnalysisTaskESDfilter("Filter"); + mgr->AddTask(filter); + // Muons + AliAnalysisTaskESDMuonFilter *esdmuonfilter = new AliAnalysisTaskESDMuonFilter("ESD Muon Filter"); + mgr->AddTask(esdmuonfilter); + + // Cuts on primary tracks + AliESDtrackCuts* esdTrackCutsL = new AliESDtrackCuts("AliESDtrackCuts", "Standard"); + esdTrackCutsL->SetMinNClustersTPC(50); + esdTrackCutsL->SetMaxChi2PerClusterTPC(3.5); + esdTrackCutsL->SetMaxCovDiagonalElements(2, 2, 0.5, 0.5, 2); + esdTrackCutsL->SetRequireTPCRefit(kTRUE); + esdTrackCutsL->SetMaxDCAToVertexXY(3.0); + esdTrackCutsL->SetMaxDCAToVertexZ(3.0); + esdTrackCutsL->SetDCAToVertex2D(kTRUE); + esdTrackCutsL->SetRequireSigmaToVertex(kFALSE); + esdTrackCutsL->SetAcceptKinkDaughters(kFALSE); + // ITS stand-alone tracks + AliESDtrackCuts* esdTrackCutsITSsa = new AliESDtrackCuts("AliESDtrackCuts", "ITS stand-alone"); + esdTrackCutsITSsa->SetRequireITSStandAlone(kTRUE); + + AliAnalysisFilter* trackFilter = new AliAnalysisFilter("trackFilter"); + trackFilter->AddCuts(esdTrackCutsL); + trackFilter->AddCuts(esdTrackCutsITSsa); + + // Cuts on V0s + AliESDv0Cuts* esdV0Cuts = new AliESDv0Cuts("AliESDv0Cuts", "Standard pp"); + esdV0Cuts->SetMinRadius(0.2); + esdV0Cuts->SetMaxRadius(200); + esdV0Cuts->SetMinDcaPosToVertex(0.05); + esdV0Cuts->SetMinDcaNegToVertex(0.05); + esdV0Cuts->SetMaxDcaV0Daughters(1.0); + esdV0Cuts->SetMinCosinePointingAngle(0.99); + AliAnalysisFilter* v0Filter = new AliAnalysisFilter("v0Filter"); + v0Filter->AddCuts(esdV0Cuts); + + +// + filter->SetTrackFilter(trackFilter); + filter->SetV0Filter(v0Filter); + + +// Create AOD Tags + AliAnalysisTaskTagCreator* tagTask = new AliAnalysisTaskTagCreator("AOD Tag Creator"); + mgr->AddTask(tagTask); + + // Pipelining + AliAnalysisDataContainer *cinput1 = mgr->GetCommonInputContainer(); + AliAnalysisDataContainer *coutput1 = mgr->GetCommonOutputContainer(); + + + AliAnalysisDataContainer *coutputT + = mgr->CreateContainer("cTag", TTree::Class(), AliAnalysisManager::kOutputContainer, "AOD.tag.root"); + + coutput1->SetSpecialOutput(); + coutputT->SetSpecialOutput(); + + if(bKineFilter) { + mgr->ConnectInput (kinefilter, 0, cinput1 ); + mgr->ConnectOutput (kinefilter, 0, coutput1 ); + } + + mgr->ConnectInput (filter, 0, cinput1 ); + mgr->ConnectOutput(filter, 0, coutput1); + + mgr->ConnectInput (esdmuonfilter, 0, cinput1 ); +// mgr->ConnectOutput(esdmuonfilter, 0, coutput1); + + mgr->ConnectInput (tagTask, 0, cinput1); + mgr->ConnectOutput(tagTask, 1, coutputT); + + // + // Run the analysis + // + mgr->InitAnalysis(); + mgr->PrintStatus(); + mgr->StartAnalysis("local", chain); +} + diff --git a/prod/LHC09d5/JDL b/prod/LHC09d5/JDL new file mode 100644 index 00000000000..db7ef4a8812 --- /dev/null +++ b/prod/LHC09d5/JDL @@ -0,0 +1,31 @@ +Executable = "aliroot_new"; +Jobtag={"comment:Early physics (2009 - FIRST DATA) pp, Pythia6 Tune D6T, 0T, 900GeV, ID #132"}; + +Packages={"VO_ALICE@AliRoot::v4-17-Rev-16","VO_ALICE@GEANT3::v1-11-4","VO_ALICE@ROOT::v20091109-1","VO_ALICE@APISCONFIG::V1.0x"}; + +TTL = "72000"; + +Requirements = member(other.GridPartitions,"fastmc"); + +Validationcommand ="/alice/cern.ch/user/a/aliprod/prod2007/configs_pbpb_hijing/validation.sh"; + +InputFile= {"LF:/alice/cern.ch/user/a/aliprod/LHC09d5/CheckESD.C", + "LF:/alice/cern.ch/user/a/aliprod/LHC09d5/Config.C", + "LF:/alice/cern.ch/user/a/aliprod/LHC09d5/rec.C", + "LF:/alice/cern.ch/user/a/aliprod/LHC09d5/sim.C", + "LF:/alice/cern.ch/user/a/aliprod/LHC09d5/simrun.C", + "LF:/alice/cern.ch/user/a/aliprod/LHC09d5/tag.C", + "LF:/alice/cern.ch/user/a/aliprod/LHC09d5/CreateAODfromESD.C"}; + + +OutputArchive={"log_archive:*.log,stdout,stderr@ALICE::JINR::SE","root_archive.zip:galice.root,Kinematics.root,TrackRefs.root,AliESDs.root,AliESDfriends.root,AliAODs.root,*QA*.root,ITS.RecPoints.root,AliITSPlaneEffSPDtracklet.root,AliITSPlaneEffSPDtrackletBkg.root,TrackletsMCpred.root,TrackleterSPDHistos.root,Run*.root@ALICE::CERN::ALICEDISK,ALICE::FZK::SE,ALICE::CNAF::SE"}; + +OutputDir="/alice/sim/LHC09d5/$1/#alien_counter_03i#"; + +JDLVariables={"Packages", "OutputDir"}; +GUIDFILE="guid.txt"; + +splitarguments="simrun.C --run $1 --event #alien_counter# --process kPythia6D6T --field kNoField --energy 900" ; +split="production:1-1000"; + +Workdirectorysize={"10000MB"}; diff --git a/prod/LHC09d5/rec.C b/prod/LHC09d5/rec.C new file mode 100644 index 00000000000..d883ceee088 --- /dev/null +++ b/prod/LHC09d5/rec.C @@ -0,0 +1,53 @@ +void rec() { + + AliReconstruction reco; + + reco.SetWriteESDfriend(); + reco.SetWriteAlignmentData(); + reco.SetRunReconstruction("ITS VZERO"); + + // ITS Efficiency + reco.SetRunPlaneEff(kTRUE); + AliITSRecoParam *itsRecoParam = AliITSRecoParam::GetPlaneEffParam(-1); + itsRecoParam->SetOptTrackletsPlaneEff(kTRUE); + + + //****** FIRST PHYSICS 2009 (same as COSMICS 2009) ********************* + //AliITSRecoParam * itsRecoParam = AliITSRecoParam::GetLowFluxParam(); + itsRecoParam->SetClusterErrorsParam(2); + // find independently ITS SA tracks + itsRecoParam->SetSAUseAllClusters(); + itsRecoParam->SetOuterStartLayerSA(AliITSgeomTGeo::GetNLayers()-2); + // to maximize efficiency + itsRecoParam->SetAllowProlongationWithEmptyRoad(); // larger seach windows for SA (in case of large misalignments) + itsRecoParam->SetNLoopsSA(33); + itsRecoParam->SetFactorSAWindowSizes(20); + // additional error due to misal (B off) + itsRecoParam->SetClusterMisalErrorY(1.0,1.0,1.0,1.0,1.0,1.0); // [cm] + itsRecoParam->SetClusterMisalErrorZ(1.0,1.0,1.0,1.0,1.0,1.0); // [cm] + // additional error due to misal (B on) + itsRecoParam->SetClusterMisalErrorYBOn(0.0,0.0,0.1,0.1,0.1,0.1); // [cm] + itsRecoParam->SetClusterMisalErrorZBOn(0.1,0.1,0.1,0.1,0.1,0.1); // [cm] + //*********************************************************************** + itsRecoParam->SetEventSpecie(AliRecoParam::kLowMult); + + reco.SetRecoParam("ITS",itsRecoParam); + + + + reco.SetDefaultStorage("alien://Folder=/alice/simulation/2008/v4-15-Release/Residual/"); +// reco.SetSpecificStorage("GRP/GRP/Data", +// Form("local://%s",gSystem->pwd())); + // We store the object in AliEn during the simulation + reco.SetSpecificStorage("GRP/GRP/Data", + "alien://Folder=/alice/simulation/2008/v4-15-Release/Ideal/"); + + + reco.SetSpecificStorage("GRP/Calib/MeanVertexSPD","alien://folder=/alice/cern.ch/user/r/rgrosso/ShiftedVertex"); + + TStopwatch timer; + timer.Start(); + reco.Run(); + timer.Stop(); + timer.Print(); +} diff --git a/prod/LHC09d5/sim.C b/prod/LHC09d5/sim.C new file mode 100644 index 00000000000..ebc41521069 --- /dev/null +++ b/prod/LHC09d5/sim.C @@ -0,0 +1,20 @@ +void sim(Int_t nev=200) { + AliSimulation simulator; + simulator.SetRunHLT(""); + simulator.SetMakeSDigits("VZERO"); + simulator.SetMakeDigitsFromHits("ITS"); + simulator.SetMakeDigits("ITS VZERO"); + +// simulator.SetWriteRawData("ALL","raw.root",kTRUE); + + simulator.SetDefaultStorage("alien://Folder=/alice/simulation/2008/v4-15-Release/Ideal/"); + +// simulator.SetSpecificStorage("GRP/GRP/Data", +// Form("local://%s",gSystem->pwd())); + + TStopwatch timer; + timer.Start(); + simulator.Run(nev); + timer.Stop(); + timer.Print(); +} diff --git a/prod/LHC09d5/simrun.C b/prod/LHC09d5/simrun.C new file mode 100644 index 00000000000..068cb6b07ee Binary files /dev/null and b/prod/LHC09d5/simrun.C differ diff --git a/prod/LHC09d5/tag.C b/prod/LHC09d5/tag.C new file mode 100644 index 00000000000..668d95ae835 --- /dev/null +++ b/prod/LHC09d5/tag.C @@ -0,0 +1,108 @@ +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/LHC09d6/CheckESD.C b/prod/LHC09d6/CheckESD.C new file mode 100644 index 00000000000..1f700d32b40 --- /dev/null +++ b/prod/LHC09d6/CheckESD.C @@ -0,0 +1,693 @@ +#if !defined( __CINT__) || defined(__MAKECINT__) +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#include "AliRunLoader.h" +#include "AliLoader.h" +#include "AliESDEvent.h" +#include "AliESDv0.h" +#include "AliESDcascade.h" +#include "AliESDMuonTrack.h" +#include "AliESDCaloCluster.h" +#include "AliRun.h" +#include "AliStack.h" +#include "AliHeader.h" +#include "AliGenEventHeader.h" +#include "AliPID.h" +#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[AliPID::kSPECIES] = + {kElectron, kMuonMinus, kPiPlus, kKPlus, kProton}; + const char* partName[AliPID::kSPECIES+1] = + {"electron", "muon", "pion", "kaon", "proton", "other"}; + Double_t partFrac[AliPID::kSPECIES] = + {0.01, 0.01, 0.85, 0.10, 0.05}; + Int_t identified[AliPID::kSPECIES+1][AliPID::kSPECIES]; + for (Int_t iGen = 0; iGen < AliPID::kSPECIES+1; iGen++) { + for (Int_t iRec = 0; iRec < AliPID::kSPECIES; 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, 50, "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 = runLoader->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 + hResPtInv->Fill(100. * (TMath::Abs(track->GetSigned1Pt()) - 1./particle->Pt()) * + particle->Pt()); + hResPhi->Fill(1000. * (track->Phi() - particle->Phi())); + hResTheta->Fill(1000. * (track->Theta() - particle->Theta())); + + // PID + if ((track->GetStatus() & AliESDtrack::kESDpid) == 0) continue; + Int_t iGen = 5; + for (Int_t i = 0; i < AliPID::kSPECIES; i++) { + if (TMath::Abs(particle->GetPdgCode()) == partCode[i]) iGen = i; + } + Double_t probability[AliPID::kSPECIES]; + track->GetESDpid(probability); + Double_t pMax = 0; + Int_t iRec = 0; + for (Int_t i = 0; i < AliPID::kSPECIES; 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[AliPID::kSPECIES]; + track->GetIntegratedTimes(time); + if (iGen == iRec) { + hDEdxRight->Fill(particle->P(), track->GetTPCsignal()); + if ((track->GetStatus() & AliESDtrack::kTOFpid) != 0) { + hResTOFRight->Fill(track->GetTOFsignal() - time[iRec]); + } + } else { + hDEdxWrong->Fill(particle->P(), 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->GetEffMassXi()); + cascade->ChangeMassHypothesis(v0q,kOmegaMinus); + hMassOmega->Fill(cascade->GetEffMassXi()); + + 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 clusters + { + for (Int_t iCluster=0; iClusterGetNumberOfCaloClusters(); iCluster++) { + AliESDCaloCluster * clust = esd->GetCaloCluster(iCluster); + if (clust->IsPHOS()) hEPHOS->Fill(clust->E()); + if (clust->IsEMCAL()) hEEMCAL->Fill(clust->E()); + } + } + + } + + // perform checks + if (nGen < checkNGenLow) { + Warning("CheckESD", "low number of generated particles: %d", Int_t(nGen)); + } + + TH1F* hEff = CreateEffHisto(hGen, hRec); + + Info("CheckESD", "%d out of %d tracks reconstructed including %d " + "fake tracks", nRec, nGen, nFake); + if (nGen > 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/LHC09d6/Config.C b/prod/LHC09d6/Config.C new file mode 100644 index 00000000000..b9cf8e85a9b --- /dev/null +++ b/prod/LHC09d6/Config.C @@ -0,0 +1,557 @@ +// +// Configuration for the first physics production 2008 +// + +// 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.C++") + +#if !defined(__CINT__) || defined(__MAKECINT__) +#include +#include +#include +#include +#include +#include +#include "STEER/AliRunLoader.h" +#include "STEER/AliRun.h" +#include "STEER/AliConfig.h" +#include "PYTHIA6/AliDecayerPythia.h" +#include "PYTHIA6/AliGenPythia.h" +#include "TDPMjet/AliGenDPMjet.h" +#include "STEER/AliMagFCheb.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/AliITSv11Hybrid.h" +#include "TPC/AliTPCv2.h" +#include "TOF/AliTOFv6T0.h" +#include "HMPID/AliHMPIDv3.h" +#include "ZDC/AliZDCv3.h" +#include "TRD/AliTRDv1.h" +#include "TRD/AliTRDgeometry.h" +#include "FMD/AliFMDv1.h" +#include "MUON/AliMUONv1.h" +#include "PHOS/AliPHOSv1.h" +#include "PHOS/AliPHOSSimParam.h" +#include "PMD/AliPMDv1.h" +#include "T0/AliT0v1.h" +#include "EMCAL/AliEMCALv2.h" +#include "ACORDE/AliACORDEv1.h" +#include "VZERO/AliVZEROv7.h" +#endif + + +enum PDC06Proc_t +{ + kPythia6, kPythia6D6T, kPhojet, kRunMax +}; + +const char * pprRunName[] = { + "kPythia6", "kPythia6D6T", "kPhojet" +}; + +enum Mag_t +{ + kNoField, k5kG, kFieldMax +}; + +const char * pprField[] = { + "kNoField", "k5kG" +}; + +//--- Functions --- +class AliGenPythia; +AliGenerator *MbPythia(); +AliGenerator *MbPythiaTuneD6T(); +AliGenerator *MbPhojet(); +void ProcessEnvironmentVars(); + +// Geterator, field, beam energy +static PDC06Proc_t proc = kPhojet; +static Mag_t mag = k5kG; +static Float_t energy = 10000; // energy in CMS +//========================// +// Set Random Number seed // +//========================// +TDatime dt; +static UInt_t seed = dt.Get(); + +// Comment line +static TString comment; + +void Config() +{ + + + // Get settings from environment variables + ProcessEnvironmentVars(); + + gRandom->SetSeed(seed); + cerr<<"Seed for random number generation= "<Load("liblhapdf"); // Parton density functions + gSystem->Load("libEGPythia6"); // TGenerator interface + if (proc != kPythia6D6T) { + gSystem->Load("libpythia6"); // Pythia 6.2 + } else { + gSystem->Load("libqpythia"); // Pythia 6.4 + } + gSystem->Load("libAliPythia6"); // ALICE specific implementations + gSystem->Load("libgeant321"); +#endif + + new TGeant3TGeo("C++ Interface to Geant3"); + + //======================================================================= + // Create the output file + + + AliRunLoader* rl=0x0; + + cout<<"Config.C: Creating Run Loader ..."<Fatal("Config.C","Can not instatiate the Run Loader"); + return; + } + rl->SetCompressionLevel(2); + rl->SetNumberOfEventsPerFile(1000); + gAlice->SetRunLoader(rl); + // gAlice->SetGeometryFromFile("geometry.root"); + // gAlice->SetGeometryFromCDB(); + + // Set the trigger configuration: proton-proton + gAlice->SetTriggerDescriptor("p-p"); + + // + //======================================================================= + // ************* STEERING parameters FOR ALICE SIMULATION ************** + // --- Specify event type to be tracked through the ALICE setup + // --- All positions are in cm, angles in degrees, and P and E in GeV + + + gMC->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(); + decayer->SetForceDecay(kAll); + decayer->Init(); + gMC->SetExternalDecayer(decayer); + + //=========================// + // Generator Configuration // + //=========================// + AliGenerator* gener = 0x0; + + if (proc == kPythia6) { + gener = MbPythia(); + } else if (proc == kPythia6D6T) { + gener = MbPythiaTuneD6T(); + } else if (proc == kPhojet) { + gener = MbPhojet(); + } + + + + // PRIMARY VERTEX + // + gener->SetOrigin(-0.07, 0.25, -1.5); // vertex position + // + // + // Size of the interaction diamond + // Longitudinal + Float_t sigmaz = 5.4 / TMath::Sqrt(2.); // [cm] + if (energy == 900) + sigmaz = 10.5 / TMath::Sqrt(2.); // [cm] + + if (energy == 7000) + sigmaz = 6.3 / TMath::Sqrt(2.); // [cm] + + // + // Transverse + Float_t betast = 10; // beta* [m] + Float_t eps = 3.75e-6; // emittance [m] + Float_t gamma = energy / 2.0 / 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(); + if (proc == kPythia6D6T) { + // F I X + AliPythia::Instance()->Pytune(109); + AliPythia::Instance()->Initialize("CMS","p","p", energy); + // F I X + } + // FIELD + // + AliMagF* field = 0x0; + if (mag == kNoField) { + comment = comment.Append(" | L3 field 0.0 T"); + field = new AliMagF("Maps","Maps", 0., 0., AliMagF::k5kGUniform,AliMagF::kBeamTypepp, energy/2.0); + } else if (mag == k5kG) { + comment = comment.Append(" | L3 field 0.5 T"); + field = new AliMagF("Maps","Maps", -1., -1., AliMagF::k5kG, AliMagF::kBeamTypepp, energy/2.0); + } + + printf("\n \n Comment: %s \n \n", comment.Data()); + + TGeoGlobalMagField::Instance()->SetField(field); + + rl->CdGAFile(); + + 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 = 1; + 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"); + FRAME->SetHoles(1); + } + + 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"); + } + + 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 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,1,7,8,9,16,17 + // 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(10,0); + geoTRD->SetSMstatus(11,0); + geoTRD->SetSMstatus(12,0); + 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"); + } + + if (iPHOS) + { + //=================== PHOS parameters =========================== + + AliPHOS *PHOS = new AliPHOSv1("PHOS", "IHEP"); + //Set simulation parameters different from the default ones. + AliPHOSSimParam* simEmc = AliPHOSSimParam::GetInstance() ; + + // APD noise of warm (+20C) PHOS: + // a2 = a1*(Y1/Y2)*(M1/M2), where a1 = 0.012 is APD noise at -25C, + // Y1 = 4.3 photo-electrons/MeV, Y2 = 1.7 p.e/MeV - light yields at -25C and +20C, + // M1 = 50, M2 = 50 - APD gain factors chosen for t1 = -25C and t2 = +20C, + // Y = MeanLightYield*APDEfficiency. + + Float_t apdNoise = 0.012*2.5; + simEmc->SetAPDNoise(apdNoise); + + //Raw Light Yield at +20C + simEmc->SetMeanLightYield(18800); + + //ADC channel width at +18C. + simEmc->SetADCchannelW(0.0125); + } + + + 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) + { + //=================== ACORDE parameters ============================ + + AliACORDE *ACORDE = new AliACORDEv1("ACORDE", "normal ACORDE"); + } + + if (iVZERO) + { + //=================== ACORDE parameters ============================ + + AliVZERO *VZERO = new AliVZEROv7("VZERO", "normal VZERO"); + } +} +// +// PYTHIA +// + +AliGenerator* MbPythia() +{ + comment = comment.Append(" pp at 14 TeV: Pythia low-pt"); +// +// 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(energy); + + return pythia; +} + +AliGenerator* MbPythiaTuneD6T() +{ + comment = comment.Append(" pp at 14 TeV: Pythia low-pt"); +// +// 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(energy); +// Tune +// 109 D6T : Rick Field's CDF Tune D6T (NB: needs CTEQ6L pdfs externally) +// pythia->SetTune(109); // F I X + pythia->SetStrucFunc(kCTEQ6l); +// + return pythia; +} + +AliGenerator* MbPhojet() +{ + comment = comment.Append(" pp at 14 TeV: Pythia low-pt"); +// +// DPMJET +#if defined(__CINT__) + gSystem->Load("libdpmjet"); // Parton density functions + gSystem->Load("libTDPMjet"); // Parton density functions +#endif + AliGenDPMjet* dpmjet = new AliGenDPMjet(-1); + dpmjet->SetMomentumRange(0, 999999.); + dpmjet->SetThetaRange(0., 180.); + dpmjet->SetYRange(-12.,12.); + dpmjet->SetPtRange(0,1000.); + dpmjet->SetProcess(kDpmMb); + dpmjet->SetEnergyCMS(energy); + + return dpmjet; +} + +void ProcessEnvironmentVars() +{ + // Run type + if (gSystem->Getenv("CONFIG_RUN_TYPE")) { + for (Int_t iRun = 0; iRun < kRunMax; iRun++) { + if (strcmp(gSystem->Getenv("CONFIG_RUN_TYPE"), pprRunName[iRun])==0) { + proc = (PDC06Proc_t)iRun; + cout<<"Run type set to "<Getenv("CONFIG_FIELD")) { + for (Int_t iField = 0; iField < kFieldMax; iField++) { + if (strcmp(gSystem->Getenv("CONFIG_FIELD"), pprField[iField])==0) { + mag = (Mag_t)iField; + cout<<"Field set to "<Getenv("CONFIG_ENERGY")) { + energy = atoi(gSystem->Getenv("CONFIG_ENERGY")); + cout<<"Energy set to "<Getenv("CONFIG_SEED")) { + seed = atoi(gSystem->Getenv("CONFIG_SEED")); + } +} diff --git a/prod/LHC09d6/CreateAODfromESD.C b/prod/LHC09d6/CreateAODfromESD.C new file mode 100644 index 00000000000..9c4148736af --- /dev/null +++ b/prod/LHC09d6/CreateAODfromESD.C @@ -0,0 +1,137 @@ +#if !defined(__CINT__) || defined(__MAKECINT__) +#include +#include +#include "AliAnalysisManager.h" +#include "AliESDInputHandler.h" +#include "AliAODHandler.h" +#include "AliAnalysisTaskESDfilter.h" +#include "AliAnalysisDataContainer.h" +#endif + +void CreateAODfromESD(const char *inFileName = "AliESDs.root", + const char *outFileName = "AliAODs.root", + Bool_t bKineFilter = kTRUE) +{ + + gSystem->Load("libTree"); + gSystem->Load("libGeom"); + gSystem->Load("libPhysics"); + gSystem->Load("libVMC"); + gSystem->Load("libSTEERBase"); + gSystem->Load("libESD"); + gSystem->Load("libAOD"); + + gSystem->Load("libANALYSIS"); + gSystem->Load("libANALYSISalice"); + gSystem->Load("libCORRFW"); + gSystem->Load("libPWG3muon"); + + TChain *chain = new TChain("esdTree"); + // Steering input chain + chain->Add(inFileName); + AliAnalysisManager *mgr = new AliAnalysisManager("ESD to AOD", "Analysis Manager"); + + // Input + AliESDInputHandler* inpHandler = new AliESDInputHandler(); + inpHandler->SetReadTags(); + mgr->SetInputEventHandler (inpHandler); + // Output + AliAODHandler* aodHandler = new AliAODHandler(); + aodHandler->SetOutputFileName(outFileName); + mgr->SetOutputEventHandler(aodHandler); + + // MC Truth + if(bKineFilter){ + AliMCEventHandler* mcHandler = new AliMCEventHandler(); + mgr->SetMCtruthEventHandler(mcHandler); + } + + + // Tasks + // 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); + + // Barrel Tracks + AliAnalysisTaskESDfilter *filter = new AliAnalysisTaskESDfilter("Filter"); + mgr->AddTask(filter); + // Muons + AliAnalysisTaskESDMuonFilter *esdmuonfilter = new AliAnalysisTaskESDMuonFilter("ESD Muon Filter"); + mgr->AddTask(esdmuonfilter); + + // Cuts on primary tracks + AliESDtrackCuts* esdTrackCutsL = new AliESDtrackCuts("AliESDtrackCuts", "Standard"); + esdTrackCutsL->SetMinNClustersTPC(50); + esdTrackCutsL->SetMaxChi2PerClusterTPC(3.5); + esdTrackCutsL->SetMaxCovDiagonalElements(2, 2, 0.5, 0.5, 2); + esdTrackCutsL->SetRequireTPCRefit(kTRUE); + esdTrackCutsL->SetMaxDCAToVertexXY(3.0); + esdTrackCutsL->SetMaxDCAToVertexZ(3.0); + esdTrackCutsL->SetDCAToVertex2D(kTRUE); + esdTrackCutsL->SetRequireSigmaToVertex(kFALSE); + esdTrackCutsL->SetAcceptKinkDaughters(kFALSE); + // ITS stand-alone tracks + AliESDtrackCuts* esdTrackCutsITSsa = new AliESDtrackCuts("AliESDtrackCuts", "ITS stand-alone"); + esdTrackCutsITSsa->SetRequireITSStandAlone(kTRUE); + + AliAnalysisFilter* trackFilter = new AliAnalysisFilter("trackFilter"); + trackFilter->AddCuts(esdTrackCutsL); + trackFilter->AddCuts(esdTrackCutsITSsa); + + // Cuts on V0s + AliESDv0Cuts* esdV0Cuts = new AliESDv0Cuts("AliESDv0Cuts", "Standard pp"); + esdV0Cuts->SetMinRadius(0.2); + esdV0Cuts->SetMaxRadius(200); + esdV0Cuts->SetMinDcaPosToVertex(0.05); + esdV0Cuts->SetMinDcaNegToVertex(0.05); + esdV0Cuts->SetMaxDcaV0Daughters(1.0); + esdV0Cuts->SetMinCosinePointingAngle(0.99); + AliAnalysisFilter* v0Filter = new AliAnalysisFilter("v0Filter"); + v0Filter->AddCuts(esdV0Cuts); + + +// + filter->SetTrackFilter(trackFilter); + filter->SetV0Filter(v0Filter); + + +// Create AOD Tags + AliAnalysisTaskTagCreator* tagTask = new AliAnalysisTaskTagCreator("AOD Tag Creator"); + mgr->AddTask(tagTask); + + // Pipelining + AliAnalysisDataContainer *cinput1 = mgr->GetCommonInputContainer(); + AliAnalysisDataContainer *coutput1 = mgr->GetCommonOutputContainer(); + + + AliAnalysisDataContainer *coutputT + = mgr->CreateContainer("cTag", TTree::Class(), AliAnalysisManager::kOutputContainer, "AOD.tag.root"); + + coutput1->SetSpecialOutput(); + coutputT->SetSpecialOutput(); + + if(bKineFilter) { + mgr->ConnectInput (kinefilter, 0, cinput1 ); + mgr->ConnectOutput (kinefilter, 0, coutput1 ); + } + + mgr->ConnectInput (filter, 0, cinput1 ); + mgr->ConnectOutput(filter, 0, coutput1); + + mgr->ConnectInput (esdmuonfilter, 0, cinput1 ); +// mgr->ConnectOutput(esdmuonfilter, 0, coutput1); + + mgr->ConnectInput (tagTask, 0, cinput1); + mgr->ConnectOutput(tagTask, 1, coutputT); + + // + // Run the analysis + // + mgr->InitAnalysis(); + mgr->PrintStatus(); + mgr->StartAnalysis("local", chain); +} + diff --git a/prod/LHC09d6/JDL b/prod/LHC09d6/JDL new file mode 100644 index 00000000000..262cb1ed0fa --- /dev/null +++ b/prod/LHC09d6/JDL @@ -0,0 +1,31 @@ +Executable = "aliroot_new"; +Jobtag={"comment:Early physics (2009 - FIRST DATA) pp, Phojet, 0T, 900GeV, ID #133"}; + +Packages={"VO_ALICE@AliRoot::v4-17-Rev-16","VO_ALICE@GEANT3::v1-11-4","VO_ALICE@ROOT::v20091109-1","VO_ALICE@APISCONFIG::V1.0x"}; + +TTL = "72000"; + +Requirements = member(other.GridPartitions,"fastmc"); + +Validationcommand ="/alice/cern.ch/user/a/aliprod/prod2007/configs_pbpb_hijing/validation.sh"; + +InputFile= {"LF:/alice/cern.ch/user/a/aliprod/LHC09d6/CheckESD.C", + "LF:/alice/cern.ch/user/a/aliprod/LHC09d6/Config.C", + "LF:/alice/cern.ch/user/a/aliprod/LHC09d6/rec.C", + "LF:/alice/cern.ch/user/a/aliprod/LHC09d6/sim.C", + "LF:/alice/cern.ch/user/a/aliprod/LHC09d6/simrun.C", + "LF:/alice/cern.ch/user/a/aliprod/LHC09d6/tag.C", + "LF:/alice/cern.ch/user/a/aliprod/LHC09d6/CreateAODfromESD.C"}; + + +OutputArchive={"log_archive:*.log,stdout,stderr@ALICE::JINR::SE","root_archive.zip:galice.root,Kinematics.root,TrackRefs.root,AliESDs.root,AliESDfriends.root,AliAODs.root,*QA*.root,ITS.RecPoints.root,AliITSPlaneEffSPDtracklet.root,AliITSPlaneEffSPDtrackletBkg.root,TrackletsMCpred.root,TrackleterSPDHistos.root,Run*.root@ALICE::CERN::ALICEDISK,ALICE::FZK::SE,ALICE::CNAF::SE"}; + +OutputDir="/alice/sim/LHC09d6/$1/#alien_counter_03i#"; + +JDLVariables={"Packages", "OutputDir"}; +GUIDFILE="guid.txt"; + +splitarguments="simrun.C --run $1 --event #alien_counter# --process kPhojet --field kNoField --energy 900" ; +split="production:1-1000"; + +Workdirectorysize={"10000MB"}; diff --git a/prod/LHC09d6/rec.C b/prod/LHC09d6/rec.C new file mode 100644 index 00000000000..d883ceee088 --- /dev/null +++ b/prod/LHC09d6/rec.C @@ -0,0 +1,53 @@ +void rec() { + + AliReconstruction reco; + + reco.SetWriteESDfriend(); + reco.SetWriteAlignmentData(); + reco.SetRunReconstruction("ITS VZERO"); + + // ITS Efficiency + reco.SetRunPlaneEff(kTRUE); + AliITSRecoParam *itsRecoParam = AliITSRecoParam::GetPlaneEffParam(-1); + itsRecoParam->SetOptTrackletsPlaneEff(kTRUE); + + + //****** FIRST PHYSICS 2009 (same as COSMICS 2009) ********************* + //AliITSRecoParam * itsRecoParam = AliITSRecoParam::GetLowFluxParam(); + itsRecoParam->SetClusterErrorsParam(2); + // find independently ITS SA tracks + itsRecoParam->SetSAUseAllClusters(); + itsRecoParam->SetOuterStartLayerSA(AliITSgeomTGeo::GetNLayers()-2); + // to maximize efficiency + itsRecoParam->SetAllowProlongationWithEmptyRoad(); // larger seach windows for SA (in case of large misalignments) + itsRecoParam->SetNLoopsSA(33); + itsRecoParam->SetFactorSAWindowSizes(20); + // additional error due to misal (B off) + itsRecoParam->SetClusterMisalErrorY(1.0,1.0,1.0,1.0,1.0,1.0); // [cm] + itsRecoParam->SetClusterMisalErrorZ(1.0,1.0,1.0,1.0,1.0,1.0); // [cm] + // additional error due to misal (B on) + itsRecoParam->SetClusterMisalErrorYBOn(0.0,0.0,0.1,0.1,0.1,0.1); // [cm] + itsRecoParam->SetClusterMisalErrorZBOn(0.1,0.1,0.1,0.1,0.1,0.1); // [cm] + //*********************************************************************** + itsRecoParam->SetEventSpecie(AliRecoParam::kLowMult); + + reco.SetRecoParam("ITS",itsRecoParam); + + + + reco.SetDefaultStorage("alien://Folder=/alice/simulation/2008/v4-15-Release/Residual/"); +// reco.SetSpecificStorage("GRP/GRP/Data", +// Form("local://%s",gSystem->pwd())); + // We store the object in AliEn during the simulation + reco.SetSpecificStorage("GRP/GRP/Data", + "alien://Folder=/alice/simulation/2008/v4-15-Release/Ideal/"); + + + reco.SetSpecificStorage("GRP/Calib/MeanVertexSPD","alien://folder=/alice/cern.ch/user/r/rgrosso/ShiftedVertex"); + + TStopwatch timer; + timer.Start(); + reco.Run(); + timer.Stop(); + timer.Print(); +} diff --git a/prod/LHC09d6/sim.C b/prod/LHC09d6/sim.C new file mode 100644 index 00000000000..ebc41521069 --- /dev/null +++ b/prod/LHC09d6/sim.C @@ -0,0 +1,20 @@ +void sim(Int_t nev=200) { + AliSimulation simulator; + simulator.SetRunHLT(""); + simulator.SetMakeSDigits("VZERO"); + simulator.SetMakeDigitsFromHits("ITS"); + simulator.SetMakeDigits("ITS VZERO"); + +// simulator.SetWriteRawData("ALL","raw.root",kTRUE); + + simulator.SetDefaultStorage("alien://Folder=/alice/simulation/2008/v4-15-Release/Ideal/"); + +// simulator.SetSpecificStorage("GRP/GRP/Data", +// Form("local://%s",gSystem->pwd())); + + TStopwatch timer; + timer.Start(); + simulator.Run(nev); + timer.Stop(); + timer.Print(); +} diff --git a/prod/LHC09d6/simrun.C b/prod/LHC09d6/simrun.C new file mode 100644 index 00000000000..068cb6b07ee Binary files /dev/null and b/prod/LHC09d6/simrun.C differ diff --git a/prod/LHC09d6/tag.C b/prod/LHC09d6/tag.C new file mode 100644 index 00000000000..668d95ae835 --- /dev/null +++ b/prod/LHC09d6/tag.C @@ -0,0 +1,108 @@ +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; +}