#include "AlidNdEtaTask.h"
#include <TStyle.h>
#include <TSystem.h>
#include <TCanvas.h>
#include <TVector3.h>
#include <TChain.h>
#include <TFile.h>
#include <TH1F.h>
#include <TH1D.h>
#include <TH2F.h>
#include <TH3F.h>
#include <TParticle.h>
#include <TRandom.h>
#include <TNtuple.h>
#include <TObjString.h>
#include <TF1.h>
#include <TGraph.h>
#include <AliLog.h>
#include <AliESDVertex.h>
#include <AliESDEvent.h>
#include <AliStack.h>
#include <AliHeader.h>
#include <AliGenEventHeader.h>
#include <AliMultiplicity.h>
#include <AliAnalysisManager.h>
#include <AliMCEventHandler.h>
#include <AliMCEvent.h>
#include <AliESDInputHandler.h>
#include <AliESDInputHandlerRP.h>
#include <AliESDHeader.h>
#include "AliESDtrackCuts.h"
#include "AliPWG0Helper.h"
#include "AliCorrection.h"
#include "AliCorrectionMatrix3D.h"
#include "dNdEtaAnalysis.h"
#include "AliTriggerAnalysis.h"
#include "AliPhysicsSelection.h"
#ifdef FULLALIROOT
#include "../ITS/AliITSRecPoint.h"
#include "AliCDBManager.h"
#include "AliCDBEntry.h"
#include "AliGeomManager.h"
#include "TGeoManager.h"
#endif
ClassImp(AlidNdEtaTask)
AlidNdEtaTask::AlidNdEtaTask(const char* opt) :
AliAnalysisTaskSE("AlidNdEtaTask"),
fESD(0),
fOutput(0),
fOption(opt),
fAnalysisMode((AliPWG0Helper::AnalysisMode) (AliPWG0Helper::kTPC | AliPWG0Helper::kFieldOn)),
fTrigger(AliTriggerAnalysis::kMB1),
fFillPhi(kFALSE),
fDeltaPhiCut(-1),
fReadMC(kFALSE),
fUseMCVertex(kFALSE),
fOnlyPrimaries(kFALSE),
fUseMCKine(kFALSE),
fSymmetrize(kFALSE),
fMultAxisEta1(kFALSE),
fDiffTreatment(AliPWG0Helper::kMCFlags),
fEsdTrackCuts(0),
fdNdEtaAnalysisESD(0),
fMult(0),
fMultVtx(0),
fEvents(0),
fVertexResolution(0),
fdNdEtaAnalysis(0),
fdNdEtaAnalysisND(0),
fdNdEtaAnalysisNSD(0),
fdNdEtaAnalysisOnePart(0),
fdNdEtaAnalysisTr(0),
fdNdEtaAnalysisTrVtx(0),
fdNdEtaAnalysisTracks(0),
fPartPt(0),
fVertex(0),
fVertexVsMult(0),
fPhi(0),
fRawPt(0),
fEtaPhi(0),
fModuleMap(0),
fDeltaPhi(0),
fDeltaTheta(0),
fFiredChips(0),
fTrackletsVsClusters(0),
fTrackletsVsUnassigned(0),
fStats(0),
fStats2(0),
fPtMin(0.15),
fEta(0x0),
fEtaMC(0x0),
fHistEvents(0),
fHistEventsMC(0),
fTrigEffNum(0),
fTrigEffDen(0),
fVtxEffNum(0),
fVtxEffDen(0),
fVtxTrigEffNum(0),
fVtxTrigEffDen(0)
{
DefineOutput(1, TList::Class());
fZPhi[0] = 0;
fZPhi[1] = 0;
AliLog::SetClassDebugLevel("AlidNdEtaTask", AliLog::kWarning);
}
AlidNdEtaTask::~AlidNdEtaTask()
{
if (fOutput && !AliAnalysisManager::GetAnalysisManager()->IsProofMode()) {
delete fOutput;
fOutput = 0;
}
}
Bool_t AlidNdEtaTask::UserNotify()
{
static Int_t count = 0;
count++;
Printf("Processing %d. file: %s", count, ((TTree*) GetInputData(0))->GetCurrentFile()->GetName());
return kTRUE;
}
void AlidNdEtaTask::ConnectInputData(Option_t *opt)
{
AliAnalysisTaskSE::ConnectInputData(opt);
Printf("AlidNdEtaTask::ConnectInputData called");
AliLog::SetClassDebugLevel("AliMCEvent", AliLog::kWarning - AliLog::kDebug + 1);
#ifdef FULLALIROOT
AliCDBManager::Instance()->SetDefaultStorage("raw://");
AliCDBManager::Instance()->SetRun(0);
AliCDBManager* mgr = AliCDBManager::Instance();
AliCDBEntry* obj = mgr->Get(AliCDBPath("GRP", "Geometry", "Data"));
AliGeomManager::SetGeometry((TGeoManager*) obj->GetObject());
AliGeomManager::GetNalignable("ITS");
AliGeomManager::ApplyAlignObjsFromCDB("ITS");
#endif
}
void AlidNdEtaTask::UserCreateOutputObjects()
{
Printf("AlidNdEtaTask::CreateOutputObjects");
if (fOnlyPrimaries)
Printf("WARNING: Processing only primaries (MC information used). This is for systematical checks only.");
if (fUseMCKine)
Printf("WARNING: Using MC kine information. This is for systematical checks only.");
if (fUseMCVertex)
Printf("WARNING: Replacing vertex by MC vertex. This is for systematical checks only.");
fOutput = new TList;
fOutput->SetOwner();
fdNdEtaAnalysisESD = new dNdEtaAnalysis("fdNdEtaAnalysisESD", "fdNdEtaAnalysisESD", fAnalysisMode);
fOutput->Add(fdNdEtaAnalysisESD);
fMult = new TH1F("fMult", "fMult;Ntracks;Count", 201, -0.5, 200.5);
fOutput->Add(fMult);
fMultVtx = new TH1F("fMultVtx", "fMultVtx;Ntracks;Count", 201, -0.5, 200.5);
fOutput->Add(fMultVtx);
for (Int_t i=0; i<3; ++i)
{
fPartEta[i] = new TH1F(Form("dndeta_check_%d", i), Form("dndeta_check_%d", i), 60, -3, 3);
fPartEta[i]->Sumw2();
fOutput->Add(fPartEta[i]);
}
fEvents = new TH1F("dndeta_check_vertex", "dndeta_check_vertex", 800, -40, 40);
fOutput->Add(fEvents);
fVertexResolution = new TH2F("dndeta_vertex_resolution_z", ";z resolution;#Delta #phi", 1000, 0, 2, 100, 0, 0.2);
fOutput->Add(fVertexResolution);
fPhi = new TH1F("fPhi", "fPhi;#phi in rad.;count", 720, 0, 2 * TMath::Pi());
fOutput->Add(fPhi);
fEtaPhi = new TH2F("fEtaPhi", "fEtaPhi;#eta;#phi in rad.;count", 80, -4, 4, 18*20, 0, 2 * TMath::Pi());
fOutput->Add(fEtaPhi);
fStats = new TH1F("fStats", "fStats", 5, 0.5, 5.5);
fStats->GetXaxis()->SetBinLabel(1, "vertexer 3d");
fStats->GetXaxis()->SetBinLabel(2, "vertexer z");
fStats->GetXaxis()->SetBinLabel(3, "trigger");
fStats->GetXaxis()->SetBinLabel(4, "physics events");
fStats->GetXaxis()->SetBinLabel(5, "physics events after veto");
fOutput->Add(fStats);
fStats2 = new TH2F("fStats2", "fStats2", 7, -0.5, 6.5, 12, -0.5, 11.5);
fStats2->GetXaxis()->SetBinLabel(1, "No trigger");
fStats2->GetXaxis()->SetBinLabel(2, "No Vertex");
fStats2->GetXaxis()->SetBinLabel(3, "Vtx res. too bad");
fStats2->GetXaxis()->SetBinLabel(4, "|z-vtx| > 15");
fStats2->GetXaxis()->SetBinLabel(5, "|z-vtx| > 10");
fStats2->GetXaxis()->SetBinLabel(6, "0 tracklets");
fStats2->GetXaxis()->SetBinLabel(7, "Selected");
fStats2->GetYaxis()->SetBinLabel(1, "n/a");
fStats2->GetYaxis()->SetBinLabel(2, "empty");
fStats2->GetYaxis()->SetBinLabel(3, "BBA");
fStats2->GetYaxis()->SetBinLabel(4, "BBC");
fStats2->GetYaxis()->SetBinLabel(5, "BBA BBC");
fStats2->GetYaxis()->SetBinLabel(6, "BGA");
fStats2->GetYaxis()->SetBinLabel(7, "BGC");
fStats2->GetYaxis()->SetBinLabel(8, "BGA BGC");
fStats2->GetYaxis()->SetBinLabel(9, "BBA BGC");
fStats2->GetYaxis()->SetBinLabel(10, "BGA BBC");
fStats2->GetYaxis()->SetBinLabel(11, "GFO");
fStats2->GetYaxis()->SetBinLabel(12, "NO GFO");
fOutput->Add(fStats2);
fTrackletsVsClusters = new TH2F("fTrackletsVsClusters", ";tracklets;clusters in ITS", 50, -0.5, 49.5, 1000, -0.5, 999.5);
fOutput->Add(fTrackletsVsClusters);
if (fAnalysisMode & AliPWG0Helper::kSPD)
{
fDeltaPhi = new TH1F("fDeltaPhi", "fDeltaPhi;#Delta #phi;Entries", 500, -0.16, 0.16);
fOutput->Add(fDeltaPhi);
fDeltaTheta = new TH1F("fDeltaTheta", "fDeltaTheta;#Delta #theta;Entries", 500, -0.05, 0.05);
fOutput->Add(fDeltaTheta);
fFiredChips = new TH2F("fFiredChips", "fFiredChips;Chips L1 + L2;tracklets", 1201, -0.5, 1201, 50, -0.5, 49.5);
fOutput->Add(fFiredChips);
fTrackletsVsUnassigned = new TH2F("fTrackletsVsUnassigned", ";tracklets;unassigned clusters in L0", 50, -0.5, 49.5, 200, -0.5, 199.5);
fOutput->Add(fTrackletsVsUnassigned);
for (Int_t i=0; i<2; i++)
{
fZPhi[i] = new TH2F(Form("fZPhi_%d", i), Form("fZPhi Layer %d;z (cm);#phi (rad.)", i), 200, -15, 15, 200, 0, TMath::Pi() * 2);
fOutput->Add(fZPhi[i]);
}
fModuleMap = new TH1F("fModuleMap", "fModuleMap;module number;cluster count", 240, -0.5, 239.5);
fOutput->Add(fModuleMap);
}
if (fAnalysisMode & AliPWG0Helper::kTPC || fAnalysisMode & AliPWG0Helper::kTPCITS)
{
fRawPt = new TH1F("fRawPt", "raw pt;p_{T};Count", 2000, 0, 100);
fOutput->Add(fRawPt);
}
fVertex = new TH3F("vertex_check", "vertex_check;x;y;z", 100, -1, 1, 100, -1, 1, 100, -30, 30);
fOutput->Add(fVertex);
fVertexVsMult = new TH3F("fVertexVsMult", "fVertexVsMult;x;y;multiplicity", 100, -1, 1, 100, -1, 1, 100, -0.5, 99.5);
fOutput->Add(fVertexVsMult);
if (fReadMC)
{
fdNdEtaAnalysis = new dNdEtaAnalysis("dndeta", "dndeta", fAnalysisMode);
fOutput->Add(fdNdEtaAnalysis);
fdNdEtaAnalysisND = new dNdEtaAnalysis("dndetaND", "dndetaND", fAnalysisMode);
fOutput->Add(fdNdEtaAnalysisND);
fdNdEtaAnalysisNSD = new dNdEtaAnalysis("dndetaNSD", "dndetaNSD", fAnalysisMode);
fOutput->Add(fdNdEtaAnalysisNSD);
fdNdEtaAnalysisOnePart = new dNdEtaAnalysis("dndetaOnePart", "dndetaOnePart", fAnalysisMode);
fOutput->Add(fdNdEtaAnalysisOnePart);
fdNdEtaAnalysisTr = new dNdEtaAnalysis("dndetaTr", "dndetaTr", fAnalysisMode);
fOutput->Add(fdNdEtaAnalysisTr);
fdNdEtaAnalysisTrVtx = new dNdEtaAnalysis("dndetaTrVtx", "dndetaTrVtx", fAnalysisMode);
fOutput->Add(fdNdEtaAnalysisTrVtx);
fdNdEtaAnalysisTracks = new dNdEtaAnalysis("dndetaTracks", "dndetaTracks", fAnalysisMode);
fOutput->Add(fdNdEtaAnalysisTracks);
fPartPt = new TH1F("dndeta_check_pt", "dndeta_check_pt", 1000, 0, 10);
fPartPt->Sumw2();
fOutput->Add(fPartPt);
}
if (fEsdTrackCuts)
{
fEsdTrackCuts->SetName("fEsdTrackCuts");
fOutput->Add(fEsdTrackCuts);
}
fEta = new TH1D("fEta", "Eta;#eta;count", 80, -4, 4);
fOutput->Add(fEta);
fEtaMC = new TH1D("fEtaMC", "Eta, MC;#eta;count", 80, -4, 4);
fOutput->Add(fEtaMC);
fHistEvents = new TH1D("fHistEvents", "N. of Events;accepted;count", 2, 0, 2);
fOutput->Add(fHistEvents);
fHistEventsMC = new TH1D("fHistEventsMC", "N. of MC Events;accepted;count", 2, 0, 2);
fOutput->Add(fHistEventsMC);
fTrigEffNum = new TH1D("fTrigEffNum", "N. of triggered events", 100,0,100);
fOutput->Add(fTrigEffNum);
fTrigEffDen = new TH1D("fTrigEffDen", "N. of true events", 100,0,100);
fOutput->Add(fTrigEffDen);
fVtxEffNum = new TH1D("fVtxEffNum", "N. of events with vtx", 100,0,100);
fOutput->Add(fVtxEffNum);
fVtxEffDen = new TH1D("fVtxEffDen", "N. of true events", 100,0,100);
fOutput->Add(fVtxEffDen);
fVtxTrigEffNum = new TH1D("fVtxTrigEffNum", "N. of triggered events with vtx", 100,0,100);
fOutput->Add(fVtxTrigEffNum);
fVtxTrigEffDen = new TH1D("fVtxTrigEffDen", "N. of triggered true events", 100,0,100);
fOutput->Add(fVtxTrigEffDen);
}
Bool_t AlidNdEtaTask::IsEventInBinZero()
{
Bool_t isZeroBin = kTRUE;
fESD = (AliESDEvent*) fInputEvent;
AliInputEventHandler* inputHandler = static_cast<AliInputEventHandler*> (AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler());
if (!inputHandler)
{
Printf("ERROR: Could not receive input handler");
return kFALSE;
}
static AliTriggerAnalysis* triggerAnalysis = 0;
if (!triggerAnalysis)
{
AliPhysicsSelection* physicsSelection = dynamic_cast<AliPhysicsSelection*> (inputHandler->GetEventSelection());
if (physicsSelection)
triggerAnalysis = physicsSelection->GetTriggerAnalysis();
}
if (!triggerAnalysis)
{
Printf("ERROR: Could not receive trigger analysis object");
return kFALSE;
}
if (!triggerAnalysis->IsTriggerFired(fESD, fTrigger))
return kFALSE;
const AliMultiplicity* mult = fESD->GetMultiplicity();
if (!mult){
Printf("AlidNdEtaTask::IsBinZero: Can't get multiplicity object from ESDs");
return kFALSE;
}
Int_t ntracklet = mult->GetNumberOfTracklets();
const AliESDVertex * vtxESD = fESD->GetPrimaryVertexSPD();
if(vtxESD) {
if( TMath::Abs(vtxESD->GetZ()) <= 15 ) {
if (vtxESD->IsFromVertexerZ()) {
if (vtxESD->GetDispersion()<=0.02 ) {
if(ntracklet>0) isZeroBin = kFALSE;
}
} else if(ntracklet>0) isZeroBin = kFALSE;
}
}
return isZeroBin;
}
void AlidNdEtaTask::UserExec(Option_t*)
{
fESD = (AliESDEvent*) fInputEvent;
Bool_t eventTriggered = kFALSE;
const AliESDVertex* vtxESD = 0;
PostData(1, fOutput);
if (fESD)
{
AliESDHeader* esdHeader = fESD->GetHeader();
if (!esdHeader)
{
Printf("ERROR: esdHeader could not be retrieved");
return;
}
AliInputEventHandler* inputHandler = (AliInputEventHandler*) AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler();
if (!inputHandler)
{
Printf("ERROR: Could not receive input handler");
return;
}
eventTriggered = inputHandler->IsEventSelected();
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);
AliTriggerAnalysis::V0Decision v0A = triggerAnalysis->V0Trigger(fESD, AliTriggerAnalysis::kASide, kFALSE);
AliTriggerAnalysis::V0Decision v0C = triggerAnalysis->V0Trigger(fESD, AliTriggerAnalysis::kCSide, kFALSE);
Bool_t fastORHW = (triggerAnalysis->SPDFiredChips(fESD, 1) > 0);
Int_t vZero = 0;
if (v0A != AliTriggerAnalysis::kV0Invalid && v0C != AliTriggerAnalysis::kV0Invalid)
{
if (v0A == AliTriggerAnalysis::kV0Empty && v0C == AliTriggerAnalysis::kV0Empty) vZero = 1;
if (v0A == AliTriggerAnalysis::kV0BB && v0C == AliTriggerAnalysis::kV0Empty) vZero = 2;
if (v0A == AliTriggerAnalysis::kV0Empty && v0C == AliTriggerAnalysis::kV0BB) vZero = 3;
if (v0A == AliTriggerAnalysis::kV0BB && v0C == AliTriggerAnalysis::kV0BB) vZero = 4;
if (v0A == AliTriggerAnalysis::kV0BG && v0C == AliTriggerAnalysis::kV0Empty) vZero = 5;
if (v0A == AliTriggerAnalysis::kV0Empty && v0C == AliTriggerAnalysis::kV0BG) vZero = 6;
if (v0A == AliTriggerAnalysis::kV0BG && v0C == AliTriggerAnalysis::kV0BG) vZero = 7;
if (v0A == AliTriggerAnalysis::kV0BB && v0C == AliTriggerAnalysis::kV0BG) vZero = 8;
if (v0A == AliTriggerAnalysis::kV0BG && v0C == AliTriggerAnalysis::kV0BB) vZero = 9;
}
if (vZero == 0)
Printf("VZERO: %d %d %d %d", vZero, eventTriggered, v0A, v0C);
Bool_t filled = kFALSE;
if (!eventTriggered)
{
fStats2->Fill(0.0, vZero);
fStats2->Fill(0.0, (fastORHW) ? 10 : 11);
filled = kTRUE;
}
if (eventTriggered)
fStats->Fill(3);
vtxESD = AliPWG0Helper::GetVertex(fESD, fAnalysisMode);
if (!vtxESD)
{
if (!filled)
{
fStats2->Fill(1, vZero);
fStats2->Fill(1, (fastORHW) ? 10 : 11);
filled = kTRUE;
}
}
else
{
if (!AliPWG0Helper::TestVertex(vtxESD, fAnalysisMode))
{
if (!filled)
{
fStats2->Fill(2, vZero);
fStats2->Fill(2, (fastORHW) ? 10 : 11);
filled = kTRUE;
}
}
Double_t vtx[3];
vtxESD->GetXYZ(vtx);
if (TMath::Abs(vtx[2]) > 15)
{
if (!filled)
{
fStats2->Fill(3, vZero);
fStats2->Fill(3, (fastORHW) ? 10 : 11);
filled = kTRUE;
}
}
if (TMath::Abs(vtx[2]) > 10)
{
if (!filled)
{
fStats2->Fill(4, vZero);
fStats2->Fill(4, (fastORHW) ? 10 : 11);
filled = kTRUE;
}
}
const AliMultiplicity* mult = fESD->GetMultiplicity();
if (!mult)
{
Printf("Returning, no Multiplicity found");
return;
}
if (mult->GetNumberOfTracklets() == 0)
{
if (!filled)
{
fStats2->Fill(5, vZero);
fStats2->Fill(5, (fastORHW) ? 10 : 11);
filled = kTRUE;
}
}
}
if (!filled)
{
fStats2->Fill(6, vZero);
fStats2->Fill(6, (fastORHW) ? 10 : 11);
}
if (eventTriggered)
fStats->Fill(3);
fStats->Fill(5);
vtxESD = AliPWG0Helper::GetVertex(fESD, fAnalysisMode);
Double_t vtx[3];
if (vtxESD)
{
if (strcmp(vtxESD->GetTitle(), "vertexer: Z") == 0)
fVertexResolution->Fill(vtxESD->GetZRes(), vtxESD->GetDispersion());
else
fVertexResolution->Fill(vtxESD->GetZRes(), 0);
{
fVertex->Fill(vtxESD->GetX(), vtxESD->GetY(), vtxESD->GetZ());
}
if (AliPWG0Helper::TestVertex(vtxESD, fAnalysisMode))
{
vtxESD->GetXYZ(vtx);
if (strcmp(vtxESD->GetTitle(), "vertexer: 3D") == 0)
fStats->Fill(1);
if (strcmp(vtxESD->GetTitle(), "vertexer: Z") == 0)
fStats->Fill(2);
if (TMath::Abs(vtx[2]) > 15)
vtxESD = 0;
}
else
vtxESD = 0;
}
AliStack* stack = 0;
TArrayF vtxMC(3);
if (fUseMCVertex || fUseMCKine || fOnlyPrimaries || fReadMC) {
if (!fReadMC) {
Printf("ERROR: fUseMCVertex or fUseMCKine or fOnlyPrimaries set without fReadMC set!");
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;
}
AliHeader* header = mcEvent->Header();
if (!header)
{
AliDebug(AliLog::kError, "Header not available");
return;
}
AliGenEventHeader* genHeader = header->GenEventHeader();
if (!genHeader)
{
AliDebug(AliLog::kError, "Could not retrieve genHeader from Header");
return;
}
genHeader->PrimaryVertex(vtxMC);
if (fUseMCVertex)
vtx[2] = vtxMC[2];
stack = mcEvent->Stack();
if (!stack)
{
AliDebug(AliLog::kError, "Stack not available");
return;
}
}
Int_t inputCount = 0;
Int_t* labelArr = 0;
Float_t* etaArr = 0;
Float_t* thirdDimArr = 0;
if (fAnalysisMode & AliPWG0Helper::kSPD)
{
if (vtxESD)
{
const AliMultiplicity* mult = fESD->GetMultiplicity();
if (!mult)
{
AliDebug(AliLog::kError, "AliMultiplicity not available");
return;
}
Int_t arrayLength = mult->GetNumberOfTracklets();
if (fAnalysisMode & AliPWG0Helper::kSPDOnlyL0)
arrayLength += mult->GetNumberOfSingleClusters();
labelArr = new Int_t[arrayLength];
etaArr = new Float_t[arrayLength];
thirdDimArr = new Float_t[arrayLength];
for (Int_t i=0; i<mult->GetNumberOfTracklets(); ++i)
{
if (fOnlyPrimaries)
if (mult->GetLabel(i, 0) < 0 || mult->GetLabel(i, 0) != mult->GetLabel(i, 1) || !stack->IsPhysicalPrimary(mult->GetLabel(i, 0)))
continue;
Float_t deltaPhi = mult->GetDeltaPhi(i);
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);
Int_t label = mult->GetLabel(i, 0);
Float_t eta = mult->GetEta(i);
Float_t phi = mult->GetPhi(i);
if (phi < 0)
phi += TMath::Pi() * 2;
fPhi->Fill(phi);
if (vtxESD && TMath::Abs(vtx[2]) < 10)
{
fEtaPhi->Fill(eta, phi);
fDeltaPhi->Fill(deltaPhi);
fDeltaTheta->Fill(mult->GetDeltaTheta(i));
}
if (fDeltaPhiCut > 0 && (TMath::Abs(deltaPhi) > fDeltaPhiCut || TMath::Abs(mult->GetDeltaTheta(i)) > fDeltaPhiCut / 0.08 * 0.025))
continue;
if (fUseMCKine)
{
if (label > 0)
{
TParticle* particle = stack->Particle(label);
eta = particle->Eta();
phi = particle->Phi();
}
else
Printf("WARNING: fUseMCKine set without fOnlyPrimaries and no label found");
}
if (fSymmetrize)
eta = TMath::Abs(eta);
etaArr[inputCount] = eta;
labelArr[inputCount] = label;
thirdDimArr[inputCount] = phi;
++inputCount;
}
if (fAnalysisMode & AliPWG0Helper::kSPDOnlyL0)
{
for (Int_t i=0; i<mult->GetNumberOfSingleClusters(); ++i)
{
etaArr[inputCount] = -TMath::Log(TMath::Tan(mult->GetThetaSingle(i)/2.));
labelArr[inputCount] = -1;
thirdDimArr[inputCount] = mult->GetPhiSingle(i);
++inputCount;
}
}
if (!fFillPhi)
{
for (Int_t i=0; i<inputCount; ++i)
thirdDimArr[i] = inputCount;
}
Int_t firedChips = mult->GetNumberOfFiredChips(0) + mult->GetNumberOfFiredChips(1);
fFiredChips->Fill(firedChips, inputCount);
Printf("Accepted %d tracklets (%d fired chips)", inputCount, firedChips);
fTrackletsVsUnassigned->Fill(inputCount, mult->GetNumberOfSingleClusters());
}
}
else if (fAnalysisMode & AliPWG0Helper::kTPC || fAnalysisMode & AliPWG0Helper::kTPCITS)
{
if (!fEsdTrackCuts)
{
AliDebug(AliLog::kError, "fESDTrackCuts not available");
return;
}
Bool_t foundInEta10 = kFALSE;
if (vtxESD)
{
TObjArray* list = fEsdTrackCuts->GetAcceptedTracks(fESD, fAnalysisMode & AliPWG0Helper::kTPC);
Int_t nGoodTracks = list->GetEntries();
Printf("Accepted %d tracks out of %d total ESD tracks", nGoodTracks, fESD->GetNumberOfTracks());
labelArr = new Int_t[nGoodTracks];
etaArr = new Float_t[nGoodTracks];
thirdDimArr = new Float_t[nGoodTracks];
for (Int_t i=0; i<nGoodTracks; i++)
{
AliESDtrack* esdTrack = dynamic_cast<AliESDtrack*> (list->At(i));
if (!esdTrack)
{
AliDebug(AliLog::kError, Form("ERROR: Could not retrieve track %d.", i));
continue;
}
Float_t phi = esdTrack->Phi();
if (phi < 0)
phi += TMath::Pi() * 2;
Float_t eta = esdTrack->Eta();
Int_t label = TMath::Abs(esdTrack->GetLabel());
Float_t pT = esdTrack->Pt();
if (!(fAnalysisMode & AliPWG0Helper::kFieldOn))
pT = 1;
fPhi->Fill(phi);
fEtaPhi->Fill(eta, phi);
if (eventTriggered && vtxESD)
fRawPt->Fill(pT);
if (esdTrack->Pt() < fPtMin)
continue;
if (fOnlyPrimaries)
{
if (label == 0)
continue;
if (stack->IsPhysicalPrimary(label) == kFALSE)
continue;
}
if (TMath::Abs(esdTrack->Eta()) < 0.8 && esdTrack->Pt() > fPtMin){
foundInEta10 = kTRUE;
}
if (fUseMCKine)
{
if (label > 0)
{
TParticle* particle = stack->Particle(label);
eta = particle->Eta();
pT = particle->Pt();
if (TMath::Abs(eta) >=0.8 || pT <= fPtMin){
AliDebug(2,Form("WARNING ************* USING KINE: eta = %f, pt = %f",eta,pT));
}
}
else
Printf("WARNING: fUseMCKine set without fOnlyPrimaries and no label found");
}
if (fSymmetrize)
eta = TMath::Abs(eta);
etaArr[inputCount] = eta;
labelArr[inputCount] = TMath::Abs(esdTrack->GetLabel());
thirdDimArr[inputCount] = pT;
++inputCount;
}
delete list;
}
if (!foundInEta10){
eventTriggered = kFALSE;
fHistEvents->Fill(0.1);
AliDebug(3,"Event rejected");
}
else{
if (eventTriggered){
fHistEvents->Fill(1.1);
AliDebug(3,Form("Event Accepted, with inputcount = %d",inputCount));
}
else{
AliDebug(3,"Event has foundInEta10 but was not triggered");
}
}
}
else
return;
if (eventTriggered)
{
fMult->Fill(inputCount);
fdNdEtaAnalysisESD->FillTriggeredEvent(inputCount);
if (vtxESD)
{
if (strcmp(vtxESD->GetTitle(), "vertexer: 3D") == 0)
fVertexVsMult->Fill(vtxESD->GetX(), vtxESD->GetY(), inputCount);
Int_t part05 = 0;
Int_t part10 = 0;
for (Int_t i=0; i<inputCount; ++i)
{
Float_t eta = etaArr[i];
Float_t thirdDim = thirdDimArr[i];
fEta->Fill(eta);
fdNdEtaAnalysisESD->FillTrack(vtx[2], eta, thirdDim);
if (TMath::Abs(vtx[2]) < 10)
{
fPartEta[0]->Fill(eta);
if (vtx[2] < 0)
fPartEta[1]->Fill(eta);
else
fPartEta[2]->Fill(eta);
}
if (TMath::Abs(eta) < 0.5)
part05++;
if (TMath::Abs(eta) < 1.0)
part10++;
}
Int_t multAxis = inputCount;
if (fMultAxisEta1)
multAxis = part10;
fMultVtx->Fill(multAxis);
fdNdEtaAnalysisESD->FillEvent(vtx[2], multAxis);
if (multAxis > 0)
fEvents->Fill(vtx[2]);
if (fReadMC)
{
for (Int_t i=0; i<inputCount; ++i)
{
Int_t label = labelArr[i];
if (label < 0)
continue;
TParticle* particle = stack->Particle(label);
if (!particle)
{
AliDebug(AliLog::kError, Form("ERROR: Could not retrieve particle %d.", label));
continue;
}
Float_t thirdDim = -1;
if (fAnalysisMode & AliPWG0Helper::kSPD)
{
if (fFillPhi)
{
thirdDim = particle->Phi();
}
else
thirdDim = multAxis;
}
else
thirdDim = particle->Pt();
Float_t eta = particle->Eta();
if (fSymmetrize)
eta = TMath::Abs(eta);
fdNdEtaAnalysisTracks->FillTrack(vtxMC[2], eta, thirdDim);
}
fdNdEtaAnalysisTracks->FillEvent(vtxMC[2], multAxis);
}
}
}
if (etaArr)
delete[] etaArr;
if (labelArr)
delete[] labelArr;
if (thirdDimArr)
delete[] thirdDimArr;
}
if (fReadMC)
{
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;
}
AliGenEventHeader* genHeader = header->GenEventHeader();
if (!genHeader)
{
AliDebug(AliLog::kError, "Could not retrieve genHeader from Header");
return;
}
TArrayF vtxMC(3);
genHeader->PrimaryVertex(vtxMC);
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, Form("Unknown process type %d.", processType));
Int_t nPrim = stack->GetNprimary();
Int_t nAcceptedParticles = 0;
Bool_t oneParticleEvent = kFALSE;
Int_t nMCparticlesinRange = 0;
for (Int_t iMc = 0; iMc < nPrim; ++iMc)
{
if (!stack->IsPhysicalPrimary(iMc))
continue;
TParticle* particle = stack->Particle(iMc);
if (!particle)
continue;
if (particle->GetPDG()->Charge() == 0)
continue;
Float_t eta = particle->Eta();
AliDebug(2,Form("particle %d: eta = %f, pt = %f",iMc,particle->Eta(),particle->Pt()));
if (TMath::Abs(eta) < 0.8 && particle->Pt() > fPtMin){
oneParticleEvent = kTRUE;
nMCparticlesinRange++;
}
if (TMath::Abs(eta) < 1.5)
nAcceptedParticles++;
}
if (TMath::Abs(vtxMC[2]) < 10.){
if (vtxESD){
fVtxEffNum->Fill(nMCparticlesinRange);
}
fVtxEffDen->Fill(nMCparticlesinRange);
if (eventTriggered){
if (vtxESD){
fVtxTrigEffNum->Fill(nMCparticlesinRange);
}
fVtxTrigEffDen->Fill(nMCparticlesinRange);
fTrigEffNum->Fill(nMCparticlesinRange);
}
fTrigEffDen->Fill(nMCparticlesinRange);
}
if (oneParticleEvent){
fHistEventsMC->Fill(1.1);
}
else{
fHistEventsMC->Fill(0.1);
}
for (Int_t iMc = 0; iMc < nPrim; ++iMc)
{
if (!stack->IsPhysicalPrimary(iMc))
continue;
TParticle* particle = stack->Particle(iMc);
if (!particle)
continue;
if (particle->GetPDG()->Charge() == 0)
continue;
Float_t eta = particle->Eta();
if (fSymmetrize)
eta = TMath::Abs(eta);
Float_t thirdDim = -1;
if (fAnalysisMode & AliPWG0Helper::kSPD)
{
if (fFillPhi)
{
thirdDim = particle->Phi();
}
else
thirdDim = nAcceptedParticles;
}
else
thirdDim = particle->Pt();
fdNdEtaAnalysis->FillTrack(vtxMC[2], eta, thirdDim);
if (processType != AliPWG0Helper::kSD)
fdNdEtaAnalysisNSD->FillTrack(vtxMC[2], eta, thirdDim);
if (processType == AliPWG0Helper::kND)
fdNdEtaAnalysisND->FillTrack(vtxMC[2], eta, thirdDim);
if (oneParticleEvent){
AliDebug(3,Form("filling dNdEtaAnalysis object:: vtx = %f, eta = %f, pt = %f",vtxMC[2], eta, thirdDim));
fdNdEtaAnalysisOnePart->FillTrack(vtxMC[2], eta, thirdDim);
}
if (eventTriggered)
{
fdNdEtaAnalysisTr->FillTrack(vtxMC[2], eta, thirdDim);
if (vtxESD)
fdNdEtaAnalysisTrVtx->FillTrack(vtxMC[2], eta, thirdDim);
}
if (TMath::Abs(eta) < 1.0 && particle->Pt() > 0 && particle->P() > 0)
{
Float_t value = 1;
fPartPt->Fill(particle->Pt(), value);
if (TMath::Abs(eta) < 0.8 && particle->Pt() > fPtMin){
fEtaMC->Fill(eta);
}
}
}
fdNdEtaAnalysis->FillEvent(vtxMC[2], nAcceptedParticles);
if (processType != AliPWG0Helper::kSD)
fdNdEtaAnalysisNSD->FillEvent(vtxMC[2], nAcceptedParticles);
if (processType == AliPWG0Helper::kND)
fdNdEtaAnalysisND->FillEvent(vtxMC[2], nAcceptedParticles);
if (oneParticleEvent)
fdNdEtaAnalysisOnePart->FillEvent(vtxMC[2], nAcceptedParticles);
if (eventTriggered)
{
fdNdEtaAnalysisTr->FillEvent(vtxMC[2], nAcceptedParticles);
if (vtxESD)
fdNdEtaAnalysisTrVtx->FillEvent(vtxMC[2], nAcceptedParticles);
}
}
}
void AlidNdEtaTask::Terminate(Option_t *)
{
fOutput = dynamic_cast<TList*> (GetOutputData(1));
if (!fOutput)
Printf("ERROR: fOutput not available");
TH1D* dNdEta = new TH1D("dNdEta","#eta;counts",80, -4, 4);
TH1D* dNdEtaMC = new TH1D("dNdEtaMC","#eta,MC;counts",80, -4, 4);
if (fOutput)
{
fdNdEtaAnalysisESD = dynamic_cast<dNdEtaAnalysis*> (fOutput->FindObject("fdNdEtaAnalysisESD"));
fMult = dynamic_cast<TH1F*> (fOutput->FindObject("fMult"));
fMultVtx = dynamic_cast<TH1F*> (fOutput->FindObject("fMultVtx"));
for (Int_t i=0; i<3; ++i)
fPartEta[i] = dynamic_cast<TH1F*> (fOutput->FindObject(Form("dndeta_check_%d", i)));
fEvents = dynamic_cast<TH1F*> (fOutput->FindObject("dndeta_check_vertex"));
fVertexResolution = dynamic_cast<TH2F*> (fOutput->FindObject("dndeta_vertex_resolution_z"));
fVertex = dynamic_cast<TH3F*> (fOutput->FindObject("vertex_check"));
fVertexVsMult = dynamic_cast<TH3F*> (fOutput->FindObject("fVertexVsMult"));
fPhi = dynamic_cast<TH1F*> (fOutput->FindObject("fPhi"));
fRawPt = dynamic_cast<TH1F*> (fOutput->FindObject("fRawPt"));
fEtaPhi = dynamic_cast<TH2F*> (fOutput->FindObject("fEtaPhi"));
for (Int_t i=0; i<2; ++i)
fZPhi[i] = dynamic_cast<TH2F*> (fOutput->FindObject(Form("fZPhi_%d", i)));
fModuleMap = dynamic_cast<TH1F*> (fOutput->FindObject("fModuleMap"));
fDeltaPhi = dynamic_cast<TH1F*> (fOutput->FindObject("fDeltaPhi"));
fDeltaTheta = dynamic_cast<TH1F*> (fOutput->FindObject("fDeltaTheta"));
fFiredChips = dynamic_cast<TH2F*> (fOutput->FindObject("fFiredChips"));
fTrackletsVsClusters = dynamic_cast<TH2F*> (fOutput->FindObject("fTrackletsVsClusters"));
fTrackletsVsUnassigned = dynamic_cast<TH2F*> (fOutput->FindObject("fTrackletsVsUnassigned"));
fStats = dynamic_cast<TH1F*> (fOutput->FindObject("fStats"));
fStats2 = dynamic_cast<TH2F*> (fOutput->FindObject("fStats2"));
fEsdTrackCuts = dynamic_cast<AliESDtrackCuts*> (fOutput->FindObject("fEsdTrackCuts"));
fEta = dynamic_cast<TH1D*>(fOutput->FindObject("fEta"));
fEtaMC = dynamic_cast<TH1D*>(fOutput->FindObject("fEtaMC"));
fHistEvents = dynamic_cast<TH1D*>(fOutput->FindObject("fHistEvents"));
fHistEventsMC = dynamic_cast<TH1D*>(fOutput->FindObject("fHistEventsMC"));
fVtxEffDen = dynamic_cast<TH1D*>(fOutput->FindObject("fVtxEffDen"));
fVtxEffNum = dynamic_cast<TH1D*>(fOutput->FindObject("fVtxEffNum"));
fTrigEffDen = dynamic_cast<TH1D*>(fOutput->FindObject("fTrigEffDen"));
fTrigEffNum = dynamic_cast<TH1D*>(fOutput->FindObject("fTrigEffNum"));
fVtxTrigEffDen = dynamic_cast<TH1D*>(fOutput->FindObject("fVtxTrigEffDen"));
fVtxTrigEffNum = dynamic_cast<TH1D*>(fOutput->FindObject("fVtxTrigEffNum"));
if (!fHistEvents){
AliError("fHistEvents not found, impossible to determine the corresponding dNdEta distribution for the control file");
}
else Printf("Events selected = %f", fHistEvents->GetBinContent(fHistEvents->GetXaxis()->FindBin(1.1)));
if (!fHistEventsMC){
AliError("fHistEventsMC not found, impossible to determine the corresponding dNdEta distribution for the control file");
}
else Printf("Events selected frm MC = %f", fHistEventsMC->GetBinContent(fHistEventsMC->GetXaxis()->FindBin(1.1)));
if (!fEta){
AliError("fEta not found, impossible to determine the corresponding dNdEta distribution for the control file");
}
if (!fEtaMC){
AliError("fEtaMC not found, impossible to determine the corresponding dNdEta distribution for the control file");
}
Float_t nevents = 0;
if (fHistEvents) nevents = fHistEvents->GetBinContent(fHistEvents->GetXaxis()->FindBin(1.1));
Float_t neventsMC = 0;
if (fHistEventsMC) neventsMC = fHistEventsMC->GetBinContent(fHistEventsMC->GetXaxis()->FindBin(1.1));
if (fEta && fEtaMC){
for (Int_t ibin = 1; ibin <= fEta->GetNbinsX(); ibin++){
Float_t eta =0;
Float_t etaerr =0;
Float_t etaMC =0;
Float_t etaerrMC =0;
if (fEta && nevents > 0) {
eta = fEta->GetBinContent(ibin)/nevents/fEta->GetBinWidth(ibin);
etaerr = fEta->GetBinError(ibin)/nevents/fEta->GetBinWidth(ibin);
dNdEta->SetBinContent(ibin,eta);
dNdEta->SetBinError(ibin,etaerr);
}
if (fEtaMC && neventsMC > 0) {
etaMC = fEtaMC->GetBinContent(ibin)/neventsMC/fEtaMC->GetBinWidth(ibin);
etaerrMC = fEtaMC->GetBinError(ibin)/neventsMC/fEtaMC->GetBinWidth(ibin);
dNdEtaMC->SetBinContent(ibin,etaMC);
dNdEtaMC->SetBinError(ibin,etaerrMC);
}
}
}
new TCanvas("eta", " eta ",50, 50, 550, 550) ;
if (fEta) fEta->Draw();
new TCanvas("etaMC", " etaMC ",50, 50, 550, 550) ;
if (fEtaMC) fEtaMC->Draw();
new TCanvas("dNdEta", "#eta;dNdEta ",50, 50, 550, 550) ;
dNdEta->Draw();
new TCanvas("dNdEtaMC", "#eta,MC;dNdEta ",50, 50, 550, 550) ;
dNdEtaMC->Draw();
new TCanvas("Events", "Events;Events ",50, 50, 550, 550) ;
if (fHistEvents) fHistEvents->Draw();
new TCanvas("Events, MC", "Events, MC;Events ",50, 50, 550, 550) ;
if (fHistEventsMC) fHistEventsMC->Draw();
TFile* outputFileCheck = new TFile("histogramsCheck.root", "RECREATE");
if (fEta) fEta->Write();
if (fEtaMC) fEtaMC->Write();
dNdEta->Write();
dNdEtaMC->Write();
outputFileCheck->Write();
outputFileCheck->Close();
}
if (!fdNdEtaAnalysisESD)
{
AliDebug(AliLog::kError, "ERROR: fdNdEtaAnalysisESD not available");
}
else
{
if (fMult && fMultVtx)
{
new TCanvas;
fMult->Draw();
fMultVtx->SetLineColor(2);
fMultVtx->Draw("SAME");
}
if (fFiredChips)
{
new TCanvas;
fFiredChips->Draw("COLZ");
}
if (fPartEta[0] && fPartEta[1] && fPartEta[2] && fEvents && fMultVtx && fMult)
{
Int_t events1 = (Int_t) fEvents->Integral(fEvents->GetXaxis()->FindBin(-9.9), fEvents->GetXaxis()->FindBin(-0.001));
Int_t events2 = (Int_t) fEvents->Integral(fEvents->GetXaxis()->FindBin(0.001), fEvents->GetXaxis()->FindBin(9.9));
Printf("%d events with vertex used", events1 + events2);
Printf("%d total events with vertex, %d total triggered events, %d triggered events with 0 multiplicity", (Int_t) fMultVtx->Integral(), (Int_t) fMult->Integral(), (Int_t) fMult->GetBinContent(1));
if (events1 > 0 && events2 > 0)
{
fPartEta[0]->Scale(1.0 / (events1 + events2));
fPartEta[1]->Scale(1.0 / events1);
fPartEta[2]->Scale(1.0 / events2);
for (Int_t i=0; i<3; ++i)
fPartEta[i]->Scale(1.0 / fPartEta[i]->GetBinWidth(1));
new TCanvas("control", "control", 500, 500);
for (Int_t i=0; i<3; ++i)
{
fPartEta[i]->SetLineColor(i+1);
fPartEta[i]->Draw((i==0) ? "" : "SAME");
}
}
}
if (fEvents)
{
new TCanvas("control3", "control3", 500, 500);
fEvents->Draw();
}
if (fStats2)
{
new TCanvas;
fStats2->Draw("TEXT");
}
TFile* fout = new TFile("analysis_esd_raw.root", "RECREATE");
if (fdNdEtaAnalysisESD)
fdNdEtaAnalysisESD->SaveHistograms();
if (fEsdTrackCuts)
fEsdTrackCuts->SaveHistograms("esd_track_cuts");
if (fMult)
fMult->Write();
if (fMultVtx)
fMultVtx->Write();
for (Int_t i=0; i<3; ++i)
if (fPartEta[i])
fPartEta[i]->Write();
if (fEvents)
fEvents->Write();
if (fVertexResolution)
{
fVertexResolution->Write();
fVertexResolution->ProjectionX();
fVertexResolution->ProjectionY();
}
if (fDeltaPhi)
fDeltaPhi->Write();
if (fDeltaTheta)
fDeltaTheta->Write();
if (fPhi)
fPhi->Write();
if (fRawPt)
fRawPt->Write();
if (fEtaPhi)
fEtaPhi->Write();
for (Int_t i=0; i<2; ++i)
if (fZPhi[i])
fZPhi[i]->Write();
if (fModuleMap)
fModuleMap->Write();
if (fFiredChips)
fFiredChips->Write();
if (fTrackletsVsClusters)
fTrackletsVsClusters->Write();
if (fTrackletsVsUnassigned)
fTrackletsVsUnassigned->Write();
if (fStats)
fStats->Write();
if (fStats2)
fStats2->Write();
if (fVertex)
fVertex->Write();
if (fVertexVsMult)
fVertexVsMult->Write();
if (fVtxEffDen) fVtxEffDen->Write();
else Printf("fVtxEffDen not found");
if (fVtxEffNum) fVtxEffNum->Write();
else Printf("fVtxEffNum not found");
if (fVtxTrigEffNum) fVtxTrigEffNum->Write();
else Printf("fVtxTrigEffNum not found");
if (fVtxTrigEffDen) fVtxTrigEffDen->Write();
else Printf("fVtxTrigEffDen not found");
if (fTrigEffNum) fTrigEffNum->Write();
else Printf("fTrigEffNum not found");
if (fTrigEffDen) fTrigEffDen->Write();
else Printf("fTrigEffDen not found");
fout->Write();
fout->Close();
Printf("Writting result to analysis_esd_raw.root");
}
if (fReadMC)
{
if (fOutput)
{
fdNdEtaAnalysis = dynamic_cast<dNdEtaAnalysis*> (fOutput->FindObject("dndeta"));
fdNdEtaAnalysisND = dynamic_cast<dNdEtaAnalysis*> (fOutput->FindObject("dndetaND"));
fdNdEtaAnalysisNSD = dynamic_cast<dNdEtaAnalysis*> (fOutput->FindObject("dndetaNSD"));
fdNdEtaAnalysisOnePart = dynamic_cast<dNdEtaAnalysis*> (fOutput->FindObject("dndetaOnePart"));
fdNdEtaAnalysisTr = dynamic_cast<dNdEtaAnalysis*> (fOutput->FindObject("dndetaTr"));
fdNdEtaAnalysisTrVtx = dynamic_cast<dNdEtaAnalysis*> (fOutput->FindObject("dndetaTrVtx"));
fdNdEtaAnalysisTracks = dynamic_cast<dNdEtaAnalysis*> (fOutput->FindObject("dndetaTracks"));
fPartPt = dynamic_cast<TH1F*> (fOutput->FindObject("dndeta_check_pt"));
}
if (!fdNdEtaAnalysis || !fdNdEtaAnalysisTr || !fdNdEtaAnalysisTrVtx || !fPartPt)
{
AliDebug(AliLog::kError, Form("ERROR: Histograms not available %p %p", (void*) fdNdEtaAnalysis, (void*) fPartPt));
return;
}
fdNdEtaAnalysis->Finish(0, -1, AlidNdEtaCorrection::kNone);
fdNdEtaAnalysisND->Finish(0, -1, AlidNdEtaCorrection::kNone);
fdNdEtaAnalysisNSD->Finish(0, -1, AlidNdEtaCorrection::kNone);
fdNdEtaAnalysisOnePart->Finish(0, -1, AlidNdEtaCorrection::kNone);
fdNdEtaAnalysisTr->Finish(0, -1, AlidNdEtaCorrection::kNone);
fdNdEtaAnalysisTrVtx->Finish(0, -1, AlidNdEtaCorrection::kNone);
fdNdEtaAnalysisTracks->Finish(0, -1, AlidNdEtaCorrection::kNone);
Int_t events = (Int_t) fdNdEtaAnalysis->GetData()->GetEventCorrection()->GetMeasuredHistogram()->Integral();
fPartPt->Scale(1.0/events);
fPartPt->Scale(1.0/fPartPt->GetBinWidth(1));
TFile* fout = new TFile("analysis_mc.root","RECREATE");
fdNdEtaAnalysis->SaveHistograms();
fdNdEtaAnalysisND->SaveHistograms();
fdNdEtaAnalysisNSD->SaveHistograms();
fdNdEtaAnalysisOnePart->SaveHistograms();
fdNdEtaAnalysisTr->SaveHistograms();
fdNdEtaAnalysisTrVtx->SaveHistograms();
fdNdEtaAnalysisTracks->SaveHistograms();
if (fPartPt)
fPartPt->Write();
fout->Write();
fout->Close();
Printf("Writting result to analysis_mc.root");
if (fPartPt)
{
new TCanvas("control2", "control2", 500, 500);
fPartPt->Draw();
}
}
}