#include <TParticle.h>
#include <TParticlePDG.h>
#include <TDatabasePDG.h>
+#include <TRandom.h>
#include <AliLog.h>
#include <AliESDVertex.h>
#include <AliESDEvent.h>
+#include <AliESDRun.h>
#include <AliStack.h>
#include <AliHeader.h>
#include <AliGenEventHeader.h>
#include "AliPWG0Helper.h"
#include "dNdEta/dNdEtaAnalysis.h"
#include "dNdEta/AlidNdEtaCorrection.h"
+#include "AliTriggerAnalysis.h"
+#include "AliPhysicsSelection.h"
ClassImp(AlidNdEtaCorrectionTask)
fESD(0),
fOutput(0),
fOption(),
- fAnalysisMode(AliPWG0Helper::kTPC),
- fTrigger(AliPWG0Helper::kMB1),
+ fAnalysisMode((AliPWG0Helper::AnalysisMode) (AliPWG0Helper::kTPC | AliPWG0Helper::kFieldOn)),
+ fTrigger(AliTriggerAnalysis::kMB1),
fFillPhi(kFALSE),
fDeltaPhiCut(-1),
+ fSymmetrize(kFALSE),
+ fMultAxisEta1(kFALSE),
+ fDiffTreatment(AliPWG0Helper::kMCFlags),
fSignMode(0),
fOnlyPrimaries(kFALSE),
fStatError(0),
+ fSystSkipParticles(kFALSE),
fEsdTrackCuts(0),
fdNdEtaCorrection(0),
fdNdEtaAnalysisMC(0),
fMultAll(0),
fMultTr(0),
fMultVtx(0),
- fEventStats(0)
+ fEventStats(0),
+ fEsdTrackCutsCheck(0),
+ fEtaCorrelationAllESD(0),
+ fpTCorrelation(0),
+ fpTCorrelationShift(0),
+ fpTCorrelationAllESD(0),
+ fpTCorrelationShiftAllESD(0),
+ fPtMin(0.15),
+ fPtMC(0),
+ fEtaMC(0),
+ fPtESD(0),
+ fEtaESD(0),
+ fVtxMC(0),
+ fNumberEventMC(0),
+ fNumberEvent(0),
+ fEventNumber(-1),
+ fWeightSecondaries(kFALSE)
{
//
// Constructor. Initialization of pointers
for (Int_t i=0; i<8; i++)
fDeltaPhi[i] = 0;
+
+ AliLog::SetClassDebugLevel("AlidNdEtaCorrectionTask", AliLog::kWarning);
}
AlidNdEtaCorrectionTask::AlidNdEtaCorrectionTask(const char* opt) :
fESD(0),
fOutput(0),
fOption(opt),
- fAnalysisMode(AliPWG0Helper::kTPC),
- fTrigger(AliPWG0Helper::kMB1),
+ fAnalysisMode((AliPWG0Helper::AnalysisMode) (AliPWG0Helper::kTPC | AliPWG0Helper::kFieldOn)),
+ fTrigger(AliTriggerAnalysis::kMB1),
fFillPhi(kFALSE),
fDeltaPhiCut(0),
+ fSymmetrize(kFALSE),
+ fMultAxisEta1(kFALSE),
+ fDiffTreatment(AliPWG0Helper::kMCFlags),
fSignMode(0),
fOnlyPrimaries(kFALSE),
fStatError(0),
+ fSystSkipParticles(kFALSE),
fEsdTrackCuts(0),
fdNdEtaCorrection(0),
fdNdEtaAnalysisMC(0),
fMultAll(0),
fMultTr(0),
fMultVtx(0),
- fEventStats(0)
+ fEventStats(0),
+ fEsdTrackCutsCheck(0),
+ fEtaCorrelationAllESD(0),
+ fpTCorrelation(0),
+ fpTCorrelationShift(0),
+ fpTCorrelationAllESD(0),
+ fpTCorrelationShiftAllESD(0),
+ fPtMin(0.15),
+ fPtMC(0),
+ fEtaMC(0),
+ fPtESD(0),
+ fEtaESD(0),
+ fVtxMC(0),
+ fNumberEventMC(0),
+ fNumberEvent(0),
+ fEventNumber(-1),
+ fWeightSecondaries(kFALSE)
{
//
// Constructor. Initialization of pointers
for (Int_t i=0; i<8; i++)
fDeltaPhi[i] = 0;
+
+ AliLog::SetClassDebugLevel("AlidNdEtaCorrectionTask", AliLog::kWarning);
}
AlidNdEtaCorrectionTask::~AlidNdEtaCorrectionTask()
// histograms are in the output list and deleted when the output
// list is deleted by the TSelector dtor
-
- if (fOutput) {
+
+ if (fOutput && !AliAnalysisManager::GetAnalysisManager()->IsProofMode()) {
delete fOutput;
fOutput = 0;
}
+
}
//________________________________________________________________________
// Enable only the needed branches
esdH->SetActiveBranches("AliESDHeader Vertex");
- if (fAnalysisMode == AliPWG0Helper::kSPD)
+ if (fAnalysisMode & AliPWG0Helper::kSPD)
esdH->SetActiveBranches("AliESDHeader Vertex AliMultiplicity");
- if (fAnalysisMode == AliPWG0Helper::kTPC || fAnalysisMode == AliPWG0Helper::kTPCITS) {
+ if (fAnalysisMode & AliPWG0Helper::kTPC || fAnalysisMode & AliPWG0Helper::kTPCITS) {
esdH->SetActiveBranches("AliESDHeader Vertex Tracks");
}
}
{
// create result objects and add to output list
+ AliDebug(2,Form("*********************************** fOption = %s", fOption.Data()));
if (fOption.Contains("only-positive"))
{
Printf("INFO: Processing only positive particles.");
{
fEsdTrackCutsPrim = dynamic_cast<AliESDtrackCuts*> (fEsdTrackCuts->Clone("fEsdTrackCutsPrim"));
fEsdTrackCutsSec = dynamic_cast<AliESDtrackCuts*> (fEsdTrackCuts->Clone("fEsdTrackCutsSec"));
+ fEsdTrackCutsCheck = dynamic_cast<AliESDtrackCuts*> (fEsdTrackCuts->Clone("fEsdTrackCutsCheck"));
+ fEsdTrackCutsCheck->SetPtRange(0.15);
+ fEsdTrackCutsCheck->SetEtaRange(-1.,1.);
fOutput->Add(fEsdTrackCutsPrim);
fOutput->Add(fEsdTrackCutsSec);
}
}
- /*
- fTemp1 = new TH2F("fTemp1", "fTemp1", 200, -0.08, 0.08, 200, -0.08, 0.08);
+ //fTemp1 = new TH2F("fTemp1", "fTemp1", 4, 0.5, 4.5, 101, -1.5, 99.5); // nsd study
+ fTemp1 = new TH2F("fTemp1", "fTemp1", 300, -15, 15, 80, -2.0, 2.0);
fOutput->Add(fTemp1);
- fTemp2 = new TH1F("fTemp2", "fTemp2", 2000, -5, 5);
+
+ fTemp2 = new TH2F("fTemp2", "fTemp2", 300, -15, 15, 80, -2.0, 2.0);
fOutput->Add(fTemp2);
- */
fVertexCorrelation = new TH2F("fVertexCorrelation", "fVertexCorrelation;MC z-vtx;ESD z-vtx", 120, -30, 30, 120, -30, 30);
fOutput->Add(fVertexCorrelation);
- fVertexCorrelationShift = new TH2F("fVertexCorrelationShift", "fVertexCorrelationShift;MC z-vtx;MC z-vtx - ESD z-vtx", 120, -30, 30, 100, -1, 1);
+ fVertexCorrelationShift = new TH3F("fVertexCorrelationShift", "fVertexCorrelationShift;MC z-vtx;MC z-vtx - ESD z-vtx;rec. tracks", 120, -30, 30, 100, -1, 1, 100, -0.5, 99.5);
fOutput->Add(fVertexCorrelationShift);
fVertexProfile = new TProfile("fVertexProfile", "fVertexProfile;MC z-vtx;MC z-vtx - ESD z-vtx", 40, -20, 20);
fOutput->Add(fVertexProfile);
fVertexShift = new TH1F("fVertexShift", "fVertexShift;(MC z-vtx - ESD z-vtx);Entries", 201, -2, 2);
fOutput->Add(fVertexShift);
- fVertexShiftNorm = new TH1F("fVertexShiftNorm", "fVertexShiftNorm;(MC z-vtx - ESD z-vtx) / #sigma_{ESD z-vtx};Entries", 200, -100, 100);
+ fVertexShiftNorm = new TH2F("fVertexShiftNorm", "fVertexShiftNorm;(MC z-vtx - ESD z-vtx);rec. tracks;Entries", 200, -100, 100, 100, -0.5, 99.5);
fOutput->Add(fVertexShiftNorm);
fEtaCorrelation = new TH2F("fEtaCorrelation", "fEtaCorrelation;MC #eta;ESD #eta", 120, -3, 3, 120, -3, 3);
fOutput->Add(fEtaCorrelation);
+ fEtaCorrelationAllESD = new TH2F("fEtaCorrelationAllESD", "fEtaCorrelationAllESD;MC #eta;ESD #eta", 120, -3, 3, 120, -3, 3);
+ fOutput->Add(fEtaCorrelationAllESD);
fEtaCorrelationShift = new TH2F("fEtaCorrelationShift", "fEtaCorrelationShift;MC #eta;MC #eta - ESD #eta", 120, -3, 3, 100, -0.1, 0.1);
fOutput->Add(fEtaCorrelationShift);
fEtaProfile = new TProfile("fEtaProfile", "fEtaProfile;MC #eta;MC #eta - ESD #eta", 120, -3, 3);
fpTResolution = new TH2F("fpTResolution", ";MC p_{T} (GeV/c);(MC p_{T} - ESD p_{T}) / MC p_{T}", 160, 0, 20, 201, -0.2, 0.2);
fOutput->Add(fpTResolution);
+ fpTCorrelation = new TH2F("fpTCorrelation", "fpTCorrelation;MC p_{T} (GeV/c);ESD p_{T}", 160, 0, 20, 160, 0, 20);
+ fOutput->Add(fpTCorrelation);
+ fpTCorrelationShift = new TH2F("fpTCorrelationShift", "fpTCorrelationShift;MC p_{T} (GeV/c);MC p_{T} - ESD p_{T}", 160, 0, 20, 100, -1, 1);
+ fOutput->Add(fpTCorrelationShift);
+ fpTCorrelationAllESD = new TH2F("fpTCorrelationAllESD", "fpTCorrelationAllESD;MC p_{T} (GeV/c);ESD p_{T}", 160, 0, 20, 160, 0, 20);
+ fOutput->Add(fpTCorrelationAllESD);
+ fpTCorrelationShiftAllESD = new TH2F("fpTCorrelationShiftAllESD", "fpTCorrelationShiftAllESD;MC p_{T} (GeV/c);MC p_{T} - ESD p_{T}", 160, 0, 20, 100, -1, 1);
+ fOutput->Add(fpTCorrelationShiftAllESD);
+
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(fDeltaPhi[i]);
}
- fEventStats = new TH2F("fEventStats", "fEventStats;event type;status;count", 106, -5.5, 100.5, 4, -0.5, 3.5);
+ fEventStats = new TH2F("fEventStats", "fEventStats;event type;status;count", 109, -6.5, 102.5, 4, -0.5, 3.5);
fOutput->Add(fEventStats);
fEventStats->GetXaxis()->SetBinLabel(1, "INEL"); // x = -5
fEventStats->GetXaxis()->SetBinLabel(2, "NSD"); // x = -4
fEventStats->GetXaxis()->SetBinLabel(3, "ND"); // x = -3
fEventStats->GetXaxis()->SetBinLabel(4, "SD"); // x = -2
fEventStats->GetXaxis()->SetBinLabel(5, "DD"); // x = -1
+
+ fEventStats->GetXaxis()->SetBinLabel(108, "INEL=0"); // x = -101
+ fEventStats->GetXaxis()->SetBinLabel(109, "INEL>0"); // x = -102
- for (Int_t i=0; i<100; i++)
+ for (Int_t i=-1; i<100; i++)
fEventStats->GetXaxis()->SetBinLabel(7+i, Form("%d", i));
fEventStats->GetYaxis()->SetBinLabel(1, "nothing");
if (fEsdTrackCuts)
{
fEsdTrackCuts->SetName("fEsdTrackCuts");
- fOutput->Add(fEsdTrackCuts);
+ // TODO like this we send an empty object back...
+ fOutput->Add(fEsdTrackCuts->Clone());
}
+ fPtMC = new TH1F("fPtMC", "Pt from MC for selected tracks;MC p_{T} (GeV/c)", 160, 0, 20);
+ fOutput->Add(fPtMC);
+ fEtaMC = new TH1F("fEtaMC", "Eta from MC for selected tracks;MC #eta", 120, -3, 3);
+ fOutput->Add(fEtaMC);
+ fPtESD = new TH1F("fPtESD", "Pt from ESD for selected tracks;ESD p_{T} (GeV/c)", 160, 0, 20);
+ fOutput->Add(fPtESD);
+ fEtaESD = new TH1F("fEtaESD", "Eta from ESD for selected tracks;ESD #eta", 120, -3, 3);
+ fOutput->Add(fEtaESD);
+
+ fVtxMC = new TH1F("fVtxMC", "Vtx,z from MC for all events;MC vtx_z (cm)", 100, -30, 30);
+ fOutput->Add(fVtxMC);
+
+ fNumberEventMC = new TH1F("fNumberEventMC","Number of event accepted at MC level",600000,-0.5,600000-0.5);
+ fOutput->Add(fNumberEventMC);
+
+ fNumberEvent = new TH1F("fNumberEvent","Number of event accepted at Reco level",600000,-0.5,600000-0.5);
+ fOutput->Add(fNumberEvent);
}
void AlidNdEtaCorrectionTask::Exec(Option_t*)
{
// process the event
+ fEventNumber++;
// Check prerequisites
if (!fESD)
{
if (fStatError > 0)
Printf("WARNING: Statistical error evaluation active!");
+
+ AliInputEventHandler* inputHandler = dynamic_cast<AliInputEventHandler*> (AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler());
+ if (!inputHandler)
+ {
+ Printf("ERROR: Could not receive input handler");
+ return;
+ }
+
+ Bool_t eventTriggered = inputHandler->IsEventSelected();
- // trigger definition
- Bool_t eventTriggered = AliPWG0Helper::IsEventTriggered(fESD, fTrigger);
-
+ static AliTriggerAnalysis* triggerAnalysis = 0;
+ if (!triggerAnalysis)
+ {
+ AliPhysicsSelection* physicsSelection = dynamic_cast<AliPhysicsSelection*> (inputHandler->GetEventSelection());
+ if (physicsSelection)
+ triggerAnalysis = physicsSelection->GetTriggerAnalysis();
+ }
+
+ if (eventTriggered)
+ eventTriggered = triggerAnalysis->IsTriggerFired(fESD, fTrigger);
+
if (!eventTriggered)
Printf("No trigger");
return;
}
- // get process type;
- Int_t processType = AliPWG0Helper::GetEventProcessType(header);
+ // get process type
+ Int_t processType = AliPWG0Helper::GetEventProcessType(fESD, header, stack, fDiffTreatment);
+
AliDebug(AliLog::kDebug+1, Form("Found process type %d", processType));
if (processType == AliPWG0Helper::kInvalidProcess)
- AliDebug(AliLog::kError, "Unknown process type.");
+ {
+ AliDebug(AliLog::kWarning, "Unknown process type. Setting to ND");
+ processType = AliPWG0Helper::kND;
+ }
// get the MC vertex
AliGenEventHeader* genHeader = header->GenEventHeader();
TArrayF vtxMC(3);
genHeader->PrimaryVertex(vtxMC);
+ fVtxMC->Fill(vtxMC[2]);
+ AliDebug(2,Form("MC vtx: x = %.6f, y = %.6f, z = %.6f",vtxMC[0],vtxMC[1],vtxMC[2]));
// get the ESD vertex
- const AliESDVertex* vtxESD = AliPWG0Helper::GetVertex(fESD, fAnalysisMode);
Bool_t eventVertex = kFALSE;
- if (vtxESD)
+ Double_t vtx[3];
+ const AliESDVertex* vtxESD = AliPWG0Helper::GetVertex(fESD, fAnalysisMode);
+ if (vtxESD && AliPWG0Helper::TestVertex(vtxESD, fAnalysisMode))
{
- Double_t vtx[3];
vtxESD->GetXYZ(vtx);
+ eventVertex = kTRUE;
- Double_t diff = vtxMC[2] - vtx[2];
- fVertexShift->Fill(diff);
- if (vtxESD->GetZRes() > 0)
- fVertexShiftNorm->Fill(diff / vtxESD->GetZRes());
-
- if (!AliPWG0Helper::TestVertex(vtxESD, fAnalysisMode))
+ // remove vertices outside +- 15 cm
+ if (TMath::Abs(vtx[2]) > 15)
{
- vtxESD = 0;
- }
- else
- {
- eventVertex = kTRUE;
-
- if (eventTriggered)
- {
- fVertexCorrelation->Fill(vtxMC[2], vtx[2]);
- fVertexCorrelationShift->Fill(vtxMC[2], vtxMC[2] - vtx[2]);
- fVertexProfile->Fill(vtxMC[2], vtxMC[2] - vtx[2]);
- }
+ eventVertex = kFALSE;
+ vtxESD = 0;
}
}
-
- // fill process type
- Int_t biny = (Int_t) eventTriggered + 2 * (Int_t) eventVertex;
- // INEL
- fEventStats->Fill(-5, biny);
- // NSD
- if (processType != AliPWG0Helper::kSD)
- fEventStats->Fill(-4, biny);
- // SD, ND, DD
- if (processType == AliPWG0Helper::kND)
- fEventStats->Fill(-3, biny);
- if (processType == AliPWG0Helper::kSD)
- fEventStats->Fill(-2, biny);
- if (processType == AliPWG0Helper::kDD)
- fEventStats->Fill(-1, biny);
+ else
+ vtxESD = 0;
// create list of (label, eta, pt) tuples
Int_t inputCount = 0;
Float_t* etaArr = 0;
Float_t* thirdDimArr = 0;
Float_t* deltaPhiArr = 0;
- if (fAnalysisMode == AliPWG0Helper::kSPD)
+ if (fAnalysisMode & AliPWG0Helper::kSPD)
{
- // get tracklets
- const AliMultiplicity* mult = fESD->GetMultiplicity();
- if (!mult)
- {
- AliDebug(AliLog::kError, "AliMultiplicity not available");
- return;
- }
-
- labelArr = new Int_t[mult->GetNumberOfTracklets()];
- labelArr2 = new Int_t[mult->GetNumberOfTracklets()];
- etaArr = new Float_t[mult->GetNumberOfTracklets()];
- thirdDimArr = new Float_t[mult->GetNumberOfTracklets()];
- deltaPhiArr = new Float_t[mult->GetNumberOfTracklets()];
-
- // get multiplicity from SPD tracklets
- for (Int_t i=0; i<mult->GetNumberOfTracklets(); ++i)
+ if (vtxESD)
{
- //printf("%d %f %f %f\n", i, mult->GetTheta(i), mult->GetPhi(i), mult->GetDeltaPhi(i));
-
- Float_t phi = mult->GetPhi(i);
- if (phi < 0)
- phi += TMath::Pi() * 2;
- Float_t deltaPhi = mult->GetDeltaPhi(i);
-
- 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)))
+ // get tracklets
+ const AliMultiplicity* mult = fESD->GetMultiplicity();
+ if (!mult)
+ {
+ AliDebug(AliLog::kError, "AliMultiplicity not available");
+ return;
+ }
+
+ labelArr = new Int_t[mult->GetNumberOfTracklets()];
+ labelArr2 = new Int_t[mult->GetNumberOfTracklets()];
+ etaArr = new Float_t[mult->GetNumberOfTracklets()];
+ thirdDimArr = new Float_t[mult->GetNumberOfTracklets()];
+ deltaPhiArr = new Float_t[mult->GetNumberOfTracklets()];
+
+ Bool_t foundInEta10 = kFALSE;
+
+ // get multiplicity from SPD 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));
+
+ Float_t phi = mult->GetPhi(i);
+ if (phi < 0)
+ phi += TMath::Pi() * 2;
+ Float_t deltaPhi = mult->GetDeltaPhi(i);
+
+ 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;
+
+ if (fDeltaPhiCut > 0 && (TMath::Abs(deltaPhi) > fDeltaPhiCut || TMath::Abs(mult->GetDeltaTheta(i)) > fDeltaPhiCut / 0.08 * 0.025))
continue;
-
- if (fDeltaPhiCut > 0 && TMath::Abs(deltaPhi) > fDeltaPhiCut)
- continue;
+
+ if (fSystSkipParticles && gRandom->Uniform() < 0.0153)
+ {
+ Printf("Skipped tracklet!");
+ continue;
+ }
+
+ // TEST exclude potentially inefficient phi region
+ //if (phi > 5.70 || phi < 0.06)
+ // continue;
+
+ // we have to repeat the trigger here, because the tracklet might have been kicked out fSystSkipParticles
+ if (TMath::Abs(mult->GetEta(i)) < 1)
+ foundInEta10 = kTRUE;
+
+ etaArr[inputCount] = mult->GetEta(i);
+ if (fSymmetrize)
+ etaArr[inputCount] = TMath::Abs(etaArr[inputCount]);
+ labelArr[inputCount] = mult->GetLabel(i, 0);
+ labelArr2[inputCount] = mult->GetLabel(i, 1);
+ thirdDimArr[inputCount] = phi;
+ deltaPhiArr[inputCount] = deltaPhi;
+ ++inputCount;
+ }
- etaArr[inputCount] = mult->GetEta(i);
- labelArr[inputCount] = mult->GetLabel(i, 0);
- labelArr2[inputCount] = mult->GetLabel(i, 1);
- thirdDimArr[inputCount] = phi;
- deltaPhiArr[inputCount] = deltaPhi;
- ++inputCount;
+ /*
+ for (Int_t i=0; i<mult->GetNumberOfSingleClusters(); ++i)
+ {
+ if (TMath::Abs(TMath::Log(TMath::Tan(mult->GetThetaSingle(i)/2.))) < 1);
+ {
+ foundInEta10 = kTRUE;
+ break;
+ }
+ }
+ */
+
+ if (fSystSkipParticles && (fTrigger & AliTriggerAnalysis::kOneParticle) && !foundInEta10)
+ eventTriggered = kFALSE;
}
}
- else if (fAnalysisMode == AliPWG0Helper::kTPC || fAnalysisMode == AliPWG0Helper::kTPCITS)
+ else if (fAnalysisMode & AliPWG0Helper::kTPC || fAnalysisMode & AliPWG0Helper::kTPCITS)
{
if (!fEsdTrackCuts)
{
AliDebug(AliLog::kError, "fESDTrackCuts not available");
return;
}
-
+
+ Bool_t foundInEta10 = kFALSE;
+
if (vtxESD)
{
+ // control histograms on pT
+ Int_t nfromstack = stack->GetNtrack();
+ AliDebug(3,Form(" from stack we have %d tracks\n",nfromstack));
+ for (Int_t itrack = 0; itrack < fESD->GetNumberOfTracks(); itrack++){
+ AliESDtrack* esdTrackcheck = dynamic_cast<AliESDtrack*> (fESD->GetTrack(itrack));
+ if (!esdTrackcheck){
+ AliDebug(AliLog::kError, Form("ERROR: Could not retrieve track %d.", itrack));
+ continue;
+ }
+ if (fOnlyPrimaries){
+ Int_t label = TMath::Abs(esdTrackcheck->GetLabel());
+ AliDebug(4,Form("label = %d\n",label));
+ if (label == 0 || label > nfromstack) continue;
+ if (stack->IsPhysicalPrimary(label) == kFALSE) continue;
+ }
+
+ Int_t label = TMath::Abs(esdTrackcheck->GetLabel());
+ if (label == 0 || label > nfromstack) continue;
+ if (!stack->Particle(label)){
+ AliDebug(4,Form("WARNING: No particle for %d", label));
+ }
+ else{
+ if (!fEsdTrackCuts->AcceptTrack(esdTrackcheck)){
+ TParticle* particle = stack->Particle(label);
+ if (fEsdTrackCutsCheck->AcceptTrack(esdTrackcheck)){
+ //if (TMath::Abs(particle->Eta() < 0.8) && particle->Pt() > fPtMin){
+ Float_t ptMC = particle->Pt();
+ Float_t etaMC = particle->Eta();
+ Float_t ptESD = esdTrackcheck->Pt();
+ Float_t etaESD = esdTrackcheck->Eta();
+ fEtaCorrelationAllESD->Fill(etaMC,etaESD);
+ fpTCorrelationAllESD->Fill(ptMC,ptESD);
+ fpTCorrelationShiftAllESD->Fill(ptMC,ptMC-ptESD);
+ }
+ }
+ }
+ } // end loop over all ESDs
+
// get multiplicity from ESD tracks
- TObjArray* list = fEsdTrackCuts->GetAcceptedTracks(fESD, (fAnalysisMode == AliPWG0Helper::kTPC));
+ TObjArray* list = fEsdTrackCuts->GetAcceptedTracks(fESD, (fAnalysisMode & AliPWG0Helper::kTPC));
Int_t nGoodTracks = list->GetEntries();
Printf("Accepted %d tracks", nGoodTracks);
continue;
}
- // TODO fOnlyPrimaries not implemented for TPC
-
+ AliDebug(3,Form("particle %d: pt = %.6f, eta = %.6f",i,esdTrack->Pt(), esdTrack->Eta()));
+ // 2 Options for INEL>0 trigger - choose one
+ // 1. HL
+ //if (esdTrack->Pt() < 0.15)
+ // foundInEta10 = kTRUE;
+ // 2. MB Working Group definition
+ if (esdTrack->Pt() < fPtMin)
+ continue;
+
+ if (fOnlyPrimaries)
+ {
+ Int_t label = TMath::Abs(esdTrack->GetLabel());
+ if (label == 0)
+ continue;
+
+ if (stack->IsPhysicalPrimary(label) == kFALSE)
+ continue;
+ }
+
+ Int_t label = TMath::Abs(esdTrack->GetLabel());
+ if (!stack->Particle(label)){
+ AliDebug(3,Form("WARNING: No particle for %d", label));
+ }
+ else{
+ TParticle* particle = stack->Particle(label);
+ Float_t ptMC = particle->Pt();
+ Float_t etaMC = particle->Eta();
+ fPtMC->Fill(ptMC);
+ fEtaMC->Fill(etaMC);
+ fPtESD->Fill(esdTrack->Pt());
+ fEtaESD->Fill(esdTrack->Eta());
+ }
+
+ // 2 Options for INEL>0 trigger - choose one
+ // 1. HL
+ //if (TMath::Abs(esdTrack->Eta()) < 1 && esdTrack->Pt() > 0.15)
+ // foundInEta10 = kTRUE;
+ // 2. MB Working Group definition
+ if (TMath::Abs(esdTrack->Eta()) < 0.8 && esdTrack->Pt() > fPtMin)
+ foundInEta10 = kTRUE;
+
etaArr[inputCount] = esdTrack->Eta();
+ if (fSymmetrize)
+ etaArr[inputCount] = TMath::Abs(etaArr[inputCount]);
labelArr[inputCount] = TMath::Abs(esdTrack->GetLabel());
labelArr2[inputCount] = labelArr[inputCount]; // no second label for tracks
thirdDimArr[inputCount] = esdTrack->Pt();
delete list;
- if (eventTriggered)
+ // TODO this code crashes for TPCITS because particles are requested from the stack for some labels that are out of bound
+ if (0 && eventTriggered)
{
// collect values for primaries and secondaries
for (Int_t iTrack = 0; iTrack < fESD->GetNumberOfTracks(); iTrack++)
{
AliESDtrack* track = 0;
- if (fAnalysisMode == AliPWG0Helper::kTPC)
+ if (fAnalysisMode & AliPWG0Helper::kTPC)
track = AliESDtrackCuts::GetTPCOnlyTrack(fESD, iTrack);
- else if (fAnalysisMode == AliPWG0Helper::kTPCITS)
+ else if (fAnalysisMode & AliPWG0Helper::kTPCITS)
track = fESD->GetTrack(iTrack);
if (!track)
if (stack->Particle(label)->GetPDG()->Charge() == 0)
continue;
- if (TMath::Abs(track->Eta()) < 1)
+ if (TMath::Abs(track->Eta()) < 0.8 && track->Pt() > 0.15)
{
if (stack->IsPhysicalPrimary(label))
{
// primary
if (fEsdTrackCutsPrim->AcceptTrack(track))
{
- if (AliESDtrackCuts::GetSigmaToVertex(track) > 900)
- {
- Printf("Track %d has nsigma of %f. Printing track and vertex...", iTrack, AliESDtrackCuts::GetSigmaToVertex(track));
- Float_t b[2];
- Float_t r[3];
- track->GetImpactParameters(b, r);
- Printf("Impact parameter %f %f and resolution: %f %f %f", b[0], b[1], r[0], r[1], r[2]);
- track->Print("");
- if (vtxESD)
- vtxESD->Print();
- }
+// if (AliESDtrackCuts::GetSigmaToVertex(track) > 900)
+// {
+// Printf("Track %d has nsigma of %f. Printing track and vertex...", iTrack, AliESDtrackCuts::GetSigmaToVertex(track));
+// Float_t b[2];
+// Float_t r[3];
+// track->GetImpactParameters(b, r);
+// Printf("Impact parameter %f %f and resolution: %f %f %f", b[0], b[1], r[0], r[1], r[2]);
+// track->Print("");
+// if (vtxESD)
+// vtxESD->Print();
+// }
}
}
else
}
// TODO mem leak in the continue statements above
- if (fAnalysisMode == AliPWG0Helper::kTPC)
+ if (fAnalysisMode & AliPWG0Helper::kTPC)
delete track;
}
}
}
+
+ if (!foundInEta10)
+ eventTriggered = kFALSE;
+ else{
+ //Printf("The event triggered. Its number in file is %d",fESD->GetEventNumberInFile());
+ fNumberEvent->Fill(fESD->GetEventNumberInFile());
+ }
}
else
return;
+ // fill process type
+ Int_t biny = (Int_t) eventTriggered + 2 * (Int_t) eventVertex;
+ // INEL
+ fEventStats->Fill(-6, biny);
+ // NSD
+ if (processType != AliPWG0Helper::kSD)
+ fEventStats->Fill(-5, biny);
+ // SD, ND, DD
+ if (processType == AliPWG0Helper::kND)
+ fEventStats->Fill(-4, biny);
+ if (processType == AliPWG0Helper::kSD)
+ fEventStats->Fill(-3, biny);
+ if (processType == AliPWG0Helper::kDD)
+ fEventStats->Fill(-2, biny);
+
// loop over mc particles
Int_t nPrim = stack->GetNprimary();
Int_t nAccepted = 0;
+ Bool_t oneParticleEvent = kFALSE;
+ for (Int_t iMc = 0; iMc < nPrim; ++iMc)
+ {
+ //Printf("Getting particle %d", iMc);
+ TParticle* particle = stack->Particle(iMc);
+
+ if (!particle)
+ continue;
+
+ if (AliPWG0Helper::IsPrimaryCharged(particle, nPrim) == kFALSE)
+ continue;
+
+ // for INEL > 0, MB Working Group definition use the second option
+ // 1. standard
+ //if (TMath::Abs(particle->Eta()) < 1.0)
+ // 2. MB Working Group definition
+ if (TMath::Abs(particle->Eta()) < 0.8 && particle->Pt() > fPtMin)
+ {
+ oneParticleEvent = kTRUE;
+ fNumberEventMC->Fill(fESD->GetEventNumberInFile());
+ break;
+ }
+ }
+
+ if (TMath::Abs(vtxMC[2]) < 5.5)
+ {
+ if (oneParticleEvent)
+ fEventStats->Fill(102, biny);
+ else
+ fEventStats->Fill(101, biny);
+ }
+
for (Int_t iMc = 0; iMc < nPrim; ++iMc)
{
//Printf("Getting particle %d", iMc);
if (SignOK(particle->GetPDG()) == kFALSE)
continue;
-
- if (fPIDParticles && TMath::Abs(particle->Eta()) < 1.0)
+
+ // for INEL > 0, MB Working Group definition use the second option
+ // 1. standard
+ //if (fPIDParticles && TMath::Abs(particle->Eta()) < 1.0)
+ // 2. MB Working Group definition
+ if (fPIDParticles && TMath::Abs(particle->Eta()) < 0.8 && particle->Pt() > fPtMin)
fPIDParticles->Fill(particle->GetPdgCode());
Float_t eta = particle->Eta();
+ if (fSymmetrize)
+ eta = TMath::Abs(eta);
Float_t thirdDim = -1;
- if (fAnalysisMode == AliPWG0Helper::kSPD)
+ if (fAnalysisMode & AliPWG0Helper::kSPD)
{
if (fFillPhi)
{
//fTemp1->Fill(eta);
//fTemp2->Fill(y);
- fdNdEtaCorrection->FillMCParticle(vtxMC[2], eta, thirdDim, eventTriggered, eventVertex, processType);
+ Int_t processType2 = processType;
+ if (oneParticleEvent)
+ processType2 |= AliPWG0Helper::kOnePart;
+
+ fdNdEtaCorrection->FillMCParticle(vtxMC[2], eta, thirdDim, eventTriggered, eventVertex, processType2);
if (fOption.Contains("process-types"))
{
// non diffractive
if (processType==AliPWG0Helper::kND)
- fdNdEtaCorrectionSpecial[0]->FillMCParticle(vtxMC[2], eta, thirdDim, eventTriggered, eventVertex, processType);
+ fdNdEtaCorrectionSpecial[0]->FillMCParticle(vtxMC[2], eta, thirdDim, eventTriggered, eventVertex, processType2);
// single diffractive
if (processType==AliPWG0Helper::kSD)
- fdNdEtaCorrectionSpecial[1]->FillMCParticle(vtxMC[2], eta, thirdDim, eventTriggered, eventVertex, processType);
+ fdNdEtaCorrectionSpecial[1]->FillMCParticle(vtxMC[2], eta, thirdDim, eventTriggered, eventVertex, processType2);
// double diffractive
if (processType==AliPWG0Helper::kDD)
- fdNdEtaCorrectionSpecial[2]->FillMCParticle(vtxMC[2], eta, thirdDim, eventTriggered, eventVertex, processType);
+ fdNdEtaCorrectionSpecial[2]->FillMCParticle(vtxMC[2], eta, thirdDim, eventTriggered, eventVertex, processType2);
}
if (fOption.Contains("particle-species"))
case 2212: id = 2; break;
default: id = 3; break;
}
- fdNdEtaCorrectionSpecial[id]->FillMCParticle(vtxMC[2], eta, thirdDim, eventTriggered, eventVertex, processType);
+ fdNdEtaCorrectionSpecial[id]->FillMCParticle(vtxMC[2], eta, thirdDim, eventTriggered, eventVertex, processType2);
}
if (eventTriggered)
nAccepted++;
}
- fEventStats->Fill(AliPWG0Helper::GetLastProcessType(), biny);
+ if (AliPWG0Helper::GetLastProcessType() >= -1)
+ fEventStats->Fill(AliPWG0Helper::GetLastProcessType(), biny);
fMultAll->Fill(nAccepted);
if (eventTriggered) {
primCount[i] = kFALSE;
}
+ Int_t nEta05 = 0;
+ Int_t nEta10 = 0;
+ for (Int_t i=0; i<inputCount; ++i)
+ {
+ if (TMath::Abs(etaArr[i]) < 0.5)
+ nEta05++;
+ if (TMath::Abs(etaArr[i]) < 1.0)
+ nEta10++;
+ }
+
for (Int_t i=0; i<inputCount; ++i)
{
Int_t label = labelArr[i];
continue;
}
- fPIDTracks->Fill(particle->GetPdgCode());
+ // for INEL > 0, MB Working Group definition use the second option
+ // 1. standard
+ //if (TMath::Abs(particle->Eta()) < 1.0)
+ // 2. INEL > 0 MB Working Group definition
+ if (TMath::Abs(particle->Eta()) < 0.8 && particle->Pt()>fPtMin)
+ fPIDTracks->Fill(particle->GetPdgCode());
// find particle that is filled in the correction map
// this should be the particle that has been reconstructed
processed++;
// resolutions
- fEtaResolution->Fill(particle->Eta() - etaArr[i]);
-
- if (fAnalysisMode == AliPWG0Helper::kTPC || fAnalysisMode == AliPWG0Helper::kTPCITS)
- if (TMath::Abs(particle->Eta() < 0.9) && particle->Pt() > 0)
- fpTResolution->Fill(particle->Pt(), (particle->Pt() - thirdDimArr[i]) / particle->Pt());
+ if (fSymmetrize)
+ fEtaResolution->Fill(TMath::Abs(particle->Eta()) - etaArr[i]);
+ else
+ fEtaResolution->Fill(particle->Eta() - etaArr[i]);
+
+ if (fAnalysisMode & AliPWG0Helper::kTPC || fAnalysisMode & AliPWG0Helper::kTPCITS){
+ // for INEL > 0, MB Working Group definition use the second option
+ // 1. standard
+ //if (TMath::Abs(particle->Eta() < 0.9) && particle->Pt() > 0)
+ // 2. INEL > 0 MB WOrking Group definition
+ if (TMath::Abs(particle->Eta() < 0.8) && particle->Pt() > 0){
+ fpTResolution->Fill(particle->Pt(), (particle->Pt() - thirdDimArr[i]) / particle->Pt());
+ fpTCorrelation->Fill(particle->Pt(),thirdDimArr[i]);
+ fpTCorrelationShift->Fill(particle->Pt(),particle->Pt()-thirdDimArr[i]);
+ }
+ }
Float_t eta = -999;
Float_t thirdDim = -1;
Bool_t firstIsPrim = stack->IsPhysicalPrimary(label);
- // in case of primary the MC values are filled, otherwise (background) the reconstructed values
- if (label == label2 && firstIsPrim)
+ // in case of same label the MC values are filled, otherwise (background) the reconstructed values
+ if (label == label2)
{
eta = particle->Eta();
+ if (fSymmetrize)
+ eta = TMath::Abs(eta);
- if (fAnalysisMode == AliPWG0Helper::kSPD)
+ if (fAnalysisMode & AliPWG0Helper::kSPD)
{
if (fFillPhi)
{
}
else
{
- if (fAnalysisMode == AliPWG0Helper::kSPD && !fFillPhi)
+ if (fAnalysisMode & AliPWG0Helper::kSPD && !fFillPhi)
{
- thirdDim = inputCount;
+ thirdDim = (fMultAxisEta1) ? nEta10 : inputCount;
}
else
thirdDim = thirdDimArr[i];
}
if (fillTrack)
- fdNdEtaCorrection->FillTrackedParticle(vtxMC[2], eta, thirdDim);
-
+ {
+ Double_t weight = 1.;
+ if (fWeightSecondaries){
+ if (!firstIsPrim){
+ weight = GetSecondaryCorrection(thirdDim);
+ }
+ }
+ fdNdEtaCorrection->FillTrackedParticle(vtxMC[2], eta, thirdDim, weight);
+ fTemp2->Fill(vtxMC[2], eta);
+ }
+
// eta comparison for tracklets with the same label (others are background)
if (label == label2)
{
- fEtaProfile->Fill(particle->Eta(), particle->Eta() - etaArr[i]);
- fEtaCorrelation->Fill(etaArr[i], particle->Eta());
- fEtaCorrelationShift->Fill(particle->Eta(), particle->Eta() - etaArr[i]);
+ Float_t eta2 = particle->Eta();
+ if (fSymmetrize)
+ eta2 = TMath::Abs(eta2);
+
+ fEtaProfile->Fill(eta2, eta2 - etaArr[i]);
+ fEtaCorrelation->Fill(eta2, etaArr[i]);
+ fEtaCorrelationShift->Fill(eta2, eta2 - etaArr[i]);
}
- fdNdEtaAnalysisESD->FillTrack(vtxMC[2], particle->Eta(), thirdDim);
+ if (fSymmetrize)
+ fdNdEtaAnalysisESD->FillTrack(vtxMC[2], TMath::Abs(particle->Eta()), thirdDim);
+ else
+ fdNdEtaAnalysisESD->FillTrack(vtxMC[2], particle->Eta(), thirdDim);
if (fOption.Contains("process-types"))
{
if (eventTriggered && eventVertex)
{
- fdNdEtaAnalysisMC->FillEvent(vtxMC[2], inputCount);
- fdNdEtaAnalysisESD->FillEvent(vtxMC[2], inputCount);
+ Double_t diff = vtxMC[2] - vtx[2];
+ fVertexShift->Fill(diff);
+
+ fVertexCorrelation->Fill(vtxMC[2], vtx[2]);
+ fVertexCorrelationShift->Fill(vtxMC[2], vtxMC[2] - vtx[2], inputCount);
+ fVertexProfile->Fill(vtxMC[2], vtxMC[2] - vtx[2]);
+
+ if (vtxESD->IsFromVertexerZ() && inputCount > 0)
+ fVertexShiftNorm->Fill(diff, vtxESD->GetNContributors());
+ }
+
+ Int_t multAxis = inputCount;
+ if (fMultAxisEta1)
+ multAxis = nEta10;
+
+ if (eventTriggered && eventVertex)
+ {
+ fdNdEtaAnalysisMC->FillEvent(vtxMC[2], multAxis);
+ fdNdEtaAnalysisESD->FillEvent(vtxMC[2], multAxis);
}
- // stuff regarding the vertex reco correction and trigger bias correction
- fdNdEtaCorrection->FillEvent(vtxMC[2], inputCount, eventTriggered, eventVertex, processType);
+ Int_t processType2 = processType;
+ if (oneParticleEvent)
+ processType2 |= AliPWG0Helper::kOnePart;
+
+ // stuff regarding the vertex reco correction and trigger bias correction
+ fdNdEtaCorrection->FillEvent(vtxMC[2], multAxis, eventTriggered, eventVertex, processType2);
if (fOption.Contains("process-types"))
{
// non diffractive
- if (processType == AliPWG0Helper::kND )
- fdNdEtaCorrectionSpecial[0]->FillEvent(vtxMC[2], inputCount, eventTriggered, eventVertex, processType);
+ if (processType == AliPWG0Helper::kND)
+ fdNdEtaCorrectionSpecial[0]->FillEvent(vtxMC[2], multAxis, eventTriggered, eventVertex, processType2);
// single diffractive
if (processType == AliPWG0Helper::kSD)
- fdNdEtaCorrectionSpecial[1]->FillEvent(vtxMC[2], inputCount, eventTriggered, eventVertex, processType);
+ fdNdEtaCorrectionSpecial[1]->FillEvent(vtxMC[2], multAxis, eventTriggered, eventVertex, processType2);
// double diffractive
if (processType == AliPWG0Helper::kDD)
- fdNdEtaCorrectionSpecial[2]->FillEvent(vtxMC[2], inputCount, eventTriggered, eventVertex, processType);
+ fdNdEtaCorrectionSpecial[2]->FillEvent(vtxMC[2], multAxis, eventTriggered, eventVertex, processType2);
}
if (fOption.Contains("particle-species"))
for (Int_t id=0; id<4; id++)
- fdNdEtaCorrectionSpecial[id]->FillEvent(vtxMC[2], inputCount, eventTriggered, eventVertex, processType);
+ fdNdEtaCorrectionSpecial[id]->FillEvent(vtxMC[2], multAxis, eventTriggered, eventVertex, processType2);
if (etaArr)
delete[] etaArr;
fVertexCorrelation = dynamic_cast<TH2F*> (fOutput->FindObject("fVertexCorrelation"));
if (fVertexCorrelation)
fVertexCorrelation->Write();
- fVertexCorrelationShift = dynamic_cast<TH2F*> (fOutput->FindObject("fVertexCorrelationShift"));
+ fVertexCorrelationShift = dynamic_cast<TH3F*> (fOutput->FindObject("fVertexCorrelationShift"));
if (fVertexCorrelationShift)
+ {
+ ((TH2*) fVertexCorrelationShift->Project3D("yx"))->FitSlicesY();
fVertexCorrelationShift->Write();
+ }
fVertexProfile = dynamic_cast<TProfile*> (fOutput->FindObject("fVertexProfile"));
if (fVertexProfile)
fVertexProfile->Write();
fVertexShift = dynamic_cast<TH1F*> (fOutput->FindObject("fVertexShift"));
if (fVertexShift)
fVertexShift->Write();
- fVertexShiftNorm = dynamic_cast<TH1F*> (fOutput->FindObject("fVertexShiftNorm"));
+ fVertexShiftNorm = dynamic_cast<TH2F*> (fOutput->FindObject("fVertexShiftNorm"));
if (fVertexShiftNorm)
+ {
+ fVertexShiftNorm->ProjectionX();
fVertexShiftNorm->Write();
+ }
fEtaCorrelation = dynamic_cast<TH2F*> (fOutput->FindObject("fEtaCorrelation"));
if (fEtaCorrelation)
fEtaCorrelation->Write();
+ fEtaCorrelationAllESD = dynamic_cast<TH2F*> (fOutput->FindObject("fEtaCorrelationAllESD"));
+ if (fEtaCorrelationAllESD)
+ fEtaCorrelationAllESD->Write();
fEtaCorrelationShift = dynamic_cast<TH2F*> (fOutput->FindObject("fEtaCorrelationShift"));
if (fEtaCorrelationShift)
+ {
+ fEtaCorrelationShift->FitSlicesY();
fEtaCorrelationShift->Write();
+ }
fEtaProfile = dynamic_cast<TProfile*> (fOutput->FindObject("fEtaProfile"));
if (fEtaProfile)
fEtaProfile->Write();
fEtaResolution = dynamic_cast<TH1F*> (fOutput->FindObject("fEtaResolution"));
if (fEtaResolution)
fEtaResolution->Write();
+ fpTCorrelation = dynamic_cast<TH2F*> (fOutput->FindObject("fpTCorrelation"));
+ if (fpTCorrelation)
+ fpTCorrelation->Write();
+ fpTCorrelationShift = dynamic_cast<TH2F*> (fOutput->FindObject("fpTCorrelationShift"));
+ if (fpTCorrelationShift)
+ {
+ fpTCorrelationShift->FitSlicesY();
+ fpTCorrelationShift->Write();
+ }
fpTResolution = dynamic_cast<TH2F*> (fOutput->FindObject("fpTResolution"));
if (fpTResolution)
{
fpTResolution->Write();
}
+ fpTCorrelationAllESD = dynamic_cast<TH2F*> (fOutput->FindObject("fpTCorrelationAllESD"));
+ if (fpTCorrelationAllESD)
+ fpTCorrelationAllESD->Write();
+ fpTCorrelationShiftAllESD = dynamic_cast<TH2F*> (fOutput->FindObject("fpTCorrelationShiftAllESD"));
+ if (fpTCorrelationShiftAllESD)
+ {
+ fpTCorrelationShiftAllESD->FitSlicesY();
+ fpTCorrelationShiftAllESD->Write();
+ }
fMultAll = dynamic_cast<TH1F*> (fOutput->FindObject("fMultAll"));
if (fMultAll)
fMultAll->Write();
if (fPIDTracks)
fPIDTracks->Write();
+ fPtMC = dynamic_cast<TH1F*> (fOutput->FindObject("fPtMC"));
+ if (fPtMC)
+ fPtMC->Write();
+ fEtaMC = dynamic_cast<TH1F*> (fOutput->FindObject("fEtaMC"));
+ if (fEtaMC)
+ fEtaMC->Write();
+ fPtESD = dynamic_cast<TH1F*> (fOutput->FindObject("fPtESD"));
+ if (fPtESD)
+ fPtESD->Write();
+ fEtaESD = dynamic_cast<TH1F*> (fOutput->FindObject("fEtaESD"));
+ if (fEtaESD)
+ fEtaESD->Write();
+ fVtxMC = dynamic_cast<TH1F*> (fOutput->FindObject("fVtxMC"));
+ if (fVtxMC)
+ fVtxMC->Write();
+
+ fNumberEventMC = dynamic_cast<TH1F*> (fOutput->FindObject("fNumberEventMC"));
+ if (fNumberEventMC)
+ fNumberEventMC->Write();
+ fNumberEvent = dynamic_cast<TH1F*> (fOutput->FindObject("fNumberEvent"));
+ if (fNumberEvent)
+ fNumberEvent->Write();
+
//fdNdEtaCorrection->DrawHistograms();
- Printf("Writting result to %s", fileName.Data());
+ Printf("Writing result to %s", fileName.Data());
if (fPIDParticles && fPIDTracks)
{
return kTRUE;
}
+Double_t AlidNdEtaCorrectionTask::GetSecondaryCorrection(Double_t pt){
+
+ // getting the data driven correction factor to correct for
+ // the underestimate of secondaries in Pythia
+ // (A. Dainese + J. Otwinowski
+
+ if (pt <= 0.17) return 1.0;
+ if (pt <= 0.4) return GetLinearInterpolationValue(0.17,1.0,0.4,1.07, pt);
+ if (pt <= 0.6) return GetLinearInterpolationValue(0.4,1.07,0.6,1.25, pt);
+ if (pt <= 1.2) return GetLinearInterpolationValue(0.6,1.25,1.2,1.5, pt);
+ return 1.5;
+
+}
+
+Double_t AlidNdEtaCorrectionTask::GetLinearInterpolationValue(Double_t const x1, Double_t const y1, Double_t const x2, Double_t const y2, const Double_t pt)
+{
+
+ //
+ // linear interpolation to be used to get the secondary correction (see AlidNdEtaCorrectionTask::GetSecondaryCorrection)
+ //
+
+ return ((y2-y1)/(x2-x1))*pt+(y2-(((y2-y1)/(x2-x1))*x2));
+}
+
git://git.uio.no / u / mrichter / AliRoot.git / blobdiff