From: loizides Date: Tue, 30 Nov 2010 22:43:32 +0000 (+0000) Subject: The code we used for the analysis of the first PbPb mult paper X-Git-Url: http://git.uio.no/git/?a=commitdiff_plain;h=ccfafc804d37e049de870bb4b63498bc0ac69846;p=u%2Fmrichter%2FAliRoot.git The code we used for the analysis of the first PbPb mult paper (rotation and MC labels + systematics): the task, the macro to apply corrections and calculate the scaling factors for combinatorial (now all in one macro) and the macro to run the task with all the parameters used for the first paper set by default. (From Maria) --- diff --git a/PWG2/EVCHAR/AliAnalysisTaskSPDdNdEta.cxx b/PWG2/EVCHAR/AliAnalysisTaskSPDdNdEta.cxx index df3e61e35fd..f2b56b8ca51 100644 --- a/PWG2/EVCHAR/AliAnalysisTaskSPDdNdEta.cxx +++ b/PWG2/EVCHAR/AliAnalysisTaskSPDdNdEta.cxx @@ -37,6 +37,7 @@ #include "AliESDInputHandler.h" #include "AliESDInputHandlerRP.h" +#include "AliESDVZERO.h" #include "../ITS/AliITSRecPoint.h" #include "AliCDBPath.h" #include "AliCDBManager.h" @@ -50,6 +51,7 @@ #include "AliTrackReference.h" #include "AliGenHijingEventHeader.h" +#include "AliGenDPMjetEventHeader.h" #include "AliLog.h" @@ -76,8 +78,12 @@ AliAnalysisTaskSPDdNdEta::AliAnalysisTaskSPDdNdEta(const char *name) fMCCentralityBin(AliAnalysisTaskSPDdNdEta::kall), fCentrLowLim(0), fCentrUpLim(0), - fCentrEst(""), + fCentrEst(kFALSE), fMinClMultLay2(0), + fMaxClMultLay2(0), + fMinMultV0(0), + fVtxLim(0), + fPartSpecies(0), fPhiWindow(0), fThetaWindow(0), @@ -90,6 +96,7 @@ AliAnalysisTaskSPDdNdEta::AliAnalysisTaskSPDdNdEta(const char *name) fMultReco(0), + fV0Ampl(0), fHistSPDRAWMultvsZ(0), fHistSPDRAWMultvsZTriggCentrEvts(0), fHistSPDRAWMultvsZCentrEvts(0), @@ -117,12 +124,22 @@ AliAnalysisTaskSPDdNdEta::AliAnalysisTaskSPDdNdEta(const char *name) 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), @@ -144,9 +161,9 @@ AliAnalysisTaskSPDdNdEta::AliAnalysisTaskSPDdNdEta(const char *name) fHistRecvsGenImpactPar(0), fHistMCNpart(0), - fHistdPhiPP(0), - fHistdPhiSS(0), - fHistdPhiComb(0), + fHistdPhidThetaPP(0), + fHistdPhidThetaSS(0), + fHistdPhidThetaComb(0), fHistDeVtx(0) @@ -185,7 +202,7 @@ void AliAnalysisTaskSPDdNdEta::UserCreateOutputObjects() 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"); @@ -198,7 +215,7 @@ void AliAnalysisTaskSPDdNdEta::UserCreateOutputObjects() fOutput = new TList(); fOutput->SetOwner(); - Int_t nBinVtx = 20; + Int_t nBinVtx = 40; Double_t MaxVtx = 20.; Int_t nBinMultCorr = 200; @@ -220,6 +237,10 @@ void AliAnalysisTaskSPDdNdEta::UserCreateOutputObjects() // 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]"); @@ -332,10 +353,22 @@ void AliAnalysisTaskSPDdNdEta::UserCreateOutputObjects() // 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); @@ -350,7 +383,19 @@ void AliAnalysisTaskSPDdNdEta::UserCreateOutputObjects() 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); @@ -414,12 +459,12 @@ void AliAnalysisTaskSPDdNdEta::UserCreateOutputObjects() 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]"); @@ -451,6 +496,7 @@ void AliAnalysisTaskSPDdNdEta::UserExec(Option_t *) 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; @@ -462,8 +508,7 @@ void AliAnalysisTaskSPDdNdEta::UserExec(Option_t *) // Centrality selection Bool_t eventInCentralityBin = kFALSE; - // Centrality selection - AliESDCentrality *centrality = fmyESD->GetCentrality(); +/* AliESDCentrality *centrality = fmyESD->GetCentrality(); if (fCentrEst=="") eventInCentralityBin = kTRUE; else { if(!centrality) { @@ -472,22 +517,38 @@ void AliAnalysisTaskSPDdNdEta::UserExec(Option_t *) 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])GetMultiplicity(); Int_t multSPD = multESD->GetNumberOfTracklets(); Int_t nSingleCl1 = multESD->GetNumberOfSingleClusters(); Int_t multSPDcl1 = nSingleCl1 + multSPD; @@ -499,9 +560,10 @@ void AliAnalysisTaskSPDdNdEta::UserExec(Option_t *) 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&&multSPDcl2Fill(multSPDcl1); fHistSPDmultcl2->Fill(multSPDcl2); fHistSPDmultcl1vscl2->Fill(multSPDcl1,multSPDcl2); @@ -610,12 +672,15 @@ void AliAnalysisTaskSPDdNdEta::UserExec(Option_t *) return; } + Float_t impactParameter = 0.; + Int_t npart = 0; + AliGenHijingEventHeader* hijingHeader = dynamic_cast(fMCEvent->Header()->GenEventHeader()); - if (!hijingHeader) { - Printf("Unknown header type. \n"); - return ; - } - Float_t impactParameter = hijingHeader->ImpactParameter(); + AliGenDPMjetEventHeader* dpmHeader = dynamic_cast(fMCEvent->Header()->GenEventHeader()); + + if (hijingHeader) impactParameter = hijingHeader->ImpactParameter(); + else if (dpmHeader) impactParameter = dpmHeader->ImpactParameter(); + Bool_t IsEventInMCCentralityBin = kFALSE; switch (fMCCentralityBin) { @@ -643,8 +708,11 @@ void AliAnalysisTaskSPDdNdEta::UserExec(Option_t *) 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(); @@ -656,14 +724,14 @@ void AliAnalysisTaskSPDdNdEta::UserExec(Option_t *) // 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) { @@ -674,8 +742,9 @@ void AliAnalysisTaskSPDdNdEta::UserExec(Option_t *) // Loop over MC particles for (Int_t imc=0; imcGetTrack(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; @@ -688,7 +757,6 @@ void AliAnalysisTaskSPDdNdEta::UserExec(Option_t *) detectedPrimaryPartLay1[multMCCharged]=kFALSE; detectedPrimaryPartLay2[multMCCharged]=kFALSE; detectablePrimaryPart[multMCCharged]=kFALSE; - primCounted[multMCCharged]=kFALSE; if (fTR) { Int_t nref = mcpart->GetNumberOfTrackReferences(); @@ -704,6 +772,12 @@ void AliAnalysisTaskSPDdNdEta::UserExec(Option_t *) } } } + 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 @@ -711,39 +785,61 @@ void AliAnalysisTaskSPDdNdEta::UserExec(Option_t *) fHistMCmultEtacut->Fill(multMCChargedEtacut); // Event selection: in centrality bin, triggered with vertex - if (eventTriggered&&eventWithVertex&&eventInCentralityBin&&multSPDcl2>fMinClMultLay2) { + if (eventTriggered&&eventWithVertex&&eventInCentralityBin&&multSPDcl2>fMinClMultLay2&&multSPDcl2Fill(vtxMC[2],vtxMC[2]-esdvtx[2]); for (Int_t itracklet=0; itrackletFill(trackletCoord[itracklet][4],esdvtx[2]); - if (trackletLab[itracklet][0]==trackletLab[itracklet][1]) { - Bool_t trakletByPrim = kFALSE; - for (Int_t imc=0; imcFill(trackletCoord[itracklet][0]); + if (trackletLab[itracklet][0]==trackletLab[itracklet][1]) { + Bool_t trakletByPrim = kFALSE; + for (Int_t imc=0; imcFill(trackletCoord[itracklet][0],trackletCoord[itracklet][1]); + if (TMath::Abs(trackletCoord[itracklet][0])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])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])Fill(trackletCoord[itracklet][4],esdvtx[2]); + fHistBkgCombLabels->Fill(trackletCoord[itracklet][4],esdvtx[2]); } - } // cut dphi + } } // end loop tracklets for (Int_t imc=0; imcGet("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: "<FindObject("fHistSPDRAWEtavsZ"))); hSPDEta_2Draw->SetNameTitle("SPDEta_2Draw","Reconstructed tracklets"); @@ -91,8 +109,7 @@ void correctSPDdNdEta (Char_t* fileRaw, Char_t* fileRawBkgCorr, Char_t* fileAlph 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); @@ -102,47 +119,119 @@ void correctSPDdNdEta (Char_t* fileRaw, Char_t* fileRawBkgCorr, Char_t* fileAlph 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; ietaSetBinContent(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; ietaSetBinContent(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; ietaSetBinError(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; ietaSetBinContent(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; ietaSetBinContent(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; ietaGetXaxis()->SetTitle("#eta"); + hAlphaCorr->GetYaxis()->SetTitle("z_{vtx} [cm]"); + + for (Int_t ieta=0; ietaSetBinContent(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: "<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 "<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;iSetBinError(i+1,0.); // dNdEta in +- 0.5 - cout<<"dN/dEta in |eta|<0.5"<ProjectionX("_x",binVtxStart+1,binVtxStop,"e")->Integral(26,35))/nEvtsTotSel<Draw(); -///////////////////////// + //Corrected distributions TH2F *hSPDEta_2D_bkgCorr = new TH2F("SPDEta_2D_bkgCorr", "",nBinEtaCorr,binLimEtaCorr,nBinVtx,-MaxVtx,MaxVtx); @@ -219,8 +326,6 @@ void correctSPDdNdEta (Char_t* fileRaw, Char_t* fileRawBkgCorr, Char_t* fileAlph 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); @@ -239,9 +344,8 @@ void correctSPDdNdEta (Char_t* fileRaw, Char_t* fileRawBkgCorr, Char_t* fileAlph //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(); @@ -269,22 +373,53 @@ void correctSPDdNdEta (Char_t* fileRaw, Char_t* fileRawBkgCorr, Char_t* fileAlph 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 "<GetBinError(1)<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 "< "< "<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; jetaGetBinContent(jeta+1,kvtx+1)!=0) hMask->SetBinContent(jeta+1,kvtx+1,1.); @@ -295,39 +430,35 @@ void correctSPDdNdEta (Char_t* fileRaw, Char_t* fileRawBkgCorr, Char_t* fileAlph hdNdEta_test->Divide(hnorm_fullyCorr); // hdNdEta_test->Scale(1./nEvtsTotSel); // -> number of MC events hdNdEta_test->Scale(nBinsPerPseudorapidityUnit); + for (Int_t i=0;iSetBinError(i+1,0.); TH1F* hRatiotest=new TH1F("ratiotest","",nBinEtaCorr,binLimEtaCorr); hRatiotest->GetXaxis()->SetTitle("#eta"); hRatiotest->GetYaxis()->SetTitle("Generated/Corrected"); - - for (Int_t i=0;iSetBinError(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;iSetBinError(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"); @@ -347,10 +478,11 @@ void correctSPDdNdEta (Char_t* fileRaw, Char_t* fileRawBkgCorr, Char_t* fileAlph 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); @@ -368,21 +500,95 @@ void correctSPDdNdEta (Char_t* fileRaw, Char_t* fileRawBkgCorr, Char_t* fileAlph 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<Integral(1,2000)/hproall->Integral(1,2000)<Integral(921,1080)/hproall->Integral(921,1080)<Integral(941,1060)/hproall->Integral(941,1060)<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<Integral(1,2000)/hproall->Integral(1,2000))<Integral(921,1080)/hproall->Integral(921,1080))<SetLineColor(kGreen); + hproall->Draw("histo,same"); +/* + TFile *fout= new TFile("scaling.root","RECREATE"); + hproall->Write(); + hprocomb->Write(); + hprolabel->Write(); + fout->Write(); + fout->Close(); +*/ +} + + diff --git a/PWG2/EVCHAR/runProofSPDdNdEta.C b/PWG2/EVCHAR/runProofSPDdNdEta.C index 3c85bdbe502..9fb3f010ca0 100644 --- a/PWG2/EVCHAR/runProofSPDdNdEta.C +++ b/PWG2/EVCHAR/runProofSPDdNdEta.C @@ -1,33 +1,35 @@ -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 @@ -64,7 +66,7 @@ void runProofSPDdNdEta(TString proofCluster = "delia@alice-caf.cern.ch",//skaf // 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); } @@ -73,7 +75,7 @@ void runProofSPDdNdEta(TString proofCluster = "delia@alice-caf.cern.ch",//skaf 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) { @@ -109,6 +111,10 @@ void Analysis(TString dataset, TString outFileCorr, TString outFileData, TString 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); @@ -125,16 +131,16 @@ void Analysis(TString dataset, TString outFileCorr, TString outFileData, TString 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 { @@ -273,4 +279,4 @@ Bool_t InputHandlerSetup(TString format = "esd", Bool_t useKine = kTRUE, Bool_t /alice/sim/LHC10h1a_000137045 | 1022 | /esdTree | 1.006e+04| 124 GB | 98 % */ -//________________________________________________________________________ +//________________________________________________________________________ \ No newline at end of file