#include "AlidNdEtaCorrectionTask.h"
+#include <TROOT.h>
#include <TCanvas.h>
#include <TChain.h>
#include <TFile.h>
#include <TH1F.h>
+#include <TH2F.h>
+#include <TProfile.h>
#include <TParticle.h>
#include <AliLog.h>
#include <AliMCEvent.h>
#include <AliESDInputHandler.h>
-#include "esdTrackCuts/AliESDtrackCuts.h"
+#include "AliESDtrackCuts.h"
#include "AliPWG0Helper.h"
//#include "AliCorrection.h"
//#include "AliCorrectionMatrix3D.h"
#include "dNdEta/dNdEtaAnalysis.h"
#include "dNdEta/AlidNdEtaCorrection.h"
+#include "AliVertexerTracks.h"
ClassImp(AlidNdEtaCorrectionTask)
fOutput(0),
fOption(opt),
fAnalysisMode(AliPWG0Helper::kTPC),
+ fTrigger(AliPWG0Helper::kMB1),
fSignMode(0),
+ fOnlyPrimaries(kFALSE),
fEsdTrackCuts(0),
fdNdEtaCorrection(0),
fdNdEtaAnalysisMC(0),
fdNdEtaAnalysisESD(0),
fPIDParticles(0),
fPIDTracks(0),
+ fVertexCorrelation(0),
+ fVertexProfile(0),
+ fVertexShiftNorm(0),
+ fEtaCorrelation(0),
+ fEtaProfile(0),
fSigmaVertexTracks(0),
- fSigmaVertexPrim(0)
+ fSigmaVertexPrim(0),
+ fMultAll(0),
+ fMultTr(0),
+ fMultVtx(0),
+ fEventStats(0)
{
//
// Constructor. Initialization of pointers
for (Int_t i=0; i<2; i++)
fdNdEtaCorrectionProcessType[i] = 0;
+
+ for (Int_t i=0; i<8; i++)
+ fDeltaPhi[i] = 0;
}
AlidNdEtaCorrectionTask::~AlidNdEtaCorrectionTask()
fOutput->Add(fdNdEtaAnalysisESD);
if (fOption.Contains("process-types")) {
- fdNdEtaCorrectionProcessType[0] = new AlidNdEtaCorrection("dndeta_correction_ND", "dndeta_correction_ND");
- fdNdEtaCorrectionProcessType[1] = new AlidNdEtaCorrection("dndeta_correction_SD", "dndeta_correction_SD");
- fdNdEtaCorrectionProcessType[2] = new AlidNdEtaCorrection("dndeta_correction_DD", "dndeta_correction_DD");
+ fdNdEtaCorrectionProcessType[0] = new AlidNdEtaCorrection("dndeta_correction_ND", "dndeta_correction_ND", fAnalysisMode);
+ fdNdEtaCorrectionProcessType[1] = new AlidNdEtaCorrection("dndeta_correction_SD", "dndeta_correction_SD", fAnalysisMode);
+ fdNdEtaCorrectionProcessType[2] = new AlidNdEtaCorrection("dndeta_correction_DD", "dndeta_correction_DD", fAnalysisMode);
fOutput->Add(fdNdEtaCorrectionProcessType[0]);
fOutput->Add(fdNdEtaCorrectionProcessType[1]);
fOutput->Add(fSigmaVertexPrim);
Printf("WARNING: sigma-vertex analysis enabled. This will produce weird results in the AliESDtrackCuts histograms");
}
+
+ fVertexCorrelation = new TH2F("fVertexCorrelation", "fVertexCorrelation;MC z-vtx;ESD z-vtx", 120, -30, 30, 120, -30, 30);
+ fOutput->Add(fVertexCorrelation);
+ fVertexProfile = new TProfile("fVertexProfile", "fVertexProfile;MC z-vtx;MC z-vtx - ESD z-vtx", 40, -20, 20);
+ fOutput->Add(fVertexProfile);
+ fVertexShiftNorm = new TH1F("fVertexShiftNorm", "fVertexShiftNorm;(MC z-vtx - ESD z-vtx) / #sigma_{ESD z-vtx};Entries", 200, -100, 100);
+ fOutput->Add(fVertexShiftNorm);
+
+ fEtaCorrelation = new TH2F("fEtaCorrelation", "fEtaCorrelation;MC #eta;ESD #eta", 120, -3, 3, 120, -3, 3);
+ fOutput->Add(fEtaCorrelation);
+ fEtaProfile = new TProfile("fEtaProfile", "fEtaProfile;MC #eta;MC #eta - ESD #eta", 120, -3, 3);
+ fOutput->Add(fEtaProfile);
+
+ fMultAll = new TH1F("fMultAll", "fMultAll", 500, -0.5, 499.5);
+ fOutput->Add(fMultAll);
+ fMultVtx = new TH1F("fMultVtx", "fMultVtx", 500, -0.5, 499.5);
+ fOutput->Add(fMultVtx);
+ fMultTr = new TH1F("fMultTr", "fMultTr", 500, -0.5, 499.5);
+ fOutput->Add(fMultTr);
+
+ for (Int_t i=0; i<8; i++)
+ {
+ fDeltaPhi[i] = new TH1F(Form("fDeltaPhi_%d", i), ";#Delta phi;Entries", 2000, -0.1, 0.1);
+ fOutput->Add(fDeltaPhi[i]);
+ }
+
+ fEventStats = new TH2F("fEventStats", "fEventStats;event type;status;count", 2, -0.5, 1.5, 4, -0.5, 3.5);
+ fOutput->Add(fEventStats);
+ fEventStats->GetXaxis()->SetBinLabel(1, "INEL");
+ fEventStats->GetXaxis()->SetBinLabel(2, "NSD");
+
+ fEventStats->GetYaxis()->SetBinLabel(1, "nothing");
+ fEventStats->GetYaxis()->SetBinLabel(2, "trg");
+ fEventStats->GetYaxis()->SetBinLabel(3, "vtx");
+ fEventStats->GetYaxis()->SetBinLabel(4, "trgvtx");
}
void AlidNdEtaCorrectionTask::Exec(Option_t*)
return;
}
+ if (fOnlyPrimaries)
+ Printf("WARNING: Processing only primaries. For systematical studies only!");
+
// trigger definition
- Bool_t eventTriggered = AliPWG0Helper::IsEventTriggered(fESD->GetTriggerMask(), AliPWG0Helper::kMB1);
+ Bool_t eventTriggered = AliPWG0Helper::IsEventTriggered(fESD->GetTriggerMask(), fTrigger);
- Bool_t eventVertex = kFALSE;
- if (AliPWG0Helper::GetVertex(fESD, fAnalysisMode))
- eventVertex = kTRUE;
+ if (!eventTriggered)
+ Printf("No trigger");
// post the data already here
PostData(0, fOutput);
+ // MC info
+ AliMCEventHandler* eventHandler = dynamic_cast<AliMCEventHandler*> (AliAnalysisManager::GetAnalysisManager()->GetMCtruthEventHandler());
+ if (!eventHandler) {
+ Printf("ERROR: Could not retrieve MC event handler");
+ return;
+ }
+
+ AliMCEvent* mcEvent = eventHandler->MCEvent();
+ if (!mcEvent) {
+ Printf("ERROR: Could not retrieve MC event");
+ return;
+ }
+
+ AliStack* stack = mcEvent->Stack();
+ if (!stack)
+ {
+ AliDebug(AliLog::kError, "Stack not available");
+ return;
+ }
+
+ AliHeader* header = mcEvent->Header();
+ if (!header)
+ {
+ AliDebug(AliLog::kError, "Header not available");
+ return;
+ }
+
+ // get process type; NB: this only works for Pythia
+ Int_t processType = AliPWG0Helper::GetPythiaEventProcessType(header);
+ AliDebug(AliLog::kDebug+1, Form("Found pythia process type %d", processType));
+
+ if (processType<0)
+ AliDebug(AliLog::kError, Form("Unknown Pythia process type %d.", processType));
+
+ // get the MC vertex
+ AliGenEventHeader* genHeader = header->GenEventHeader();
+ TArrayF vtxMC(3);
+ genHeader->PrimaryVertex(vtxMC);
+
+ // get the ESD vertex
+ const AliESDVertex* vtxESD = AliPWG0Helper::GetVertex(fESD, fAnalysisMode);
+
+ Bool_t eventVertex = kFALSE;
+ Double_t vtx[3];
+ if (vtxESD)
+ {
+ vtxESD->GetXYZ(vtx);
+ eventVertex = kTRUE;
+
+ Double_t diff = vtxMC[2] - vtx[2];
+ if (vtxESD->GetZRes() > 0)
+ fVertexShiftNorm->Fill(diff / vtxESD->GetZRes());
+ }
+ else
+ Printf("No vertex found");
+
+ // fill process type
+ Int_t biny = (Int_t) eventTriggered + 2 * (Int_t) eventVertex;
+ // INEL
+ fEventStats->Fill(0.0, biny);
+ // NDS
+ if (processType != 92 && processType != 93)
+ fEventStats->Fill(1.0, biny);
+
// create list of (label, eta, pt) tuples
Int_t inputCount = 0;
Int_t* labelArr = 0;
+ Int_t* labelArr2 = 0; // only for case of SPD
Float_t* etaArr = 0;
Float_t* ptArr = 0;
+ Float_t* deltaPhiArr = 0;
if (fAnalysisMode == AliPWG0Helper::kSPD)
{
// get tracklets
}
labelArr = new Int_t[mult->GetNumberOfTracklets()];
+ labelArr2 = new Int_t[mult->GetNumberOfTracklets()];
etaArr = new Float_t[mult->GetNumberOfTracklets()];
ptArr = new Float_t[mult->GetNumberOfTracklets()];
+ deltaPhiArr = new Float_t[mult->GetNumberOfTracklets()];
// get multiplicity from ITS tracklets
for (Int_t i=0; i<mult->GetNumberOfTracklets(); ++i)
{
//printf("%d %f %f %f\n", i, mult->GetTheta(i), mult->GetPhi(i), mult->GetDeltaPhi(i));
- // This removes non-tracklets in PDC06 data. Very bad solution. New solution is implemented for newer data. Keeping this for compatibility.
- if (mult->GetDeltaPhi(i) < -1000)
- continue;
+ Float_t deltaPhi = mult->GetDeltaPhi(i);
+ // prevent values to be shifted by 2 Pi()
+ if (deltaPhi < -TMath::Pi())
+ deltaPhi += TMath::Pi() * 2;
+ if (deltaPhi > TMath::Pi())
+ deltaPhi -= TMath::Pi() * 2;
+
+ if (TMath::Abs(deltaPhi) > 1)
+ printf("WARNING: Very high Delta Phi: %d %f %f %f\n", i, mult->GetTheta(i), mult->GetPhi(i), deltaPhi);
+
+ if (fOnlyPrimaries)
+ if (mult->GetLabel(i, 0) < 0 || mult->GetLabel(i, 0) != mult->GetLabel(i, 1) || !stack->IsPhysicalPrimary(mult->GetLabel(i, 0)))
+ continue;
etaArr[inputCount] = mult->GetEta(i);
- labelArr[inputCount] = mult->GetLabel(i);
+ labelArr[inputCount] = mult->GetLabel(i, 0);
+ labelArr2[inputCount] = mult->GetLabel(i, 1);
ptArr[inputCount] = 0; // no pt for tracklets
+ deltaPhiArr[inputCount] = deltaPhi;
++inputCount;
}
}
TObjArray* list = fEsdTrackCuts->GetAcceptedTracks(fESD);
Int_t nGoodTracks = list->GetEntries();
+ Printf("Accepted %d tracks", nGoodTracks);
+
labelArr = new Int_t[nGoodTracks];
+ labelArr2 = new Int_t[nGoodTracks];
etaArr = new Float_t[nGoodTracks];
ptArr = new Float_t[nGoodTracks];
+ deltaPhiArr = new Float_t[nGoodTracks];
// loop over esd tracks
for (Int_t i=0; i<nGoodTracks; i++)
etaArr[inputCount] = esdTrack->Eta();
labelArr[inputCount] = TMath::Abs(esdTrack->GetLabel());
+ labelArr2[inputCount] = labelArr[inputCount]; // no second label for tracks
ptArr[inputCount] = esdTrack->Pt();
+ deltaPhiArr[inputCount] = 0; // no delta phi for tracks
++inputCount;
}
- }
- else
- return;
- AliMCEventHandler* eventHandler = dynamic_cast<AliMCEventHandler*> (AliAnalysisManager::GetAnalysisManager()->GetMCtruthEventHandler());
- if (!eventHandler) {
- Printf("ERROR: Could not retrieve MC event handler");
- return;
- }
-
- AliMCEvent* mcEvent = eventHandler->MCEvent();
- if (!mcEvent) {
- Printf("ERROR: Could not retrieve MC event");
- return;
+ delete list;
}
-
- AliStack* stack = mcEvent->Stack();
- if (!stack)
- {
- AliDebug(AliLog::kError, "Stack not available");
+ else
return;
- }
- AliHeader* header = mcEvent->Header();
- if (!header)
+ if (eventTriggered && eventVertex)
{
- AliDebug(AliLog::kError, "Header not available");
- return;
+ fVertexCorrelation->Fill(vtxMC[2], vtx[2]);
+ fVertexProfile->Fill(vtxMC[2], vtxMC[2] - vtx[2]);
}
- // get process type; NB: this only works for Pythia
- Int_t processType = AliPWG0Helper::GetPythiaEventProcessType(header);
- AliDebug(AliLog::kDebug+1, Form("Found pythia process type %d", processType));
-
- if (processType<0)
- AliDebug(AliLog::kError, Form("Unknown Pythia process type %d.", processType));
-
- // get the MC vertex
- AliGenEventHeader* genHeader = header->GenEventHeader();
- TArrayF vtxMC(3);
- genHeader->PrimaryVertex(vtxMC);
// loop over mc particles
Int_t nPrim = stack->GetNprimary();
+ Int_t nAccepted = 0;
for (Int_t iMc = 0; iMc < nPrim; ++iMc)
{
if (eventTriggered)
if (eventVertex)
fdNdEtaAnalysisMC->FillTrack(vtxMC[2], eta, pt);
+
+ // TODO this value might be needed lower for the SPD study (only used for control histograms anyway)
+ if (TMath::Abs(eta) < 1 && pt > 0.2)
+ nAccepted++;
}
+ fMultAll->Fill(nAccepted);
+ if (eventTriggered) {
+ fMultTr->Fill(nAccepted);
+ if (eventVertex)
+ fMultVtx->Fill(nAccepted);
+ }
+
+ Int_t processed = 0;
+
for (Int_t i=0; i<inputCount; ++i)
{
Int_t label = labelArr[i];
+ Int_t label2 = labelArr2[i];
if (label < 0)
{
- AliDebug(AliLog::kWarning, Form("WARNING: cannot find corresponding mc part for track %d.", label));
+ Printf("WARNING: cannot find corresponding mc particle for track(let) %d with label %d.", i, label);
continue;
}
AliDebug(AliLog::kError, Form("UNEXPECTED: particle with label %d not found in stack (track loop).", label));
continue;
}
-
- if (SignOK(particle->GetPDG()) == kFALSE)
- continue;
+
+ // find particle that is filled in the correction map
+ // this should be the particle that has been reconstructed
+ // for TPC this is clear and always identified by the track's label
+ // for SPD the labels might not be identical. In this case the particle closest to the beam line that is a primary is taken:
+ // L1 L2 (P = primary, S = secondary)
+ // P P' : --> P
+ // P S : --> P
+ // S P : --> P
+ // S S' : --> S
+
+ if (label != label2 && label2 >= 0)
+ {
+ if (stack->IsPhysicalPrimary(label) == kFALSE && stack->IsPhysicalPrimary(label2))
+ {
+ particle = stack->Particle(label2);
+ if (!particle)
+ {
+ AliDebug(AliLog::kError, Form("UNEXPECTED: particle with label %d not found in stack (track loop).", label2));
+ continue;
+ }
+ }
+ }
if (eventTriggered && eventVertex)
{
- fdNdEtaCorrection->FillTrackedParticle(vtxMC[2], particle->Eta(), particle->Pt());
+ if (SignOK(particle->GetPDG()) == kFALSE)
+ continue;
+
+ processed++;
+
+ Bool_t firstIsPrim = stack->IsPhysicalPrimary(label);
+ // in case of primary the MC values are filled, otherwise (background) the reconstructed values
+ if (label == label2 && firstIsPrim)
+ {
+ fdNdEtaCorrection->FillTrackedParticle(vtxMC[2], particle->Eta(), particle->Pt());
+ }
+ else
+ {
+ fdNdEtaCorrection->FillTrackedParticle(vtxMC[2], etaArr[i], ptArr[i]);
+ }
+
+ fEtaProfile->Fill(particle->Eta(), particle->Eta() - etaArr[i]);
fdNdEtaAnalysisESD->FillTrack(vtxMC[2], particle->Eta(), particle->Pt());
+
+ fEtaCorrelation->Fill(etaArr[i], particle->Eta());
+
if (particle->Pt() > 0.1 && particle->Pt() < 0.2)
{
fPIDTracks->Fill(particle->GetPdgCode());
if (processType==94)
fdNdEtaCorrectionProcessType[2]->FillTrackedParticle(vtxMC[2], particle->Eta(), particle->Pt());
}
+
+ // control histograms
+ Int_t hist = -1;
+ if (label == label2)
+ {
+ if (firstIsPrim)
+ {
+ hist = 0;
+ }
+ else
+ hist = 1;
+ }
+ else if (label2 >= 0)
+ {
+ Bool_t secondIsPrim = stack->IsPhysicalPrimary(label2);
+ if (firstIsPrim && secondIsPrim)
+ {
+ hist = 2;
+ }
+ else if (firstIsPrim && !secondIsPrim)
+ {
+ hist = 3;
+
+ // check if secondary is caused by the primary or it is a fake combination
+ //Printf("PS case --> Label 1 is %d, label 2 is %d", label, label2);
+ Int_t mother = label2;
+ while (stack->Particle(mother) && stack->Particle(mother)->GetMother(0) >= 0)
+ {
+ //Printf(" %d created by %d, %d", mother, stack->Particle(mother)->GetMother(0), stack->Particle(mother)->GetMother(1));
+ if (stack->Particle(mother)->GetMother(0) == label)
+ {
+ /*Printf("The untraceable decay was:");
+ Printf(" from:");
+ particle->Print();
+ Printf(" to (status code %d):", stack->Particle(mother)->GetStatusCode());
+ stack->Particle(mother)->Print();*/
+ hist = 4;
+ }
+ mother = stack->Particle(mother)->GetMother(0);
+ }
+ }
+ else if (!firstIsPrim && secondIsPrim)
+ {
+ hist = 5;
+ }
+ else if (!firstIsPrim && !secondIsPrim)
+ {
+ hist = 6;
+ }
+
+ }
+ else
+ hist = 7;
+
+ fDeltaPhi[hist]->Fill(deltaPhiArr[i]);
}
}
+ if (processed < inputCount)
+ Printf("Only %d out of %d track(let)s were used", processed, inputCount);
+
if (eventTriggered && eventVertex)
{
fdNdEtaAnalysisMC->FillEvent(vtxMC[2], inputCount);
delete[] etaArr;
delete[] labelArr;
+ delete[] labelArr2;
delete[] ptArr;
+ delete[] deltaPhiArr;
// fills the fSigmaVertex histogram (systematic study)
if (fSigmaVertexTracks)
if (fSigmaVertexPrim)
fSigmaVertexPrim->Write();
+ fVertexCorrelation = dynamic_cast<TH2F*> (fOutput->FindObject("fVertexCorrelation"));
+ if (fVertexCorrelation)
+ fVertexCorrelation->Write();
+ fVertexProfile = dynamic_cast<TProfile*> (fOutput->FindObject("fVertexProfile"));
+ if (fVertexProfile)
+ fVertexProfile->Write();
+ fVertexShiftNorm = dynamic_cast<TH1F*> (fOutput->FindObject("fVertexShiftNorm"));
+ if (fVertexShiftNorm)
+ fVertexShiftNorm->Write();
+
+ fEtaCorrelation = dynamic_cast<TH2F*> (fOutput->FindObject("fEtaCorrelation"));
+ if (fEtaCorrelation)
+ fEtaCorrelation->Write();
+ fEtaProfile = dynamic_cast<TProfile*> (fOutput->FindObject("fEtaProfile"));
+ if (fEtaProfile)
+ fEtaProfile->Write();
+
+ fMultAll = dynamic_cast<TH1F*> (fOutput->FindObject("fMultAll"));
+ if (fMultAll)
+ fMultAll->Write();
+
+ fMultTr = dynamic_cast<TH1F*> (fOutput->FindObject("fMultTr"));
+ if (fMultTr)
+ fMultTr->Write();
+
+ fMultVtx = dynamic_cast<TH1F*> (fOutput->FindObject("fMultVtx"));
+ if (fMultVtx)
+ fMultVtx->Write();
+
+ for (Int_t i=0; i<8; ++i)
+ {
+ fDeltaPhi[i] = dynamic_cast<TH1F*> (fOutput->FindObject(Form("fDeltaPhi_%d", i)));
+ if (fDeltaPhi[i])
+ fDeltaPhi[i]->Write();
+ }
+
+ fEventStats = dynamic_cast<TH2F*> (fOutput->FindObject("fEventStats"));
+ if (fEventStats)
+ fEventStats->Write();
+
fout->Write();
fout->Close();
- fdNdEtaCorrection->DrawHistograms();
+ //fdNdEtaCorrection->DrawHistograms();
Printf("Writting result to %s", fileName.Data());