#include "AliESDInputHandler.h"
#include "AliESDInputHandlerRP.h"
+#include "AliESDVZERO.h"
#include "../ITS/AliITSRecPoint.h"
#include "AliCDBPath.h"
#include "AliCDBManager.h"
#include "AliTrackReference.h"
#include "AliGenHijingEventHeader.h"
+#include "AliGenDPMjetEventHeader.h"
#include "AliLog.h"
fMCCentralityBin(AliAnalysisTaskSPDdNdEta::kall),
fCentrLowLim(0),
fCentrUpLim(0),
- fCentrEst(""),
+ fCentrEst(kFALSE),
fMinClMultLay2(0),
+ fMaxClMultLay2(0),
+ fMinMultV0(0),
+ fVtxLim(0),
+ fPartSpecies(0),
fPhiWindow(0),
fThetaWindow(0),
fMultReco(0),
+ fV0Ampl(0),
fHistSPDRAWMultvsZ(0),
fHistSPDRAWMultvsZTriggCentrEvts(0),
fHistSPDRAWMultvsZCentrEvts(0),
fHistSPDphicl2(0),
fHistBkgCorrDen(0),
+ fHistReconstructedProtons(0),
+ fHistReconstructedKaons(0),
+ fHistReconstructedPions(0),
+ fHistReconstructedOthers(0),
+ fHistReconstructedSec(0),
fHistBkgCorrDenPrimGen(0),
+ fHistBkgCombLabels(0),
fHistBkgCorrNum(0),
fHistAlgEffNum(0),
fHistNonDetectableCorrDen(0),
fHistNonDetectableCorrNum(0),
+ fHistProtons(0),
+ fHistKaons(0),
+ fHistPions(0),
+ fHistOthers(0),
fHistAllPrimaries(0),
fHistTrackCentrEvts(0),
fHistTrackTrigCentrEvts(0),
fHistRecvsGenImpactPar(0),
fHistMCNpart(0),
- fHistdPhiPP(0),
- fHistdPhiSS(0),
- fHistdPhiComb(0),
+ fHistdPhidThetaPP(0),
+ fHistdPhidThetaSS(0),
+ fHistdPhidThetaComb(0),
fHistDeVtx(0)
man->SetDefaultStorage("alien://Folder=/alice/data/2010/OCDB");
if (fUseMC) man->SetDefaultStorage("alien://Folder=/alice/simulation/2008/v4-15-Release/Residual");
else man->SetDefaultStorage("alien://Folder=/alice/data/2010/OCDB");
- man->SetRun(137045);
+ man->SetRun(137161);
AliCDBEntry* obj = man->Get(AliCDBPath("GRP", "Geometry", "Data"));
AliGeomManager::SetGeometry((TGeoManager*) obj->GetObject());
AliGeomManager::GetNalignable("ITS");
fOutput = new TList();
fOutput->SetOwner();
- Int_t nBinVtx = 20;
+ Int_t nBinVtx = 40;
Double_t MaxVtx = 20.;
Int_t nBinMultCorr = 200;
// Data to be corrected
// ...event level
+
+ fV0Ampl = new TH1F("fV0Ampl","",500,0.,30000);
+ fOutput->Add(fV0Ampl);
+
fHistSPDRAWMultvsZ = new TH2F("fHistSPDRAWMultvsZ", "",nBinMultCorr,binLimMultCorr,nBinVtx,-MaxVtx,MaxVtx);
fHistSPDRAWMultvsZ->GetXaxis()->SetTitle("Tracklet multiplicity");
fHistSPDRAWMultvsZ->GetYaxis()->SetTitle("z_{SPDvtx} [cm]");
// Track level correction histograms
fHistBkgCorrDen = new TH2F("fHistBkgCorrDen","",nBinEtaCorr,binLimEtaCorr,nBinVtx,-MaxVtx,MaxVtx);
fOutput->Add(fHistBkgCorrDen);
-
+
+ fHistReconstructedProtons = new TH2F("fHistReconstructedProtons","",nBinEtaCorr,binLimEtaCorr,nBinVtx,-MaxVtx,MaxVtx);
+ fOutput->Add(fHistReconstructedProtons);
+ fHistReconstructedKaons = new TH2F("fHistReconstructedKaons","",nBinEtaCorr,binLimEtaCorr,nBinVtx,-MaxVtx,MaxVtx);
+ fOutput->Add(fHistReconstructedKaons);
+ fHistReconstructedPions = new TH2F("fHistReconstructedPions","",nBinEtaCorr,binLimEtaCorr,nBinVtx,-MaxVtx,MaxVtx);
+ fOutput->Add(fHistReconstructedPions);
+ fHistReconstructedOthers = new TH2F("fHistReconstructedOthers","",nBinEtaCorr,binLimEtaCorr,nBinVtx,-MaxVtx,MaxVtx);
+ fOutput->Add(fHistReconstructedOthers);
+
fHistBkgCorrDenPrimGen = new TH2F("fHistBkgCorrDenPrimGen","",nBinEtaCorr,binLimEtaCorr,nBinVtx,-MaxVtx,MaxVtx);
fOutput->Add(fHistBkgCorrDenPrimGen);
+ fHistBkgCombLabels = new TH2F("fHistBkgCombLabels","",nBinEtaCorr,binLimEtaCorr,nBinVtx,-MaxVtx,MaxVtx);
+ fOutput->Add(fHistBkgCombLabels);
+
if (fTR) {
fHistBkgCorrNum = new TH2F("fHistBkgCorrNum","",nBinEtaCorr,binLimEtaCorr,nBinVtx,-MaxVtx,MaxVtx);
fOutput->Add(fHistBkgCorrNum);
fHistNonDetectableCorrNum = new TH2F("fHistNonDetectableCorrNum","",nBinEtaCorr,binLimEtaCorr,nBinVtx,-MaxVtx,MaxVtx);
fOutput->Add(fHistNonDetectableCorrNum);
-
+
+ fHistProtons = new TH2F("fHistProtons","",nBinEtaCorr,binLimEtaCorr,nBinVtx,-MaxVtx,MaxVtx);
+ fOutput->Add(fHistProtons);
+ fHistKaons = new TH2F("fHistKaons","",nBinEtaCorr,binLimEtaCorr,nBinVtx,-MaxVtx,MaxVtx);
+ fOutput->Add(fHistKaons);
+ fHistPions = new TH2F("fHistPions","",nBinEtaCorr,binLimEtaCorr,nBinVtx,-MaxVtx,MaxVtx);
+ fOutput->Add(fHistPions);
+ fHistOthers = new TH2F("fHistOthers","",nBinEtaCorr,binLimEtaCorr,nBinVtx,-MaxVtx,MaxVtx);
+ fOutput->Add(fHistOthers);
+ fHistReconstructedSec = new TH2F("fHistReconstructedSec","",nBinEtaCorr,binLimEtaCorr,nBinVtx,-MaxVtx,MaxVtx);
+ fOutput->Add(fHistReconstructedSec);
+
+
fHistAllPrimaries = new TH2F("fHistAllPrimaries","",nBinEtaCorr,binLimEtaCorr,nBinVtx,-MaxVtx,MaxVtx);
fOutput->Add(fHistAllPrimaries);
fHistMCNpart->GetYaxis()->SetTitle("Entries");
fOutput->Add(fHistMCNpart);
- fHistdPhiPP = new TH1F("fHistdPhiPP","",400,-0.1,.1);
- fOutput->Add(fHistdPhiPP);
- fHistdPhiSS = new TH1F("fHistdPhiSS","",400,-0.1,.1);
- fOutput->Add(fHistdPhiSS);
- fHistdPhiComb = new TH1F("fHistdPhiComb","",400,-0.1,.1);
- fOutput->Add(fHistdPhiComb);
+ fHistdPhidThetaPP = new TH2F("fHistdPhidThetaPP","",2000,-1.,1.,1000,-0.25,.25);
+ fOutput->Add(fHistdPhidThetaPP);
+ fHistdPhidThetaSS = new TH2F("fHistdPhidThetaSS","",2000,-1.,1.,1000,-0.25,.25);
+ fOutput->Add(fHistdPhidThetaSS);
+ fHistdPhidThetaComb = new TH2F("fHistdPhidThetaComb","",2000,-1.,1.,1000,-0.25,.25);
+ fOutput->Add(fHistdPhidThetaComb);
fHistDeVtx = new TH2F("fHistDeVtx","",80,-20.,20.,5000,-0.5,0.5);
fHistDeVtx->GetXaxis()->SetTitle("z_{MCvtx} [cm]");
AliMagF* field = (AliMagF*)TGeoGlobalMagField::Instance()->GetField();
if (!field && !fmyESD->InitMagneticField()) {Printf("Failed to initialize the B field\n");return;}
+
// Trigger selection
Bool_t eventTriggered = kTRUE;
static AliTriggerAnalysis* triggerAnalysis = 0;
// Centrality selection
Bool_t eventInCentralityBin = kFALSE;
- // Centrality selection
- AliESDCentrality *centrality = fmyESD->GetCentrality();
+/* AliESDCentrality *centrality = fmyESD->GetCentrality();
if (fCentrEst=="") eventInCentralityBin = kTRUE;
else {
if(!centrality) {
if (centrality->IsEventInCentralityClass(fCentrLowLim,fCentrUpLim,fCentrEst.Data())) eventInCentralityBin = kTRUE;
}
}
+*/
+ if (fCentrEst) {
+ AliESDVZERO* esdV0 = fmyESD->GetVZEROData();
+ Float_t multV0A=esdV0->GetMTotV0A();
+ Float_t multV0C=esdV0->GetMTotV0C();
+ fV0Ampl->Fill(multV0A+multV0C);
+ if (multV0A+multV0C>=fMinMultV0) eventInCentralityBin = kTRUE;
+ } else if (!fCentrEst) {
+ eventInCentralityBin = kTRUE;
+ }
+ const AliMultiplicity* multESD = fmyESD->GetMultiplicity();
// ESD vertex
Bool_t eventWithVertex = kFALSE;
const AliESDVertex* vtxESD = fmyESD->GetVertex();
+ const AliESDVertex* vtxTPC = fmyESD->GetPrimaryVertexTPC();
Double_t esdvtx[3];
vtxESD->GetXYZ(esdvtx);
Int_t nContrib = vtxESD->GetNContributors();
- if (nContrib>0) {
-// if (strcmp(vtxESD->GetTitle(),"vertexer: 3D") == 0)
- fHistSPDvtx->Fill(esdvtx[2]);
- eventWithVertex = kTRUE;
+ Int_t nContribTPC = vtxTPC->GetNContributors();
+ if (nContrib>0&&nContribTPC>0) {
+ if (vtxESD->GetDispersion()<0.04)
+ if (vtxESD->GetZRes()<0.25)
+ if (nContribTPC>(-10.+0.25*multESD->GetNumberOfITSClusters(0))) {
+ fHistSPDvtx->Fill(esdvtx[2]);
+ if (TMath::Abs(esdvtx[2])<fVtxLim)
+ eventWithVertex = kTRUE;
+ }
}
// Reconstructing or loading tracklets...
- const AliMultiplicity* multESD = fmyESD->GetMultiplicity();
Int_t multSPD = multESD->GetNumberOfTracklets();
Int_t nSingleCl1 = multESD->GetNumberOfSingleClusters();
Int_t multSPDcl1 = nSingleCl1 + multSPD;
Float_t trackletCoord[multSPDcl1][5];
Float_t trackletLab[multSPDcl1][2];
-
+ Bool_t eventSelected = kFALSE;
// Event selection: in centrality bin, triggered with vertex
- if (eventTriggered&&eventWithVertex&&eventInCentralityBin&&multSPDcl2>fMinClMultLay2) {
+ if (eventTriggered&&eventWithVertex&&eventInCentralityBin&&multSPDcl2>fMinClMultLay2&&multSPDcl2<fMaxClMultLay2) {
+ eventSelected = kTRUE;
fHistSPDmultcl1->Fill(multSPDcl1);
fHistSPDmultcl2->Fill(multSPDcl2);
fHistSPDmultcl1vscl2->Fill(multSPDcl1,multSPDcl2);
return;
}
+ Float_t impactParameter = 0.;
+ Int_t npart = 0;
+
AliGenHijingEventHeader* hijingHeader = dynamic_cast<AliGenHijingEventHeader*>(fMCEvent->Header()->GenEventHeader());
- if (!hijingHeader) {
- Printf("Unknown header type. \n");
- return ;
- }
- Float_t impactParameter = hijingHeader->ImpactParameter();
+ AliGenDPMjetEventHeader* dpmHeader = dynamic_cast<AliGenDPMjetEventHeader*>(fMCEvent->Header()->GenEventHeader());
+
+ if (hijingHeader) impactParameter = hijingHeader->ImpactParameter();
+ else if (dpmHeader) impactParameter = dpmHeader->ImpactParameter();
+
Bool_t IsEventInMCCentralityBin = kFALSE;
switch (fMCCentralityBin) {
fHistMCvtxx->Fill(vtxMC[0]);
fHistMCvtxy->Fill(vtxMC[1]);
fHistMCvtxz->Fill(vtxMC[2]);
+
// Printf("Impact parameter gen: %f", impactParameter);
- Int_t npart = hijingHeader->TargetParticipants()+hijingHeader->ProjectileParticipants();
+ if (hijingHeader) npart = hijingHeader->TargetParticipants()+hijingHeader->ProjectileParticipants();
+ else if (dpmHeader)npart = dpmHeader->TargetParticipants()+dpmHeader->ProjectileParticipants();
+
//Rec centrality vs gen centrality
AliESDZDC *zdcRec = fmyESD->GetESDZDC();
Double_t impactParameterZDC = zdcRec->GetImpactParameter();
// Tracks from MC
Int_t multMCCharged = 0;
Int_t multMCChargedEtacut = 0;
- Int_t nMCPart = stack->GetNprimary();
+// Int_t nMCPart = stack->GetNprimary();
+ Int_t nMCPart = stack->GetNtrack(); // decay products of D and B mesons are also primaries and produced in HIJING during transport
Float_t* etagen = new Float_t[nMCPart];
Int_t* stackIndexOfPrimaryParts = new Int_t[nMCPart];
Bool_t* reconstructedPrimaryPart = new Bool_t[nMCPart];
Bool_t* detectedPrimaryPartLay1 = new Bool_t[nMCPart];
Bool_t* detectedPrimaryPartLay2 = new Bool_t[nMCPart];
Bool_t* detectablePrimaryPart = new Bool_t[nMCPart];
- Bool_t* primCounted = new Bool_t[nMCPart];
TTree* tRefTree;
if (fTR) {
// Loop over MC particles
for (Int_t imc=0; imc<nMCPart; imc++) {
AliMCParticle *mcpart = (AliMCParticle*)fMCEvent->GetTrack(imc);
+
Bool_t isPrimary = stack->IsPhysicalPrimary(imc);
- if (!isPrimary) continue;
+ if (!isPrimary) continue;
if (mcpart->Charge() == 0) continue;
Float_t theta = mcpart->Theta();
if (theta==0 || theta==TMath::Pi()) continue;
detectedPrimaryPartLay1[multMCCharged]=kFALSE;
detectedPrimaryPartLay2[multMCCharged]=kFALSE;
detectablePrimaryPart[multMCCharged]=kFALSE;
- primCounted[multMCCharged]=kFALSE;
if (fTR) {
Int_t nref = mcpart->GetNumberOfTrackReferences();
}
}
}
+ if (eventSelected&&fPartSpecies) {
+ if (TMath::Abs(mcpart->PdgCode())==2212) fHistProtons->Fill(etagen[multMCCharged],vtxMC[2]);
+ else if (TMath::Abs(mcpart->PdgCode())==321) fHistKaons->Fill(etagen[multMCCharged],vtxMC[2]);
+ else if (TMath::Abs(mcpart->PdgCode())==211) fHistPions->Fill(etagen[multMCCharged],vtxMC[2]);
+ else fHistOthers->Fill(etagen[multMCCharged],vtxMC[2]); //includes leptons pdg->11,13
+ }
multMCCharged++;
if (TMath::Abs(eta)<1.4) multMCChargedEtacut++;
} // end of MC particle loop
fHistMCmultEtacut->Fill(multMCChargedEtacut);
// Event selection: in centrality bin, triggered with vertex
- if (eventTriggered&&eventWithVertex&&eventInCentralityBin&&multSPDcl2>fMinClMultLay2) {
+ if (eventTriggered&&eventWithVertex&&eventInCentralityBin&&multSPDcl2>fMinClMultLay2&&multSPDcl2<fMaxClMultLay2) {
fHistDeVtx->Fill(vtxMC[2],vtxMC[2]-esdvtx[2]);
for (Int_t itracklet=0; itracklet<multSPD; ++itracklet) {
- if (TMath::Abs(trackletCoord[itracklet][0])<fPhiWindowAna) {
+ if (TMath::Abs(trackletCoord[itracklet][0])<fPhiWindowAna)
fHistBkgCorrDen->Fill(trackletCoord[itracklet][4],esdvtx[2]);
- if (trackletLab[itracklet][0]==trackletLab[itracklet][1]) {
- Bool_t trakletByPrim = kFALSE;
- for (Int_t imc=0; imc<multMCCharged; imc++) {
- if (trackletLab[itracklet][0]==stackIndexOfPrimaryParts[imc]) {
- if (!primCounted[imc]) {
- primCounted[imc] = kTRUE;
- }
- fHistdPhiPP->Fill(trackletCoord[itracklet][0]);
+ if (trackletLab[itracklet][0]==trackletLab[itracklet][1]) {
+ Bool_t trakletByPrim = kFALSE;
+ for (Int_t imc=0; imc<multMCCharged; imc++) {
+ if (trackletLab[itracklet][0]==stackIndexOfPrimaryParts[imc]) {
+ fHistdPhidThetaPP->Fill(trackletCoord[itracklet][0],trackletCoord[itracklet][1]);
+ if (TMath::Abs(trackletCoord[itracklet][0])<fPhiWindowAna) {
fHistBkgCorrDenPrimGen->Fill(etagen[imc],vtxMC[2]);
- trakletByPrim = kTRUE;
- if (detectedPrimaryPartLay1[imc]&&detectedPrimaryPartLay2[imc]) {
- reconstructedPrimaryPart[imc]=kTRUE; // tracklet by prim and tr (=rp) on both layers
- }
- break;
- }
- }
- if (!trakletByPrim) {
- fHistdPhiSS->Fill(trackletCoord[itracklet][0]);
+ if (fPartSpecies) {
+ if (TMath::Abs(((AliMCParticle*)fMCEvent->GetTrack(stackIndexOfPrimaryParts[imc]))->PdgCode())==2212)
+ fHistReconstructedProtons->Fill(trackletCoord[itracklet][4],esdvtx[2]);
+ else if (TMath::Abs(((AliMCParticle*)fMCEvent->GetTrack(stackIndexOfPrimaryParts[imc]))->PdgCode())==321)
+ fHistReconstructedKaons->Fill(trackletCoord[itracklet][4],esdvtx[2]);
+ else if (TMath::Abs(((AliMCParticle*)fMCEvent->GetTrack(stackIndexOfPrimaryParts[imc]))->PdgCode())==211)
+ fHistReconstructedPions->Fill(trackletCoord[itracklet][4],esdvtx[2]);
+ else fHistReconstructedOthers->Fill(trackletCoord[itracklet][4],esdvtx[2]);
+ }
+ }
+ trakletByPrim = kTRUE;
+ if (detectedPrimaryPartLay1[imc]&&detectedPrimaryPartLay2[imc]) {
+ reconstructedPrimaryPart[imc]=kTRUE; // tracklet by prim and tr (=rp) on both layers
+ }
+ break;
+ }
+ }
+ if (!trakletByPrim) {
+ fHistdPhidThetaSS->Fill(trackletCoord[itracklet][0],trackletCoord[itracklet][1]);
+ if (TMath::Abs(trackletCoord[itracklet][0])<fPhiWindowAna) {
+
fHistBkgCorrDenPrimGen->Fill(trackletCoord[itracklet][4],esdvtx[2]);
+ if (fPartSpecies) {
+ Int_t motherlab = ((AliMCParticle*)fMCEvent->GetTrack((Int_t)trackletLab[itracklet][0]))->GetMother();
+ if (motherlab>-1) {
+ if (TMath::Abs(fMCEvent->GetTrack(motherlab)->PdgCode())==2212) fHistReconstructedProtons->Fill(trackletCoord[itracklet][4],esdvtx[2]);
+ else if (TMath::Abs(fMCEvent->GetTrack(motherlab)->PdgCode())==321) fHistReconstructedKaons->Fill(trackletCoord[itracklet][4],esdvtx[2]);
+ else if (TMath::Abs(fMCEvent->GetTrack(motherlab)->PdgCode())==211) fHistReconstructedPions->Fill(trackletCoord[itracklet][4],esdvtx[2]);
+ else fHistReconstructedOthers->Fill(trackletCoord[itracklet][4],esdvtx[2]);
+ } else fHistReconstructedSec->Fill(trackletCoord[itracklet][4],esdvtx[2]);
+ }
}
- } else {
- fHistdPhiComb->Fill(trackletCoord[itracklet][0]);
+ }
+ } else {
+ fHistdPhidThetaComb->Fill(trackletCoord[itracklet][0],trackletCoord[itracklet][1]);
+ if (TMath::Abs(trackletCoord[itracklet][0])<fPhiWindowAna) {
fHistBkgCorrDenPrimGen->Fill(trackletCoord[itracklet][4],esdvtx[2]);
+ fHistBkgCombLabels->Fill(trackletCoord[itracklet][4],esdvtx[2]);
}
- } // cut dphi
+ }
} // end loop tracklets
for (Int_t imc=0; imc<multMCCharged; imc++) {
Printf("Terminating...");
// AliAnalysisTaskSE::Terminate();
-}
+}
\ No newline at end of file
void SetMCCentralityBin(MCCentralityBin mccentrbin) {fMCCentralityBin=mccentrbin;}
void SetCentralityLowLim(Float_t centrlowlim) {fCentrLowLim=centrlowlim;}
void SetCentralityUpLim(Float_t centruplim) {fCentrUpLim=centruplim;}
- void SetCentralityEst(TString centrest) {fCentrEst=centrest;}
+// void SetCentralityEst(TString centrest) {fCentrEst=centrest;}
+ void SetCentralityEst(Bool_t centrest) {fCentrEst=centrest;}
void SetMinClusterMultLay2(Int_t minClMultLay2=0) {fMinClMultLay2=minClMultLay2;}
+ void SetMaxClusterMultLay2(Int_t maxClMultLay2=0) {fMaxClMultLay2=maxClMultLay2;}
+ void SetMinV0Mult(Int_t minV0Mult=0) {fMinMultV0=minV0Mult;}
+ void SetVertexRange(Float_t vl=7.) {fVtxLim=vl;}
+ void SetPartSpecies(Bool_t partsp = kFALSE) {fPartSpecies=partsp;}
void SetPhiWindow(Float_t w=0.08) {fPhiWindow=w;}
void SetThetaWindow(Float_t w=0.025) {fThetaWindow=w;}
TList *fOutput; // ! output list send on output slot 1
Bool_t fUseMC; // flag to enable the calculation of correction histograms
- Bool_t fPbPb; // flag to analyze PbPb data
AliTriggerAnalysis::Trigger fTrigger;
Bool_t fTR; // to read track references and calculate factors of track to particle correction
Bool_t fRecoTracklets; // flag to recostruct tracklets
MCCentralityBin fMCCentralityBin; // to select MC centrality bin in which corrections are calculated
Float_t fCentrLowLim; // to select centrality bin on data
Float_t fCentrUpLim; // to select centrality bin on data
- TString fCentrEst; // to select centrality estimator
- Int_t fMinClMultLay2; // to select multiplicity class
+// TString fCentrEst; // to select centrality estimator
+ Bool_t fCentrEst; // to select centrality estimator
+ Int_t fMinClMultLay2; // to select multiplicity class
+ Int_t fMaxClMultLay2; // to select multiplicity class
+ Int_t fMinMultV0; // to select centrality class
+ Float_t fVtxLim; // to select vertex range
+ Bool_t fPartSpecies; // to fill correction matrices for each part species (for syst studies)
+
Float_t fPhiWindow; // Search window in phi
Float_t fThetaWindow; // Search window in theta
Float_t fPhiShift; // Phi shift reference value (at 0.5 T)
AliTrackletAlg *fMultReco; // tracklet reconstruction class
+
+ TH1F *fV0Ampl; // ! V0 amplitudes to cut on centrality
TH2F *fHistSPDRAWMultvsZ; // ! data to be corrected
TH2F *fHistSPDRAWMultvsZTriggCentrEvts; // ! data to be corrected
TH2F *fHistSPDRAWMultvsZCentrEvts; // ! data to be corrected
TH1F *fHistSPDmultcl1; // ! cluster inner layer and tracklet check histos
TH1F *fHistSPDmultcl2; // ! cluster inner layer and tracklet check histos
TH2F *fHistSPDmultcl1vscl2; // ! cluster inner layer and tracklet check histos
- TH2F *fHistSPDmultvscl1; // ! cluster inner layer and tracklet check histos
- TH2F *fHistSPDmultvscl2; // ! cluster inner layer and tracklet check histos
+ TH2F *fHistSPDmultvscl1; // ! cluster inner layer and tracklet check histos
+ TH2F *fHistSPDmultvscl2; // ! cluster inner layer and tracklet check histos
TH1F *fHistSPDeta; // ! cluster inner layer and tracklet check histos
TH1F *fHistSPDphi; // ! cluster inner layer and tracklet check histos
TH1F *fHistSPDvtxAnalysis; // ! SPD vertex distributions
TH2F *fHistSPDdePhideTheta; // ! histogram for combinatorial background studies
- TH1F *fHistSPDphicl1; // ! cluster inner layer and tracklet check histos
- TH1F *fHistSPDphicl2; // ! cluster inner layer and tracklet check histos
+ TH1F *fHistSPDphicl1; // ! cluster inner layer and tracklet check histos
+ TH1F *fHistSPDphicl2; // ! cluster inner layer and tracklet check histos
+
+ TH2F* fHistBkgCorrDen; // ! track level correction histograms
+ TH2F* fHistReconstructedProtons; // ! track level correction histograms
+ TH2F* fHistReconstructedKaons; // ! track level correction histograms
+ TH2F* fHistReconstructedPions; // ! track level correction histograms
+ TH2F* fHistReconstructedOthers; // ! track level correction histograms
+ TH2F* fHistReconstructedSec; // ! track level correction histograms
- TH2F* fHistBkgCorrDen; // ! track level correction histograms
TH2F* fHistBkgCorrDenPrimGen; // ! track level correction histograms
+ TH2F* fHistBkgCombLabels; // ! track level correction histograms
TH2F* fHistBkgCorrNum; // ! track level correction histograms
TH2F* fHistAlgEffNum; // ! track level correction histograms
TH2F* fHistNonDetectableCorrDen; // ! track level correction histograms
TH2F* fHistNonDetectableCorrNum; // ! track level correction histograms
+ TH2F* fHistProtons; // ! track level correction histograms
+ TH2F* fHistKaons; // ! track level correction histograms
+ TH2F* fHistPions; // ! track level correction histograms
+ TH2F* fHistOthers; // ! track level correction histograms
TH2F* fHistAllPrimaries; // ! track level correction histograms
TH2F* fHistTrackCentrEvts; // ! track level correction histograms
TH2F* fHistTrackTrigCentrEvts; // ! track level correction histograms
TH2F* fHistTrigCentrEvts; // ! event level correction histograms
TH2F* fHistSelEvts; // ! event level correction histograms
- TH1F* fHistMCmultEtacut; // ! MC distributions
+ TH1F* fHistMCmultEtacut; // ! MC distributions
TH2F* fHistMCmultEtacutvsSPDmultEtacut; // ! MC distributions
- TH2F* fHistMCmultEtacutvsSPDmultcl1; // ! MC distributions
- TH2F* fHistMCmultEtacutvsSPDmultcl2; // ! MC distributions
+ TH2F* fHistMCmultEtacutvsSPDmultcl1; // ! MC distributions
+ TH2F* fHistMCmultEtacutvsSPDmultcl2; // ! MC distributions
TH1F* fHistMCvtxx; // ! MC vertex
TH1F* fHistMCvtxy; // ! MC vertex
TH2F* fHistRecvsGenImpactPar; // ! impact parameter correlation (ZDC est vs gen)
TH1F* fHistMCNpart; // ! distribution of number of participants from MC
- TH1F* fHistdPhiPP; // ! tracklet check histo
- TH1F* fHistdPhiSS; // ! tracklet check histo
- TH1F* fHistdPhiComb; // ! tracklet check histo
+ TH2F* fHistdPhidThetaPP; // ! tracklet check histo
+ TH2F* fHistdPhidThetaSS; // ! tracklet check histo
+ TH2F* fHistdPhidThetaComb; // ! tracklet check histo
TH2F* fHistDeVtx; // ! check histo
AliAnalysisTaskSPDdNdEta(const AliAnalysisTaskSPDdNdEta&); // not implemented
AliAnalysisTaskSPDdNdEta& operator=(const AliAnalysisTaskSPDdNdEta&); // not implemented
- ClassDef(AliAnalysisTaskSPDdNdEta, 3);
+ ClassDef(AliAnalysisTaskSPDdNdEta, 5);
};
-#endif
+#endif
\ No newline at end of file
#include "TLegendEntry.h"
#include "TCanvas.h"
-void correctSPDdNdEta (Char_t* fileRaw, Char_t* fileRawBkgCorr, Char_t* fileAlphaCorr, Char_t* fileMCBkgCorr, Char_t* fileMC,
+Float_t CalculateErrordNdEtaPerPartPair(Float_t a, Float_t b, Float_t da, Float_t db);
+void combscalingfactors (TList *lall, TList *lcomb, TList *lmc, Double_t *scaleNormRange_rot, Double_t *scaleNormRange_labels);
+
+void correctSPDdNdEta (Char_t* fileRaw, Char_t* fileRawBkgCorr,
+ Char_t* fileAlphaCorr, Char_t* fileMCBkgCorr, Char_t* fileMC,
Int_t binVtxStart=5, Int_t binVtxStop=15,
- Double_t cutCorr=10, Double_t cutBkg=1., Double_t cutBkgMC=1., Float_t scaleBkg=0.43, Float_t scaleBkgMC=.505) {
+ Double_t cutCorr=10, Double_t cutBkg=1., Double_t cutBkgMC=1.,
+ Bool_t bkgFromMCLabels = kFALSE,
+ Bool_t changeStrangeness = kFALSE) {
+
+// vtx lim 13-27 for +-7cm
+// vtx lim 15-25 for +-5cm
- Int_t nBinVtx = 20;
+ Int_t nBinVtx = 40;
Double_t MaxVtx = 20.;
const Int_t nBinMultCorr = 200;
TList *l_rawbkg = (TList*)f_rawbkg->Get("cOutput");
TList *l_raw = (TList*)f_raw->Get("cOutput");
+ Double_t scaleBkg = 1.;
+ Double_t scaleBkgMC = 1.;
+ Double_t scaleBkgLabels = 1.;
+
+ cout<<" Background scaling factors for MC: "<<endl;
+ combscalingfactors (l_acorr,l_mcbkgcorr,l_acorr,&scaleBkgMC,&scaleBkgLabels);
+ cout<<" Background scaling factors for data: "<<endl;
+ combscalingfactors (l_raw,l_rawbkg,l_acorr,&scaleBkg,&scaleBkgLabels);
+
//Histogram to be corrected at tracklet level
TH2F *hSPDEta_2Draw = new TH2F(*((TH2F*)l_raw->FindObject("fHistSPDRAWEtavsZ")));
hSPDEta_2Draw->SetNameTitle("SPDEta_2Draw","Reconstructed tracklets");
hCombBkgCorrData->GetXaxis()->SetTitle("#eta");
hCombBkgCorrData->GetYaxis()->SetTitle("z_{SPDvtx} [cm]");
hCombBkgCorrData->Sumw2();
- hCombBkgCorrData->Divide(hCombBkg,hSPDEta_2Draw,1.,1.);
-
+
// Combinatorial background: beta correction from MC
TH2F *hCombBkgMC = (TH2F*)l_mcbkgcorr->FindObject("fHistSPDRAWEtavsZ");
hCombBkgMC->Scale(scaleBkgMC);
hCombBkgCorrMC->GetXaxis()->SetTitle("#eta");
hCombBkgCorrMC->GetYaxis()->SetTitle("z_{SPDvtx} [cm]");
hCombBkgCorrMC->Sumw2();
- hCombBkgCorrMC->Divide(hCombBkgMC,hCombBkgMCDen,1.,1.);
-
- // 1-beta in MC
- TH2F *hCombBkgCorr2MC = new TH2F("backgroundCorr2MC","Combinatorial background correction",nBinEtaCorr,binLimEtaCorr,nBinVtx,-MaxVtx,MaxVtx);
- for (Int_t ieta=0; ieta<nBinEtaCorr; ieta++)
- for (Int_t jz=0; jz<nBinVtx; jz++)
- hCombBkgCorr2MC->SetBinContent(ieta+1,jz+1,1-hCombBkgCorrMC->GetBinContent(ieta+1,jz+1));
+
+ // Background correction from MC labels
+ TH2F *hBkgCombLabels = (TH2F*)l_acorr->FindObject("fHistBkgCombLabels");
+ if (bkgFromMCLabels) {
+ hCombBkgCorrData->Divide(hBkgCombLabels,hSPDEta_2Draw,1.,1.);
+ hCombBkgCorrData->Scale(scaleBkgLabels);
+ } else {
+ hCombBkgCorrData->Divide(hCombBkg,hSPDEta_2Draw,1.,1.);
+ }
-// new TCanvas();
-// hCombBkgCorr2MC->Draw();
+ // 1-beta in DATA
+ TH2F *hCombBkgCorr2Data = new TH2F("backgroundCorr2Data","Combinatorial background correction",nBinEtaCorr,binLimEtaCorr,nBinVtx,-MaxVtx,MaxVtx);
+ for (Int_t ieta=0; ieta<nBinEtaCorr; ieta++) {
+ for (Int_t jz=0; jz<nBinVtx; jz++) {
+ hCombBkgCorr2Data->SetBinContent(ieta+1,jz+1,1-hCombBkgCorrData->GetBinContent(ieta+1,jz+1));
+ hCombBkgCorr2Data->SetBinError(ieta+1,jz+1,hCombBkgCorrData->GetBinError(ieta+1,jz+1));
+ }
+ }
+
// Errors on beta
//..data
for (Int_t ieta=0; ieta<nBinEtaCorr; ieta++) {
for (Int_t jz=0; jz<nBinVtx; jz++) {
hCombBkgCorrData->SetBinError(ieta+1,jz+1,scaleBkg*TMath::Sqrt(hCombBkg->GetBinContent(ieta+1,jz+1))/hSPDEta_2Draw->GetBinContent(ieta+1,jz+1));
+// hCombBkgCorrData->SetBinError(ieta+1,jz+1,hCombBkgCorrData->GetBincontent(ieta+1,jz+1)*TMath::Sqrt(hSPDEta_2Draw->GetBinContent(ieta+1,jz+1))/hSPDEta_2Draw->GetBinContent(ieta+1,jz+1));
}
}
- //..MC--> no, computed in alpha
-
-
// Alpha correction: MC only
TH2F *hPrimaries_evtTrVtx = new TH2F(*((TH2F*) l_acorr->FindObject("fHistNonDetectableCorrNum"))); // all prim in sel events
-
- TH2F *hAlphaCorr = new TH2F("AlphaCorr","",nBinEtaCorr,binLimEtaCorr,nBinVtx,-MaxVtx,MaxVtx);
- hAlphaCorr->GetXaxis()->SetTitle("#eta");
- hAlphaCorr->GetYaxis()->SetTitle("z_{vtx} [cm]");
- hAlphaCorr->Add(hPrimaries_evtTrVtx);
- hAlphaCorr->Divide(hCombBkgCorr2MC);
- hAlphaCorr->Divide(hCombBkgMCDen);
new TCanvas();
- hAlphaCorr->Draw();
-
+ hPrimaries_evtTrVtx->Draw();
+
TH2F *hInvAlphaCorr = new TH2F("InvAlphaCorr","",nBinEtaCorr,binLimEtaCorr,nBinVtx,-MaxVtx,MaxVtx);
hInvAlphaCorr->Sumw2();
- hInvAlphaCorr->Multiply(hCombBkgCorr2MC,hCombBkgMCDen);
- hInvAlphaCorr->Divide(hPrimaries_evtTrVtx);
+
+ TH2F *hGenProtons = new TH2F(*((TH2F*) l_acorr->FindObject("fHistProtons")));
+ TH2F *hGenKaons = new TH2F(*((TH2F*) l_acorr->FindObject("fHistKaons")));
+ TH2F *hGenPions = new TH2F(*((TH2F*) l_acorr->FindObject("fHistPions")));
+ TH2F *hGenOthers = new TH2F(*((TH2F*) l_acorr->FindObject("fHistOthers")));
+
+ TH2F *hPrimReweighted = new TH2F("primReweighted","",nBinEtaCorr,binLimEtaCorr,nBinVtx,-MaxVtx,MaxVtx);
+ hPrimReweighted->Add(hGenKaons);
+ hPrimReweighted->Scale(2.25); //double kaon fraction
+ hPrimReweighted->Add(hGenProtons);
+ hPrimReweighted->Add(hGenPions);
+ hPrimReweighted->Add(hGenOthers); //use this in alpha if flag for syst true
+
+ TH2F *hRecProtons = new TH2F(*((TH2F*) l_acorr->FindObject("fHistReconstructedProtons")));
+ TH2F *hRecKaons = new TH2F(*((TH2F*) l_acorr->FindObject("fHistReconstructedKaons")));
+ TH2F *hRecPions = new TH2F(*((TH2F*) l_acorr->FindObject("fHistReconstructedPions")));
+ TH2F *hRecOthers = new TH2F(*((TH2F*) l_acorr->FindObject("fHistReconstructedOthers")));
+ TH2F *hRecSec = new TH2F(*((TH2F*) l_acorr->FindObject("fHistReconstructedSec")));
+
+ TH2F *hRecReweighted = new TH2F("recReweighted","",nBinEtaCorr,binLimEtaCorr,nBinVtx,-MaxVtx,MaxVtx);
+ hRecReweighted->Add(hRecKaons);
+ hRecReweighted->Scale(2.25);
+ hRecReweighted->Add(hRecProtons);
+ hRecReweighted->Add(hRecPions);
+ hRecReweighted->Add(hRecOthers); //use this in alpha if flag for syst true
+ hRecReweighted->Add(hRecSec);
+ hRecReweighted->Add(hBkgCombLabels);
+
+ // 1-beta in MC
+ TH2F *hCombBkgCorr2MC = new TH2F("backgroundCorr2MC","Combinatorial background correction",nBinEtaCorr,binLimEtaCorr,nBinVtx,-MaxVtx,MaxVtx);
+
+ if (changeStrangeness) {
+ if (bkgFromMCLabels) hCombBkgCorrMC->Divide(hBkgCombLabels,hRecReweighted,1.,1.);
+ else hCombBkgCorrMC->Divide(hCombBkgMC,hRecReweighted,1.,1.);
+
+ for (Int_t ieta=0; ieta<nBinEtaCorr; ieta++)
+ for (Int_t jz=0; jz<nBinVtx; jz++)
+ hCombBkgCorr2MC->SetBinContent(ieta+1,jz+1,1-hCombBkgCorrMC->GetBinContent(ieta+1,jz+1));
+
+ hInvAlphaCorr->Multiply(hCombBkgCorr2MC,hRecReweighted);
+ hInvAlphaCorr->Divide(hPrimReweighted);
+ } else {
+
+ if (bkgFromMCLabels) hCombBkgCorrMC->Divide(hBkgCombLabels,hCombBkgMCDen,1.,1.);
+ else hCombBkgCorrMC->Divide(hCombBkgMC,hCombBkgMCDen,1.,1.);
+
+ for (Int_t ieta=0; ieta<nBinEtaCorr; ieta++)
+ for (Int_t jz=0; jz<nBinVtx; jz++)
+ hCombBkgCorr2MC->SetBinContent(ieta+1,jz+1,1-hCombBkgCorrMC->GetBinContent(ieta+1,jz+1));
+
+ hInvAlphaCorr->Multiply(hCombBkgCorr2MC,hCombBkgMCDen);
+ hInvAlphaCorr->Divide(hPrimaries_evtTrVtx);
+
+ }
+
+// new TCanvas();
+// hCombBkgCorr2MC->Draw();
+
+
+ // Efficiency for particle species
+ TH2F *hEffProtons = new TH2F("effProtons","",nBinEtaCorr,binLimEtaCorr,nBinVtx,-MaxVtx,MaxVtx);
+ TH2F *hEffKaons = new TH2F("effKaons","",nBinEtaCorr,binLimEtaCorr,nBinVtx,-MaxVtx,MaxVtx);
+ TH2F *hEffPions = new TH2F("effPions","",nBinEtaCorr,binLimEtaCorr,nBinVtx,-MaxVtx,MaxVtx);
+ TH2F *hEffOthers = new TH2F("effOthers","",nBinEtaCorr,binLimEtaCorr,nBinVtx,-MaxVtx,MaxVtx);
+
+ hEffProtons->Divide(hRecProtons,hGenProtons,1.,1.);
+ hEffKaons->Divide(hRecKaons,hGenKaons,1.,1.);
+ hEffPions->Divide(hRecPions,hGenPions,1.,1.);
+ hEffOthers->Divide(hRecOthers,hGenOthers,1.,1.);
+
+// hEffKaons->Divide(hEffPions);
+
+
new TCanvas();
- hInvAlphaCorr->Draw();
- // Errors on alpha
+ hInvAlphaCorr->Draw();
+
+ // Errors on alpha
for (Int_t ieta=0; ieta<nBinEtaCorr; ieta++) {
for (Int_t jz=0; jz<nBinVtx; jz++) {
}
}
- // Number of events and normalization histogram
+ TH2F *hAlphaCorr = new TH2F("AlphaCorr","",nBinEtaCorr,binLimEtaCorr,nBinVtx,-MaxVtx,MaxVtx);
+ hAlphaCorr->GetXaxis()->SetTitle("#eta");
+ hAlphaCorr->GetYaxis()->SetTitle("z_{vtx} [cm]");
+
+ for (Int_t ieta=0; ieta<nBinEtaCorr; ieta++) {
+ for (Int_t jz=0; jz<nBinVtx; jz++) {
+ hAlphaCorr->SetBinContent(ieta+1,jz+1,1./hInvAlphaCorr->GetBinContent(ieta+1,jz+1));
+ }
+ }
+
+ new TCanvas();
+ hAlphaCorr->Draw();
+
+
+ //////////////////////// Factorize alpha ///////////////////////////
+
+ ////////////////////////////////////////////////////////////////////
+
+ // Number of events and normalization histograms
TH1D *hSPDvertex = (TH1D*)l_raw->FindObject("fHistSPDvtxAnalysis");
Double_t nEvtsRec=hSPDvertex->Integral(binVtxStart+1,binVtxStop);
cout<<"Number of reconstructed events in the selected vertex range: "<<nEvtsRec <<endl;
new TCanvas();
hCombBkgCorrData->Draw();
-/////////recheck!!!!!!!!
// MC distribution
//... to compare corrected distribution to MC in selected events
TH2F *hVertexMC2D = (TH2F*)l_MC->FindObject("fHistSelEvts");
- TH1D *hMCvertex = new TH1D(*(hVertexMC2D->ProjectionY("MCvertex",0,-1,"e"))); //MC vertex distrib
- TH2F *Eta = (TH2F*)l_MC->FindObject("fHistNonDetectableCorrNum");
+ TH1D *hMCvertex = new TH1D(*(hVertexMC2D->ProjectionY("MCvertex",0,-1,"e"))); //MC vertex distrib
+ new TCanvas();
+ hMCvertex->Draw();
-// TH2F *Eta = (TH2F*)l_MC->FindObject("fHistAllPrimaries"); // to compare with MC selection!!
+ TH2F *Eta = (TH2F*)l_MC->FindObject("fHistNonDetectableCorrNum");
TH1D *hdNdEta = new TH1D("dNdEta","Pseudorapidity ",nBinEtaCorr,binLimEtaCorr);
- hdNdEta->Sumw2();
- hdNdEta = Eta->ProjectionX("dNdEta",0,-1,"e");
+ hdNdEta = Eta->ProjectionX("dNdEta",0,-1,"e"); //here already all events
Double_t nEvtsTot=hMCvertex->Integral(0,hMCvertex->GetNbinsX()+1);
+
Double_t nEvtsTotSel= hMCvertex->Integral(binVtxStart+1,binVtxStop);
- cout<<"Number of generated events "<<nEvtsTot<<endl;
+ cout<<"Number of generated events: "<<nEvtsTot<<endl;
cout<<"Number of generated events in the selected vertex range: "<<nEvtsTotSel <<endl;
- Float_t nprimcentraleta = Eta->ProjectionX("dNdEta",0,-1,"e")->Integral(26,35);
- hdNdEta->Scale(nBinsPerPseudorapidityUnit/nEvtsTot);
+ Float_t nprimcentraleta = Eta->ProjectionX("dNdEta",0,-1,"e")->Integral(26,35); // project all or only in the sel vertex range?
+ hdNdEta->Scale(nBinsPerPseudorapidityUnit/nEvtsTot); // if so change number of events
hdNdEta->GetXaxis()->SetTitle("#eta");
hdNdEta->GetYaxis()->SetTitle("dN_{ch}/d#eta");
-
+ for (Int_t i=0;i<nBinEtaCorr;i++) hdNdEta->SetBinError(i+1,0.);
// dNdEta in +- 0.5
- cout<<"dN/dEta in |eta|<0.5"<<nprimcentraleta/nEvtsTot<<endl;
-
+ cout<<"Monte Carlo dN/dEta in |eta|<0.5 (selected events for analysis and all events!): "<<nprimcentraleta/nEvtsTot<<endl;
+ cout<<"Monte Carlo dN/dEta in |eta|<0.5 (selected events for analysis and events in vtx range!): "<<(Eta->ProjectionX("_x",binVtxStart+1,binVtxStop,"e")->Integral(26,35))/nEvtsTotSel<<endl;
new TCanvas();
hdNdEta->Draw();
-/////////////////////////
+
//Corrected distributions
TH2F *hSPDEta_2D_bkgCorr = new TH2F("SPDEta_2D_bkgCorr", "",nBinEtaCorr,binLimEtaCorr,nBinVtx,-MaxVtx,MaxVtx);
hSPDEta_2D_fullyCorr->Sumw2();
TH1F *hnorm_fullyCorr = new TH1F("norm_fullyCorr", "",nBinEtaCorr,binLimEtaCorr);
- hnorm_fullyCorr->Sumw2();
-
//dN/deta
TH1F *hdNdEta_raw = new TH1F("dNdEta_raw","",nBinEtaCorr,binLimEtaCorr);
//Apply corrections at
//...track level
- hSPDEta_2D_bkgCorr->Add(hSPDEta_2Draw,-1);
- hSPDEta_2D_bkgCorr->Multiply(hCombBkgCorrData);
hSPDEta_2D_bkgCorr->Add(hSPDEta_2Draw);
+ hSPDEta_2D_bkgCorr->Multiply(hCombBkgCorr2Data);
hSPDEta_2D_fullyCorr->Add(hSPDEta_2D_bkgCorr);
hSPDEta_2D_fullyCorr->Divide(hInvAlphaCorr);
new TCanvas();
hdNdEta_bkgCorr->Scale(1./nEvtsRec);
hdNdEta_fullyCorr->Divide(hSPDEta_2D_fullyCorr->ProjectionX("SPDEta_2D_fullyCorr_x",binVtxStart+1,binVtxStop,"e"),hnorm_fullyCorr);
+ TH1F *hCorrecteddNdEtaCentr = new TH1F("correcteddNdEtaCentr","",10,-0.5,0.5);
+ hCorrecteddNdEtaCentr->Sumw2();
+ for (Int_t jeta=0; jeta<10; jeta++) {
+ hCorrecteddNdEtaCentr->SetBinContent(jeta+1,hdNdEta_fullyCorr->GetBinContent(26+jeta));
+ hCorrecteddNdEtaCentr->SetBinError(jeta+1,hdNdEta_fullyCorr->GetBinError(26+jeta));
+ }
+ hCorrecteddNdEtaCentr->Rebin(10);
+
+ new TCanvas();
+ hCorrecteddNdEtaCentr->Draw();
+
Float_t ncorrtrackletscentr = hSPDEta_2D_fullyCorr->ProjectionX("SPDEta_2D_fullyCorr_x",binVtxStart+1,binVtxStop,"e")->Integral(26,35);
Float_t nevtscentr = hnorm_fullyCorr->Integral(26,35)/10.;
- cout<<"Corrected dN/dEta in |eta|<0.5 "<<ncorrtrackletscentr/nevtscentr<<endl;
+ cout<<""<<endl;
+ cout<<"Corrected dN/dEta in |eta|<0.5: "<<ncorrtrackletscentr/nevtscentr<<endl;
+ cout<<"Error on corrected tracklets in |eta|<0.5: +-"<<hCorrecteddNdEtaCentr->GetBinError(1)<<endl;
+ // dN/dEta per participant pairs
+
+ Float_t NpartV0Glauber = 381.188;
+ Float_t NpartCl2Glauber = 378.5;
+ Float_t NpartCl2Hij = 365.3;
+
+ Float_t errNpartV0Glauber = 18.4951;
+ Float_t errNpartCl2Glauber = 7.57; // 2% not used
+ Float_t errNpartCl2Hij = 7.306;
+
+ Float_t dNdEtaCorr = ncorrtrackletscentr/nevtscentr;
+ Float_t errdNdEtaCorr = hCorrecteddNdEtaCentr->GetBinError(1);
+
+ Float_t dNdEtaPartPairsV0G = dNdEtaCorr/(.5*NpartV0Glauber);
+ Float_t dNdEtaPartPairsCl2G = dNdEtaCorr/(.5*NpartCl2Glauber);
+ Float_t dNdEtaPartPairsCl2H = dNdEtaCorr/(.5*NpartCl2Hij);
+
+ Float_t errdNdEtaPartPairsV0G = CalculateErrordNdEtaPerPartPair(dNdEtaCorr,NpartV0Glauber,errdNdEtaCorr,errNpartV0Glauber);
+ Float_t errdNdEtaPartPairsCl2G = CalculateErrordNdEtaPerPartPair(dNdEtaCorr,NpartCl2Glauber,errdNdEtaCorr,errNpartCl2Glauber);
+ Float_t errdNdEtaPartPairsCl2H = CalculateErrordNdEtaPerPartPair(dNdEtaCorr,NpartCl2Hij,errdNdEtaCorr,errNpartCl2Hij);
+
+ cout<<"V0 Glauber dN/dEta/.5<Npart> "<<dNdEtaPartPairsV0G<<" +- "<<errdNdEtaPartPairsV0G<<endl;
+ cout<<"Cl2 Glauber dN/dEta/.5<Npart> "<<dNdEtaPartPairsCl2G<<" +- "<<errdNdEtaPartPairsCl2G<<endl;
+ cout<<"Cl2 HIJING dN/dEta/.5<Npart> "<<dNdEtaPartPairsCl2H<<" +- "<<errdNdEtaPartPairsCl2H<<endl;
hdNdEta_bkgCorr->Scale(nBinsPerPseudorapidityUnit);
hdNdEta_fullyCorr->Scale(nBinsPerPseudorapidityUnit);
- new TCanvas();
- hdNdEta_bkgCorr->Draw();
- new TCanvas();
- hdNdEta_fullyCorr->Draw();
- hdNdEta->Draw("same");
// Create mask for MC prim in SPD acceptance
TH1D *hdNdEta_test = new TH1D("dNdEta_test","Pseudorapidity ",nBinEtaCorr,binLimEtaCorr);
TH2F *hMask = new TH2F("Mask","",nBinEtaCorr,binLimEtaCorr,nBinVtx,-MaxVtx,MaxVtx);
- hMask->Sumw2();
for (Int_t jeta=0; jeta<nBinEtaCorr; jeta++) {
for (Int_t kvtx=0; kvtx<nBinVtx; kvtx++) {
if (hSPDEta_2D_fullyCorr->GetBinContent(jeta+1,kvtx+1)!=0) hMask->SetBinContent(jeta+1,kvtx+1,1.);
hdNdEta_test->Divide(hnorm_fullyCorr);
// hdNdEta_test->Scale(1./nEvtsTotSel); // -> number of MC events
hdNdEta_test->Scale(nBinsPerPseudorapidityUnit);
+ for (Int_t i=0;i<nBinEtaCorr;i++) hdNdEta_test->SetBinError(i+1,0.);
TH1F* hRatiotest=new TH1F("ratiotest","",nBinEtaCorr,binLimEtaCorr);
hRatiotest->GetXaxis()->SetTitle("#eta");
hRatiotest->GetYaxis()->SetTitle("Generated/Corrected");
-
- for (Int_t i=0;i<nBinEtaCorr;i++) hdNdEta_test->SetBinError(i+1,0.);
hRatiotest->Divide(hdNdEta_test,hdNdEta_fullyCorr,1.,1.);
new TCanvas();
hRatiotest->Draw("p");
-
- // Generated/corrected ratios
+ // Generated/corrected ratios for consistency checks
TH2F* hRatiodNdEta2D=new TH2F("ratiodNdEta2D","",nBinEtaCorr,binLimEtaCorr,nBinVtx,-MaxVtx,MaxVtx);
hRatiodNdEta2D->Divide(hSPDEta_2D_fullyCorr,Eta,1.,1.);
new TCanvas();
hRatiodNdEta2D->Draw();
- TH1F* hRatiodNdEta=new TH1F("ratiodNdEta","",nBinEtaCorr,binLimEtaCorr); // Add ratio for two corrected dN/dEta
+ TH1F* hRatiodNdEta=new TH1F("ratiodNdEta","",nBinEtaCorr,binLimEtaCorr);
hRatiodNdEta->GetXaxis()->SetTitle("#eta");
hRatiodNdEta->GetYaxis()->SetTitle("Generated/Corrected");
-
- for (Int_t i=0;i<nBinEtaCorr;i++) hdNdEta->SetBinError(i+1,0.);
hRatiodNdEta->Divide(hdNdEta,hdNdEta_fullyCorr,1.,1.);
new TCanvas();
hRatiodNdEta->Draw("p");
- // Draw dN/dEta
+ // Draw dN/dEta
new TCanvas();
-/* hdNdEta->SetLineWidth(3);
- hdNdEta->SetMinimum(0.);
- hdNdEta->SetMaximum(4000);
- hdNdEta->Draw("histo");
-*/
+// hdNdEta->SetLineWidth(3);
+// hdNdEta->Draw("histo");
+
+ hdNdEta_raw->SetMinimum(0.);
+ hdNdEta_raw->SetMaximum(3000);
hdNdEta_raw->SetLineColor(kGreen);
hdNdEta_raw->SetLineWidth(3);
hdNdEta_raw->Draw("histo");
leg1->SetTextSize(0.0304569);
TLegendEntry *entry1=leg1->AddEntry(hdNdEta_raw,"Raw","l");
// entry1=leg1->AddEntry(hdNdEta,"Generated","l");
- entry1=leg1->AddEntry(hdNdEta_bkgCorr,"Bkg corrected","p");
- entry1=leg1->AddEntry(hdNdEta_fullyCorr,"Corrected","p");
+ entry1=leg1->AddEntry(hdNdEta_bkgCorr,"Comb bkg corrected","p");
+ entry1=leg1->AddEntry(hdNdEta_fullyCorr,"Eff/acc corrected","p");
leg1->Draw();
+
/* new TCanvas();
// plot the relative stat error for this correction
TH2F* hStatErrTrackToPart = new TH2F("staterrperctracktopart","",nBinEtaCorr,binLimEtaCorr,nBinVtx,-MaxVtx,MaxVtx);
hStatErrTrackToPart->Draw();
new TCanvas();
hRatiodNdEta_trackToParticle->Draw();
-new TCanvas();
-hTrackToParticleCorr->Draw();
+ new TCanvas();
+ hTrackToParticleCorr->Draw();*/
+
// Write histos
-*/
TFile *fout= new TFile("SPDdNdEta.root","RECREATE");
-
hCombBkgCorrData->Write();
hCombBkgCorrMC->Write();
hCombBkgCorr2MC->Write();
hAlphaCorr->Write();
+ hdNdEta_raw->Write();
+ hdNdEta_bkgCorr->Write();
hdNdEta_fullyCorr->Write();
-
+
+ hCorrecteddNdEtaCentr->Write();
+
+ hEffProtons->Write();
+ hEffKaons->Write();
+ hEffPions->Write();
+ hEffOthers->Write();
+
fout->Write();
fout->Close();
}
+
+Float_t CalculateErrordNdEtaPerPartPair(Float_t a, Float_t b, Float_t da, Float_t db) {
+
+ Float_t errdndetaperpartpairs = (a/b)*TMath::Sqrt(TMath::Power(da,2)/TMath::Power(a,2)+TMath::Power(db,2)/TMath::Power(b,2));
+ return errdndetaperpartpairs;
+
+}
+
+void combscalingfactors (TList *lall, TList *lcomb, TList *lmc, Double_t *scaleNormRange_rot, Double_t *scaleNormRange_labels) {
+
+
+ TH2F *hall = (TH2F*) lall->FindObject("fHistSPDdePhideTheta");
+ TH2F *hcomb = (TH2F*) lcomb->FindObject("fHistSPDdePhideTheta");
+
+ TH2F *hlabel = (TH2F*) lmc->FindObject("fHistdPhidThetaComb");
+
+ TH1D* hproall = new TH1D("_xall","",2000,-1.,1.);
+ hproall->Sumw2();
+ hproall = hall->ProjectionX("_xa",0,-1,"e");
+
+ TH1D* hprocomb = new TH1D("_xcomb","",2000,-1.,1.);
+ hprocomb->Sumw2();
+ hprocomb = hcomb->ProjectionX("_xc",0,-1,"e");
+
+ TH1D *hprolabel = new TH1D("_xcomblabel","",2000,-1.,1.);
+ hprolabel = hlabel->ProjectionX("_xcombl",0,-1,"e");
+
+ // Normalize to integral in [0.08,0.20]
+ Double_t denNormRange = hproall->Integral(801,920) + hproall->Integral(1081,1200);
+ Double_t numNormRange_rot = hprocomb->Integral(801,920) + hprocomb->Integral(1081,1200);
+ *scaleNormRange_rot = denNormRange / numNormRange_rot;
+
+ new TCanvas();
+ hproall->Draw("histo");
+
+ hprocomb->Scale(*scaleNormRange_rot);
+ cout<<"Scaling factor normalizing to close tails: "<<*scaleNormRange_rot<<endl;
+ cout<<"Percentage of background"<< hprocomb->Integral(1,2000)/hproall->Integral(1,2000)<<endl;
+ cout<<"Percentage of background in |#D#phi|<0.08 "<< hprocomb->Integral(921,1080)/hproall->Integral(921,1080)<<endl;
+ cout<<"Percentage of background in |#D#phi|<0.06 "<< hprocomb->Integral(941,1060)/hproall->Integral(941,1060)<<endl;
+ hprocomb->SetLineColor(kBlue);
+ hprocomb->Draw("histo,same");
+
+ // Fit the label bkg comb distribution to the data distribution
+ // Normalize to integral in [0.08,0.20]
+ Double_t numNormRange_labels = hprolabel->Integral(801,920) + hprolabel->Integral(1081,1200);
+ *scaleNormRange_labels = denNormRange / numNormRange_labels;
+
+ hprolabel->Scale(*scaleNormRange_labels);
+ cout<<"Scaling factor labels: "<<*scaleNormRange_labels<<endl;
+ cout<<"Percentage of background with labels "<< (hprolabel->Integral(1,2000)/hproall->Integral(1,2000))<<endl;
+ cout<<"Percentage of background with labels in |#D#phi|<0.08 "<< (hprolabel->Integral(921,1080)/hproall->Integral(921,1080))<<endl;
+
+ hprolabel->SetLineColor(kGreen);
+ hproall->Draw("histo,same");
+/*
+ TFile *fout= new TFile("scaling.root","RECREATE");
+ hproall->Write();
+ hprocomb->Write();
+ hprolabel->Write();
+ fout->Write();
+ fout->Close();
+*/
+}
+
+
-void runProofSPDdNdEta(TString proofCluster = "delia@alice-caf.cern.ch",//skaf.saske.sk",
- TString alirootVer = "VO_ALICE@AliRoot::v4-21-02-AN",
- TString rootVer = "VO_ALICE@ROOT::v5-27-06a-1",
-// TString dataset = "/alice/data/LHC10h_000137161_p1_plusplus",
- TString dataset = "/alice/data/LHC10h_000137161_p1_plus",
-// TString dataset = "/alice/data/LHC10h_000137045_p1_plus",
-// TString dataset = "/alice/sim/LHC10h3_000137161",
-// TString dataset = "/alice/sim/LHC10h2_000137161",
-// TString dataset = "/alice/sim/LHC10h1_000137045",
+void runProofSPDdNdEta_CorrRef4(
+ TString proofCluster = "mnicassi@alice-caf.cern.ch", //skaf.saske.sk",
+ TString alirootVer = "VO_ALICE@AliRoot::v4-21-04-AN",
+ TString rootVer = "VO_ALICE@ROOT::v5-27-06a-1",
+// TString dataset = "/alice/data/LHC10h_000137161_p1_plusplusplus",
+ TString dataset = "/alice/sim/LHC10h8_000137161",
TString outFileCorr = "SPDdNdEtaCorr.root",
+// TString outFileCorr = "SPDdNdEtaCorrRot.root",
TString outFileData = "SPDdNdEtaData.root",
TString percentFile = "./AliCentralityBy1D_LHC10g2a_100.root",
TString percentFile2 = "./AliCentralityByFunction_LHC10g2a_100.root",
-// Bool_t useMC = kTRUE,
- Bool_t useMC = kFALSE,
+ Bool_t useMC = kTRUE,
Bool_t readtr = kFALSE,
Bool_t recotracklets = kTRUE,
Bool_t dataonalien = kFALSE,
Float_t centrlowlim = 0.,
Float_t centruplim = 5.,
- TString centrest = "",
- Int_t minClMultLay2 = 4500,
+ Bool_t centrest = kFALSE,
+ Int_t minClMultLay2 = 4300,
+// Int_t minClMultLay2 = -1,
+ Int_t maxClMultLay2 = 1.0*1e5,
+ Int_t minV0Mult = 14674.5,
+ Float_t vtxlim = 7.,
+ Bool_t partsp = kTRUE,
Float_t phiWindow = 0.8,
Float_t thetaWindow = 0.025,
Float_t phiShift = 0.0045,
Float_t removeClFOvl = kFALSE,
Float_t phiOvlCut = 0.005,
Float_t zetaOvlCut = 0.05,
-// Float_t phiRotAngle = 0.,
- Float_t phiRotAngle = TMath::Pi(),
+ Float_t phiRotAngle = 0.,
+// Float_t phiRotAngle = TMath::Pi(),
Float_t phiWindowAna = 0.08,
Int_t nEvents = 1.0*1e7,
Int_t nEventsSkip = 0) { //1.0*1e7
// gROOT->LoadMacro("AnalysisMacroa.C");
Analysis(dataset.Data(), outFileCorr, outFileData, percentFile, percentFile2,
useMC, readtr, recotracklets,
- nEvents, nEventsSkip, centrlowlim, centruplim, centrest, minClMultLay2,
+ nEvents, nEventsSkip, centrlowlim, centruplim, centrest, minClMultLay2, maxClMultLay2, minV0Mult, vtxlim, partsp,
phiWindow, thetaWindow, phiShift, removeClFOvl, phiOvlCut, zetaOvlCut, phiRotAngle,phiWindowAna);
}
void Analysis(TString dataset, TString outFileCorr, TString outFileData, TString percentFile, TString percentFile2,
Bool_t useMC, Bool_t readtr, Bool_t recotracklets,
Int_t nEvents, Int_t nEventsSkip,
- Float_t centrlowlim, Float_t centruplim, TString centrest, Int_t minClMultLay2,
+ Float_t centrlowlim, Float_t centruplim, Bool_t centrest, Int_t minClMultLay2, Int_t maxClMultLay2, Int_t minV0Mult, Float_t vtxlim, Bool_t partsp,
Float_t phiWindow, Float_t thetaWindow, Float_t phiShift, Bool_t removeClFOvl,
Float_t phiOvlCut, Float_t zetaOvlCut, Float_t phiRotAngle, Float_t phiWindowAna) {
task->SetCentralityUpLim(centruplim);
task->SetCentralityEst(centrest);
task->SetMinClusterMultLay2(minClMultLay2);
+ task->SetMaxClusterMultLay2(maxClMultLay2);
+ task->SetMinV0Mult(minV0Mult);
+ task->SetVertexRange(vtxlim);
+ task->SetPartSpecies(partsp);
task->SetPhiWindow(phiWindow);
task->SetThetaWindow(thetaWindow);
if (!useMC) {
AliPhysicsSelection * physSel = physSelTask->GetPhysicsSelection();
- if (dataset=="/alice/data/LHC10h_000137045_p1_plus") {
+/* if (dataset=="/alice/data/LHC10h_000137045_p1_plus") {
physSel->AddCollisionTriggerClass("+CMBS1A-B-NOPF-ALL");
physSel->AddCollisionTriggerClass("+CMBS1C-B-NOPF-ALL");
physSel->AddCollisionTriggerClass("+CMBAC-B-NOPF-ALL");
- } else {
+ } else {*/
physSel->AddCollisionTriggerClass("+CMBS2A-B-NOPF-ALL");
physSel->AddCollisionTriggerClass("+CMBS2C-B-NOPF-ALL");
physSel->AddCollisionTriggerClass("+CMBAC-B-NOPF-ALL");
- }
+// }
task->SelectCollisionCandidates(AliVEvent::kUserDefined);
} else {
/alice/sim/LHC10h1a_000137045 | 1022 | /esdTree | 1.006e+04| 124 GB | 98 %
*/
-//________________________________________________________________________
+//________________________________________________________________________
\ No newline at end of file
git://git.uio.no / u / mrichter / AliRoot.git / commitdiff