From a1c5e4a1669d3392094284375d732f2bcf2a6c40 Mon Sep 17 00:00:00 2001 From: fprino Date: Thu, 7 Aug 2014 00:31:34 +0200 Subject: [PATCH] Updates in ITSsa pi,kp task: new multiplicity estimator, storeage of DCA cut, code cleanup. Added macro to fit the DCA distributions --- PWGLF/CMakelibPWGLFspectra.pkg | 1 + PWGLF/PWGLFspectraLinkDef.h | 1 + .../PiKaPr/ITSsa/AddTaskITSsaSpectra.C | 13 +- .../ITSsa/AliAnalysisTaskSEITSsaSpectra.cxx | 403 ++++++------------ .../ITSsa/AliAnalysisTaskSEITSsaSpectra.h | 19 +- .../PiKaPr/ITSsa/FitDCADistributions.C | 303 +++++++++++++ 6 files changed, 448 insertions(+), 292 deletions(-) create mode 100644 PWGLF/SPECTRA/PiKaPr/ITSsa/FitDCADistributions.C diff --git a/PWGLF/CMakelibPWGLFspectra.pkg b/PWGLF/CMakelibPWGLFspectra.pkg index e83dea4430e..93d03994cee 100644 --- a/PWGLF/CMakelibPWGLFspectra.pkg +++ b/PWGLF/CMakelibPWGLFspectra.pkg @@ -44,6 +44,7 @@ set ( SRCS SPECTRA/ChargedHadrons/dNdPt/AlidNdPtAcceptanceCuts.cxx SPECTRA/ChargedHadrons/dNdPt/AlidNdPtTrackDumpTask.cxx SPECTRA/ChargedHadrons/dNdPt/AliPtResolAnalysis.cxx SPECTRA/ChargedHadrons/dNdPt/AliPtResolAnalysisPbPb.cxx + SPECTRA/PiKaPr/ITSsa/AliAnalysisTaskSEITSsaSpectra.cxx SPECTRA/PiKaPr/TPCTOF/AliAnalysisCombinedHadronSpectra.cxx SPECTRA/PiKaPr/TPCTOFpA/AliAnalysisTPCTOFpA.cxx SPECTRA/PiKaPr/TOF/pp7/TOFSpectrappAnalysis.cxx diff --git a/PWGLF/PWGLFspectraLinkDef.h b/PWGLF/PWGLFspectraLinkDef.h index b65c16558e4..48a5235bc7f 100644 --- a/PWGLF/PWGLFspectraLinkDef.h +++ b/PWGLF/PWGLFspectraLinkDef.h @@ -18,6 +18,7 @@ #pragma link C++ class AlidNdPtTrackDumpTask+; #pragma link C++ class AliPtResolAnalysis+; #pragma link C++ class AliPtResolAnalysisPbPb+; +#pragma link C++ class AliAnalysisTaskSEITSsaSpectra+; #pragma link C++ class AliAnalysisCombinedHadronSpectra+; #pragma link C++ class AliAnalysisTPCTOFpA+; #pragma link C++ class TOFSpectrappAnalysis+; diff --git a/PWGLF/SPECTRA/PiKaPr/ITSsa/AddTaskITSsaSpectra.C b/PWGLF/SPECTRA/PiKaPr/ITSsa/AddTaskITSsaSpectra.C index 15eb06850b3..4d092595ed5 100644 --- a/PWGLF/SPECTRA/PiKaPr/ITSsa/AddTaskITSsaSpectra.C +++ b/PWGLF/SPECTRA/PiKaPr/ITSsa/AddTaskITSsaSpectra.C @@ -55,7 +55,7 @@ AliAnalysisTaskSEITSsaSpectra *AddTaskITSsaSpectra(Int_t optNtuple=0,Int_t readM outputFileName += ":PWG2SpectraITSsa"; AliAnalysisDataContainer *coutput =0x0; - + AliAnalysisDataContainer *coutputCuts = 0x0; if(hi) { coutput = mgr->CreateContainer(Form("clistITSsaCent%ito%i",lowcut,upcut), @@ -67,11 +67,18 @@ AliAnalysisTaskSEITSsaSpectra *AddTaskITSsaSpectra(Int_t optNtuple=0,Int_t readM { coutput = mgr->CreateContainer(Form("clistITSsaMult%ito%i",lowcut,upcut), TList::Class(), - AliAnalysisManager::kOutputContainer, - outputFileName ); + AliAnalysisManager::kOutputContainer, + outputFileName ); } + + coutputCuts = mgr->CreateContainer(Form("DCACutMult%ito%i",lowcut,upcut), + TList::Class(), + AliAnalysisManager::kParamContainer, + outputFileName ); + mgr->ConnectInput(taskits, 0, mgr->GetCommonInputContainer()); mgr->ConnectOutput(taskits, 1, coutput); + mgr->ConnectOutput(taskits, 2, coutputCuts); return taskits; } diff --git a/PWGLF/SPECTRA/PiKaPr/ITSsa/AliAnalysisTaskSEITSsaSpectra.cxx b/PWGLF/SPECTRA/PiKaPr/ITSsa/AliAnalysisTaskSEITSsaSpectra.cxx index def4a83e070..8b3ef29f3a4 100644 --- a/PWGLF/SPECTRA/PiKaPr/ITSsa/AliAnalysisTaskSEITSsaSpectra.cxx +++ b/PWGLF/SPECTRA/PiKaPr/ITSsa/AliAnalysisTaskSEITSsaSpectra.cxx @@ -25,6 +25,7 @@ /////////////////////////////////////////////////////////////////////////// #include +#include #include #include #include @@ -53,10 +54,10 @@ ClassImp(AliAnalysisTaskSEITSsaSpectra) //________________________________________________________________________ AliAnalysisTaskSEITSsaSpectra::AliAnalysisTaskSEITSsaSpectra(): -AliAnalysisTaskSE("Task CFit"), +AliAnalysisTaskSE("taskITSsaSpectra"), fESD(0), - fesdTrackCutsMult(0), fOutput(0), + fListCuts(0), fHistNEvents(0), fHistMult(0), fHistCen(0), @@ -67,6 +68,8 @@ AliAnalysisTaskSE("Task CFit"), fHistDEDXdouble(0), fHistBeforeEvSel(0), fHistAfterEvSel(0), + fDCAxyCutFunc(0x0), + fDCAzCutFunc(0x0), fITSPIDResponse(0), fMinSPDPts(1), fMinNdEdxSamples(3), @@ -81,7 +84,7 @@ AliAnalysisTaskSE("Task CFit"), fUpMult(-1), fLowCentrality(-1.0), fUpCentrality(-1.0), - fSPD(0), + fMultEstimator(0), fHImode(0), fYear(2010), fMC(kFALSE), @@ -98,25 +101,9 @@ AliAnalysisTaskSE("Task CFit"), Double_t xbins[kNbins+1]={0.08,0.10,0.12,0.14,0.16,0.18,0.20,0.25,0.30,0.35,0.40,0.45,0.50,0.55,0.60,0.65,0.70,0.75,0.80,0.85,0.90,0.95,1.0}; for(Int_t iBin=0; iBinSetMinNClustersTPC(70); - fesdTrackCutsMult->SetMaxChi2PerClusterTPC(4); - fesdTrackCutsMult->SetAcceptKinkDaughters(kFALSE); - fesdTrackCutsMult->SetRequireTPCRefit(kTRUE); - // ITS - fesdTrackCutsMult->SetRequireITSRefit(kTRUE); - fesdTrackCutsMult->SetClusterRequirementITS(AliESDtrackCuts::kSPD, - AliESDtrackCuts::kAny); - fesdTrackCutsMult->SetDCAToVertex2D(kFALSE); - fesdTrackCutsMult->SetRequireSigmaToVertex(kFALSE); - fesdTrackCutsMult->SetEtaRange(-0.8,+0.8); - fesdTrackCutsMult->SetPtRange(0.15, 1e10); - SetYear(2010); - DefineInput(0, TChain::Class()); DefineOutput(1, TList::Class()); + DefineOutput(2,TList::Class()); Printf("end of AliAnalysisTaskSEITSsaSpectra"); } @@ -127,8 +114,14 @@ AliAnalysisTaskSEITSsaSpectra::~AliAnalysisTaskSEITSsaSpectra(){ delete fOutput; fOutput = 0; } + if (fListCuts) { + delete fListCuts; + fListCuts = 0; + } if(fRandGener) delete fRandGener; if(fITSPIDResponse) delete fITSPIDResponse; + delete fDCAxyCutFunc; + delete fDCAzCutFunc; } //________________________________________________________________________ @@ -170,26 +163,37 @@ Double_t AliAnalysisTaskSEITSsaSpectra::CookdEdx(Double_t *s) const { //________________________________________________________________________ -void AliAnalysisTaskSEITSsaSpectra::SetYear(Int_t year){ - // Set year dependent quantities - fYear=year; - if(fYear==2009){ - fesdTrackCutsMult->SetMaxDCAToVertexXYPtDep("0.0350+0.0420/TMath::Power(pt,0.9)"); //2009 standard cut - fesdTrackCutsMult->SetMaxDCAToVertexZ(20); //2009 standard cut - }else{ - fesdTrackCutsMult->SetMaxDCAToVertexXYPtDep("0.0182+0.0350/pt^1.01"); //2010 standard cut - fesdTrackCutsMult->SetMaxDCAToVertexZ(2); //2010 standard cut - } -} - -//________________________________________________________________________ -Bool_t AliAnalysisTaskSEITSsaSpectra::DCAcut(Double_t impactXY, Double_t impactZ, Double_t pt, Bool_t optMC) const { +Bool_t AliAnalysisTaskSEITSsaSpectra::DCAcut(Double_t impactXY, Double_t impactZ, Double_t pt) const { // cut on transverse impact parameter updaated on 20-5-2010 // from the study of L. Milano, F. Prino on the ITS standalone tracks // using the common binning of the TPC tracks + + Double_t xyMax = fDCAxyCutFunc->Eval(pt); //in micron + if((TMath::Abs(impactXY)*10000)>xyMax) return kFALSE; + Double_t zMax = fDCAzCutFunc->Eval(pt); //in micron + if((TMath::Abs(impactZ)*10000)>zMax) return kFALSE; + + return kTRUE; +} + +//________________________________________________________________________ +Double_t AliAnalysisTaskSEITSsaSpectra::Eta2y(Double_t pt, Double_t m, Double_t eta) const { + // convert eta to y + Double_t mt = TMath::Sqrt(m*m + pt*pt); + return TMath::ASinH(pt/mt*TMath::SinH(eta)); +} + + +//________________________________________________________________________ +void AliAnalysisTaskSEITSsaSpectra::Init(){ + // + // Initialization + // + fListCuts=new TList(); + fListCuts->SetOwner(); Double_t xyP[3]; Double_t zP[3]; - if(optMC){ + if(fMC){ if(fYear==2009){ xyP[0]=88.63; //MC LHC10a12 xyP[1]=19.57; @@ -205,8 +209,7 @@ Bool_t AliAnalysisTaskSEITSsaSpectra::DCAcut(Double_t impactXY, Double_t impactZ zP[1]=59.8; zP[2]=1.2; } - } - else{ + }else{ if(fYear==2009){ xyP[0]=85.28;//DATA 900 GeV pass6 xyP[1]=25.78; @@ -223,24 +226,21 @@ Bool_t AliAnalysisTaskSEITSsaSpectra::DCAcut(Double_t impactXY, Double_t impactZ zP[2]=1.2; } } - Double_t xySigma = xyP[0] + xyP[1]/TMath::Power(TMath::Abs(pt),xyP[2]); - Double_t xyMax = fNSigmaDCAxy*xySigma; //in micron - if((TMath::Abs(impactXY)*10000)>xyMax) return kFALSE; - - Double_t zSigma = zP[0] + zP[1]/TMath::Power(TMath::Abs(pt),zP[2]); - Double_t zMax = fNSigmaDCAz*zSigma; //in micron - if((TMath::Abs(impactZ)*10000)>zMax) return kFALSE; - - return kTRUE; -} + fDCAxyCutFunc = new TF1("fDCAxyCutFunc","[3]*([0]+[1]/TMath::Power(TMath::Abs(x),[2]))",0.05,10.); + for(Int_t ipar=0; ipar<3; ipar++) fDCAxyCutFunc->SetParameter(ipar,xyP[ipar]); + fDCAxyCutFunc->SetParameter(3,fNSigmaDCAxy); -//________________________________________________________________________ -Double_t AliAnalysisTaskSEITSsaSpectra::Eta2y(Double_t pt, Double_t m, Double_t eta) const { - // convert eta to y - Double_t mt = TMath::Sqrt(m*m + pt*pt); - return TMath::ASinH(pt/mt*TMath::SinH(eta)); -} + + fDCAzCutFunc = new TF1("fDCAzCutFunc","[3]*([0]+[1]/TMath::Power(TMath::Abs(x),[2]))",0.05,10.); + for(Int_t ipar=0; ipar<3; ipar++) fDCAzCutFunc->SetParameter(ipar,zP[ipar]); + fDCAzCutFunc->SetParameter(3,fNSigmaDCAz); + + fListCuts->Add(fDCAxyCutFunc); + fListCuts->Add(fDCAzCutFunc); + + PostData(2,fListCuts); +} //________________________________________________________________________ void AliAnalysisTaskSEITSsaSpectra::UserCreateOutputObjects(){ @@ -251,14 +251,26 @@ void AliAnalysisTaskSEITSsaSpectra::UserCreateOutputObjects(){ fOutput->SetOwner(); fOutput->SetName("Spiderman"); - fHistNEvents = new TH1F("fHistNEvents", "Number of processed events",8,-1.5,6.5); + fHistNEvents = new TH1F("fHistNEvents", "Number of processed events",9,-1.5,7.5); fHistNEvents->Sumw2(); fHistNEvents->SetMinimum(0); + fHistNEvents->GetXaxis()->SetBinLabel(1,"Read from ESD"); + fHistNEvents->GetXaxis()->SetBinLabel(2,"Pass Phys. Sel. + Trig"); + fHistNEvents->GetXaxis()->SetBinLabel(3,"SDD read out"); + fHistNEvents->GetXaxis()->SetBinLabel(4,"In mult. range"); + fHistNEvents->GetXaxis()->SetBinLabel(5,"With SPD vertex"); + fHistNEvents->GetXaxis()->SetBinLabel(6,"Vertex contributors >0"); + fHistNEvents->GetXaxis()->SetBinLabel(7,"|zVertex|<10"); + fHistNEvents->GetXaxis()->SetBinLabel(8,"Error on zVertex<0.5"); + fHistNEvents->GetXaxis()->SetBinLabel(9,"Good Z vertex"); fOutput->Add(fHistNEvents); fHistMult = new TH1F("fHistMult", "Event Multiplicity",3000,-0.5,2999.5); fHistMult->Sumw2(); fHistMult->SetMinimum(0); + if(fMultEstimator==0) fHistMult->GetXaxis()->SetTitle("Multiplicity |#eta|<0.8"); + else if(fMultEstimator==1) fHistMult->GetXaxis()->SetTitle("Tracklets |#eta|<0.8"); + else if(fMultEstimator==2) fHistMult->GetXaxis()->SetTitle("Clusters on SPD1"); fOutput->Add(fHistMult); fHistCen = new TH1F("fHistCen", "Event Centrality",101,-0.5,100.5); @@ -542,6 +554,8 @@ void AliAnalysisTaskSEITSsaSpectra::UserCreateOutputObjects(){ fNtupleMC = new TNtuple("fNtupleMC","fNtupleMC","ptMC:pdgcode:signMC:etaMC:yMC:isph:evSel:run"); fOutput->Add(fNtupleMC); + + // Post output data. PostData(1,fOutput); @@ -577,20 +591,17 @@ void AliAnalysisTaskSEITSsaSpectra::UserExec(Option_t *){ } fHistNEvents->Fill(-1); - if(fLowEnergypp){ // remove events without SDD in pp 2.6 TeV - Bool_t hasSDD=kFALSE; - for (Int_t iTrack=0; iTrackGetNumberOfTracks(); iTrack++) { - track = (AliESDtrack*)fESD->GetTrack(iTrack); - if (!track) continue; - UInt_t clumap = track->GetITSClusterMap(); - if(clumap&4) hasSDD=kTRUE; - if(clumap&8) hasSDD=kTRUE; - if(hasSDD) break; - } - if(!hasSDD) return; - } + UInt_t maskPhysSel = ((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->IsEventSelected(); + TString firedTriggerClasses=fESD->GetFiredTriggerClasses(); + // if(!firedTriggerClasses.Contains("CINT1B")) return; + if((maskPhysSel & AliVEvent::kMB)==0) return; fHistNEvents->Fill(0); + if(fLowEnergypp && !fMC){ // remove events without SDD in pp 2.76 TeV + if(!firedTriggerClasses.Contains("ALL")) return; + } + fHistNEvents->Fill(1); + if(fMC){ AliMCEventHandler* eventHandler = dynamic_cast (AliAnalysisManager::GetAnalysisManager()->GetMCtruthEventHandler()); @@ -618,7 +629,7 @@ void AliAnalysisTaskSEITSsaSpectra::UserExec(Option_t *){ if(stack) nTrackMC = stack->GetNtrack(); const AliESDVertex *vtx = fESD->GetPrimaryVertexSPD(); -///////////selection of the centrality or multiplicity bin + ///////////selection of the centrality or multiplicity bin //selection on the event centrality if(fHImode){ @@ -634,74 +645,39 @@ void AliAnalysisTaskSEITSsaSpectra::UserExec(Option_t *){ } //selection on the event multiplicity based on global tracks - Int_t multiplicity = fesdTrackCutsMult->CountAcceptedTracks(fESD); - if(!fSPD){ - if(fLowMult>-1) - { - if(multiplicity-1) - { - if(multiplicity>fUpMult) - return; - } - - Printf("Multiplicity of the event (global tracks) : %i",multiplicity); - Printf("Multiplicity cut (global tracks) : %i to %i",fLowMult,fUpMult); - fHistMult->Fill(multiplicity); - } - - //multipicity selection based on SPD - if(fSPD){ - Float_t spdCorr=-1.0; + Int_t multiplicity = -1; + if(fMultEstimator==0){ + // tracks+tracklets + multiplicity=AliESDtrackCuts::GetReferenceMultiplicity(fESD, AliESDtrackCuts::kTrackletsITSTPC, 0.8); + }else if(fMultEstimator==1){ + // tracklets + multiplicity=AliESDtrackCuts::GetReferenceMultiplicity(fESD, AliESDtrackCuts::kTracklets, 0.8); + }else if(fMultEstimator==2){ + // clusters in SPD1 const AliMultiplicity *mult = fESD->GetMultiplicity(); - Float_t nClusters[6]={0.0,0.0,0.0,0.0,0.0,0.0}; - for(Int_t ilay=0; ilay<6; ilay++) - { - nClusters[ilay] = (Float_t)mult->GetNumberOfITSClusters(ilay); - } - spdCorr = AliESDUtils::GetCorrSPD2(nClusters[1],vtx->GetZ()); - { - if(fLowMult>-1) - { - if(((Int_t)spdCorr)-1) - { - if(((Int_t)spdCorr)>fUpMult) - { - return; - } - } - } - - Printf("Multiplicity of the event (SPD) : %i",(Int_t)spdCorr); - Printf("Multiplicity cut (SPD) : %i to %i",fLowMult,fUpMult); - fHistMult->Fill(spdCorr); + Float_t nClu1 = (Float_t)mult->GetNumberOfITSClusters(1); + multiplicity =(Int_t)(AliESDUtils::GetCorrSPD2(nClu1,vtx->GetZ())+0.5); } - - + if(fLowMult>-1 && multiplicity-1 && multiplicity>fUpMult) return; + fHistMult->Fill(multiplicity); + fHistNEvents->Fill(2); //event selection - fHistNEvents->Fill(1); if(!vtx)evSel=0; else{ - fHistNEvents->Fill(2); + fHistNEvents->Fill(3); if(vtx->GetNContributors()<0) evSel=0; else{ - fHistNEvents->Fill(3); + fHistNEvents->Fill(4); if(TMath::Abs(vtx->GetZv())>10) evSel=0; else{ - fHistNEvents->Fill(4); + fHistNEvents->Fill(5); if(vtx->GetZRes()>0.5) evSel=0; else{ - fHistNEvents->Fill(5); + fHistNEvents->Fill(6); if(vtx->IsFromVertexerZ() && vtx->GetDispersion()>0.03) evSel=0; - else fHistNEvents->Fill(6); + else fHistNEvents->Fill(7); } } } @@ -757,26 +733,26 @@ void AliAnalysisTaskSEITSsaSpectra::UserExec(Option_t *){ if(jpart>=0){ if(stack->IsPhysicalPrimary(imc)){ - if(signMC>0) fHistPrimMCposBefEvSel[jpart]->Fill(TMath::Abs(ptMC)); - else fHistPrimMCnegBefEvSel[jpart]->Fill(TMath::Abs(ptMC)); + if(signMC>0) fHistPrimMCposBefEvSel[jpart]->Fill(ptMC); + else fHistPrimMCnegBefEvSel[jpart]->Fill(ptMC); if(evSel==1){ - if(signMC>0) fHistPrimMCpos[jpart]->Fill(TMath::Abs(ptMC)); - else fHistPrimMCneg[jpart]->Fill(TMath::Abs(ptMC)); + if(signMC>0) fHistPrimMCpos[jpart]->Fill(ptMC); + else fHistPrimMCneg[jpart]->Fill(ptMC); } }else{ if(mfl==3 && uniqueID == kPDecay){ // If a particle is not a physical primary, check if it comes from weak decay - if(signMC>0) fHistSecStrMCposBefEvSel[jpart]->Fill(TMath::Abs(ptMC)); - else fHistSecStrMCnegBefEvSel[jpart]->Fill(TMath::Abs(ptMC)); + if(signMC>0) fHistSecStrMCposBefEvSel[jpart]->Fill(ptMC); + else fHistSecStrMCnegBefEvSel[jpart]->Fill(ptMC); if(evSel==1){ - if(signMC>0) fHistSecStrMCpos[jpart]->Fill(TMath::Abs(ptMC)); - else fHistSecStrMCneg[jpart]->Fill(TMath::Abs(ptMC)); + if(signMC>0) fHistSecStrMCpos[jpart]->Fill(ptMC); + else fHistSecStrMCneg[jpart]->Fill(ptMC); } }else{ - if(signMC>0) fHistSecMatMCposBefEvSel[jpart]->Fill(TMath::Abs(ptMC)); - else fHistSecMatMCnegBefEvSel[jpart]->Fill(TMath::Abs(ptMC)); + if(signMC>0) fHistSecMatMCposBefEvSel[jpart]->Fill(ptMC); + else fHistSecMatMCnegBefEvSel[jpart]->Fill(ptMC); if(evSel==1){ - if(signMC>0) fHistSecMatMCpos[jpart]->Fill(TMath::Abs(ptMC)); - else fHistSecMatMCneg[jpart]->Fill(TMath::Abs(ptMC)); + if(signMC>0) fHistSecMatMCpos[jpart]->Fill(ptMC); + else fHistSecMatMCneg[jpart]->Fill(ptMC); } } } @@ -1086,7 +1062,7 @@ void AliAnalysisTaskSEITSsaSpectra::UserExec(Option_t *){ } //DCA cut on xy and z - if(!DCAcut(impactXY,impactZ,pt,fMC)) continue; + if(!DCAcut(impactXY,impactZ,pt)) continue; label="DCA"; fHistNTracks->Fill(countBinTrk); @@ -1119,8 +1095,8 @@ void AliAnalysisTaskSEITSsaSpectra::UserExec(Option_t *){ else signMC=-1; ptMC=part->Pt(); if(stack->IsPhysicalPrimary(track->GetLabel())){ - if(signMC>0) fHistPrimMCposReco[jpart]->Fill(TMath::Abs(ptMC)); - else fHistPrimMCnegReco[jpart]->Fill(TMath::Abs(ptMC)); + if(signMC>0) fHistPrimMCposReco[jpart]->Fill(ptMC); + else fHistPrimMCnegReco[jpart]->Fill(ptMC); }else{ Int_t indexMoth=part->GetFirstMother(); if(indexMoth>=0){ @@ -1130,11 +1106,11 @@ void AliAnalysisTaskSEITSsaSpectra::UserExec(Option_t *){ } uniqueID = part->GetUniqueID(); if(mfl==3 && uniqueID == kPDecay){ // strangeness - if(signMC>0) fHistSecStrMCposReco[jpart]->Fill(TMath::Abs(ptMC)); - else fHistSecStrMCnegReco[jpart]->Fill(TMath::Abs(ptMC)); + if(signMC>0) fHistSecStrMCposReco[jpart]->Fill(ptMC); + else fHistSecStrMCnegReco[jpart]->Fill(ptMC); }else{ - if(signMC>0) fHistSecMatMCposReco[jpart]->Fill(TMath::Abs(ptMC)); - else fHistSecMatMCnegReco[jpart]->Fill(TMath::Abs(ptMC)); + if(signMC>0) fHistSecMatMCposReco[jpart]->Fill(ptMC); + else fHistSecMatMCnegReco[jpart]->Fill(ptMC); } } } @@ -1265,7 +1241,8 @@ void AliAnalysisTaskSEITSsaSpectra::UserExec(Option_t *){ // Post output data. PostData(1,fOutput); - Printf("............. end of Exec"); + PostData(2,fListCuts); + // Printf("............. end of Exec"); } //________________________________________________________________________ @@ -1279,147 +1256,7 @@ void AliAnalysisTaskSEITSsaSpectra::Terminate(Option_t *) { return; } fHistNEvents = dynamic_cast(fOutput->FindObject("fHistNEvents")); - fHistMult = dynamic_cast(fOutput->FindObject("fHistMult")); - fHistCen = dynamic_cast(fOutput->FindObject("fHistCen")); - fHistNTracks = dynamic_cast(fOutput->FindObject("fHistNTracks")); - fHistNTracksPos = dynamic_cast(fOutput->FindObject("fHistNTracksPos")); - fHistNTracksNeg = dynamic_cast(fOutput->FindObject("fHistNTracksNeg")); - fHistDEDX = dynamic_cast(fOutput->FindObject("fHistDEDX")); - fHistDEDXdouble = dynamic_cast(fOutput->FindObject("fHistDEDXdouble")); - - - fHistBeforeEvSel = dynamic_cast(fOutput->FindObject("fHistBeforeEvSel")); - fHistAfterEvSel = dynamic_cast(fOutput->FindObject("fHistAfterEvSel")); - - - for(Int_t j=0;j<3;j++){ - fHistPosNSigmaSep[j] = dynamic_cast(fOutput->FindObject(Form("fHistPosNSigmaSep%d",j))); - fHistNegNSigmaSep[j] = dynamic_cast(fOutput->FindObject(Form("fHistNegNSigmaSep%d",j))); - if(fMC){ - fHistPrimMCpos[j] = dynamic_cast(fOutput->FindObject(Form("fHistPrimMCpos%d",j))); - fHistPrimMCneg[j] = dynamic_cast(fOutput->FindObject(Form("fHistPrimMCneg%d",j))); - fHistSecStrMCpos[j] = dynamic_cast(fOutput->FindObject(Form("fHistSecStrMCpos%d",j))); - fHistSecStrMCneg[j] = dynamic_cast(fOutput->FindObject(Form("fHistSecStrMCneg%d",j))); - fHistSecMatMCpos[j] = dynamic_cast(fOutput->FindObject(Form("fHistSecMatMCpos%d",j))); - fHistSecMatMCneg[j] = dynamic_cast(fOutput->FindObject(Form("fHistSecMatMCneg%d",j))); - // - fHistPrimMCposBefEvSel[j] = dynamic_cast(fOutput->FindObject(Form("fHistPrimMCposBefEvSel%d",j))); - fHistPrimMCnegBefEvSel[j] = dynamic_cast(fOutput->FindObject(Form("fHistPrimMCnegBefEvSel%d",j))); - fHistSecStrMCposBefEvSel[j] = dynamic_cast(fOutput->FindObject(Form("fHistSecStrMCposBefEvSel%d",j))); - fHistSecStrMCnegBefEvSel[j] = dynamic_cast(fOutput->FindObject(Form("fHistSecStrMCnegBefEvSel%d",j))); - fHistSecMatMCposBefEvSel[j] = dynamic_cast(fOutput->FindObject(Form("fHistSecMatMCposBefEvSel%d",j))); - fHistSecMatMCnegBefEvSel[j] = dynamic_cast(fOutput->FindObject(Form("fHistSecMatMCnegBefEvSel%d",j))); - // - fHistPrimMCposReco[j] = dynamic_cast(fOutput->FindObject(Form("fHistPrimMCposReco%d",j))); - fHistPrimMCnegReco[j] = dynamic_cast(fOutput->FindObject(Form("fHistPrimMCnegReco%d",j))); - fHistSecStrMCposReco[j] = dynamic_cast(fOutput->FindObject(Form("fHistSecStrMCposReco%d",j))); - fHistSecStrMCnegReco[j] = dynamic_cast(fOutput->FindObject(Form("fHistSecStrMCnegReco%d",j))); - fHistSecMatMCposReco[j] = dynamic_cast(fOutput->FindObject(Form("fHistSecMatMCposReco%d",j))); - fHistSecMatMCnegReco[j] = dynamic_cast(fOutput->FindObject(Form("fHistSecMatMCnegReco%d",j))); - } - } - - for(Int_t i=0; i<4; i++){ - fHistCharge[i] = dynamic_cast(fOutput->FindObject(Form("fHistChargeLay%d",i))); - } - - for(Int_t i=0; i(fOutput->FindObject(Form("fHistPosPi%d",i))); - fHistPosK[i] = dynamic_cast(fOutput->FindObject(Form("fHistPosK%d",i))); - fHistPosP[i] = dynamic_cast(fOutput->FindObject(Form("fHistPosP%d",i))); - fHistNegPi[i] = dynamic_cast(fOutput->FindObject(Form("fHistNegPi%d",i))); - fHistNegK[i] = dynamic_cast(fOutput->FindObject(Form("fHistNegK%d",i))); - fHistNegP[i] = dynamic_cast(fOutput->FindObject(Form("fHistNegP%d",i))); - - fHistDCAPosPi[i] = dynamic_cast(fOutput->FindObject(Form("fHistDCAPosPi%d",i))); - fHistDCAPosK[i] = dynamic_cast(fOutput->FindObject(Form("fHistDCAPosK%d",i))); - fHistDCAPosP[i] = dynamic_cast(fOutput->FindObject(Form("fHistDCAPosP%d",i))); - fHistDCANegPi[i] = dynamic_cast(fOutput->FindObject(Form("fHistDCANegPi%d",i))); - fHistDCANegK[i] = dynamic_cast(fOutput->FindObject(Form("fHistDCANegK%d",i))); - fHistDCANegP[i] = dynamic_cast(fOutput->FindObject(Form("fHistDCANegP%d",i))); - - - if(fMC){ - - fHistMCPrimDCAPosPi[i] = dynamic_cast(fOutput->FindObject(Form("fHistMCPrimDCAPosPi%d",i))); - fHistMCPrimDCAPosK[i] = dynamic_cast(fOutput->FindObject(Form("fHistMCPrimDCAPosK%d",i))); - fHistMCPrimDCAPosP[i] = dynamic_cast(fOutput->FindObject(Form("fHistMCPrimDCAPosP%d",i))); - fHistMCPrimDCANegPi[i] = dynamic_cast(fOutput->FindObject(Form("fHistMCPrimDCANegPi%d",i))); - fHistMCPrimDCANegK[i] = dynamic_cast(fOutput->FindObject(Form("fHistMCPrimDCANegK%d",i))); - fHistMCPrimDCANegP[i] = dynamic_cast(fOutput->FindObject(Form("fHistMCPrimDCANegP%d",i))); - - fHistMCSecStDCAPosPi[i] = dynamic_cast(fOutput->FindObject(Form("fHistMCSecStDCAPosPi%d",i))); - fHistMCSecStDCAPosK[i] = dynamic_cast(fOutput->FindObject(Form("fHistMCSecStDCAPosK%d",i))); - fHistMCSecStDCAPosP[i] = dynamic_cast(fOutput->FindObject(Form("fHistMCSecStDCAPosP%d",i))); - fHistMCSecStDCANegPi[i] = dynamic_cast(fOutput->FindObject(Form("fHistMCSecStDCANegPi%d",i))); - fHistMCSecStDCANegK[i] = dynamic_cast(fOutput->FindObject(Form("fHistMCSecStDCANegK%d",i))); - fHistMCSecStDCANegP[i] = dynamic_cast(fOutput->FindObject(Form("fHistMCSecStDCANegP%d",i))); - - fHistMCSecMatDCAPosPi[i] = dynamic_cast(fOutput->FindObject(Form("fHistMCSecMatDCAPosPi%d",i))); - fHistMCSecMatDCAPosK[i] = dynamic_cast(fOutput->FindObject(Form("fHistMCSecMatDCAPosK%d",i))); - fHistMCSecMatDCAPosP[i] = dynamic_cast(fOutput->FindObject(Form("fHistMCSecMatDCAPosP%d",i))); - fHistMCSecMatDCANegPi[i] = dynamic_cast(fOutput->FindObject(Form("fHistMCSecMatDCANegPi%d",i))); - fHistMCSecMatDCANegK[i] = dynamic_cast(fOutput->FindObject(Form("fHistMCSecMatDCANegK%d",i))); - fHistMCSecMatDCANegP[i] = dynamic_cast(fOutput->FindObject(Form("fHistMCSecMatDCANegP%d",i))); - - fHistMCPosOtherHypPion[i] = dynamic_cast(fOutput->FindObject(Form("fHistMCPosOtherHypPion%d",i))); - fHistMCPosOtherHypKaon[i] = dynamic_cast(fOutput->FindObject(Form("fHistMCPosOtherHypKaon%d",i))); - fHistMCPosOtherHypProton[i] = dynamic_cast(fOutput->FindObject(Form("fHistMCPosOtherHypProton%d",i))); - fHistMCPosElHypPion[i] = dynamic_cast(fOutput->FindObject(Form("fHistMCPosElHypPion%d",i))); - fHistMCPosElHypKaon[i] = dynamic_cast(fOutput->FindObject(Form("fHistMCPosElHypKaon%d",i))); - fHistMCPosElHypProton[i] = dynamic_cast(fOutput->FindObject(Form("fHistMCPosElHypProton%d",i))); - fHistMCPosPiHypPion[i] = dynamic_cast(fOutput->FindObject(Form("fHistMCPosPiHypPion%d",i))); - fHistMCPosPiHypKaon[i] = dynamic_cast(fOutput->FindObject(Form("fHistMCPosPiHypKaon%d",i))); - fHistMCPosPiHypProton[i] = dynamic_cast(fOutput->FindObject(Form("fHistMCPosPiHypProton%d",i))); - fHistMCPosKHypPion[i] = dynamic_cast(fOutput->FindObject(Form("fHistMCPosKHypPion%d",i))); - fHistMCPosKHypKaon[i] = dynamic_cast(fOutput->FindObject(Form("fHistMCPosKHypKaon%d",i))); - fHistMCPosKHypProton[i] = dynamic_cast(fOutput->FindObject(Form("fHistMCPosKHypProton%d",i))); - fHistMCPosPHypPion[i] = dynamic_cast(fOutput->FindObject(Form("fHistMCPosPHypPion%d",i))); - fHistMCPosPHypKaon[i] = dynamic_cast(fOutput->FindObject(Form("fHistMCPosPHypKaon%d",i))); - fHistMCPosPHypProton[i] = dynamic_cast(fOutput->FindObject(Form("fHistMCPosPHypProton%d",i))); - - fHistMCNegOtherHypPion[i] = dynamic_cast(fOutput->FindObject(Form("fHistMCNegOtherHypPion%d",i))); - fHistMCNegOtherHypKaon[i] = dynamic_cast(fOutput->FindObject(Form("fHistMCNegOtherHypKaon%d",i))); - fHistMCNegOtherHypProton[i] = dynamic_cast(fOutput->FindObject(Form("fHistMCNegOtherHypProton%d",i))); - fHistMCNegElHypPion[i] = dynamic_cast(fOutput->FindObject(Form("fHistMCNegElHypPion%d",i))); - fHistMCNegElHypKaon[i] = dynamic_cast(fOutput->FindObject(Form("fHistMCNegElHypKaon%d",i))); - fHistMCNegElHypProton[i] = dynamic_cast(fOutput->FindObject(Form("fHistMCNegElHypProton%d",i))); - fHistMCNegPiHypPion[i] = dynamic_cast(fOutput->FindObject(Form("fHistMCNegPiHypPion%d",i))); - fHistMCNegPiHypKaon[i] = dynamic_cast(fOutput->FindObject(Form("fHistMCNegPiHypKaon%d",i))); - fHistMCNegPiHypProton[i] = dynamic_cast(fOutput->FindObject(Form("fHistMCNegPiHypProton%d",i))); - fHistMCNegKHypPion[i] = dynamic_cast(fOutput->FindObject(Form("fHistMCNegKHypPion%d",i))); - fHistMCNegKHypKaon[i] = dynamic_cast(fOutput->FindObject(Form("fHistMCNegKHypKaon%d",i))); - fHistMCNegKHypProton[i] = dynamic_cast(fOutput->FindObject(Form("fHistMCNegKHypProton%d",i))); - fHistMCNegPHypPion[i] = dynamic_cast(fOutput->FindObject(Form("fHistMCNegPHypPion%d",i))); - fHistMCNegPHypKaon[i] = dynamic_cast(fOutput->FindObject(Form("fHistMCNegPHypKaon%d",i))); - fHistMCNegPHypProton[i] = dynamic_cast(fOutput->FindObject(Form("fHistMCNegPHypProton%d",i))); - - } - } - - for(Int_t j=0;j<3;j++){ - fHistPosNSigmaMean[j] = dynamic_cast(fOutput->FindObject(Form("hHistPosNSigmaMean%d",j))); - fHistNegNSigmaMean[j] = dynamic_cast(fOutput->FindObject(Form("hHistNegNSigmaMean%d",j))); - fHistPosNSigmaMCMean[j] = dynamic_cast(fOutput->FindObject(Form("hHistPosNSigmaMCMean%d",j))); - fHistNegNSigmaMCMean[j] = dynamic_cast(fOutput->FindObject(Form("hHistNegNSigmaMCMean%d",j))); - fHistPosNSigmaPrimMean[j] = dynamic_cast(fOutput->FindObject(Form("hHistPosNSigmaPrimMean%d",j))); - fHistNegNSigmaPrimMean[j] = dynamic_cast(fOutput->FindObject(Form("hHistNegNSigmaPrimMean%d",j))); - fHistPosNSigmaPrimMCMean[j] = dynamic_cast(fOutput->FindObject(Form("hHistPosNSigmaPrimMCMean%d",j))); - fHistNegNSigmaPrimMCMean[j] = dynamic_cast(fOutput->FindObject(Form("hHistNegNSigmaPrimMCMean%d",j))); - - fHistPosNSigma[j] = dynamic_cast(fOutput->FindObject(Form("hHistPosNSigma%d",j))); - fHistNegNSigma[j] = dynamic_cast(fOutput->FindObject(Form("hHistNegNSigma%d",j))); - fHistPosNSigmaMC[j] = dynamic_cast(fOutput->FindObject(Form("hHistPosNSigmaMC%d",j))); - fHistNegNSigmaMC[j] = dynamic_cast(fOutput->FindObject(Form("hHistNegNSigmaMC%d",j))); - fHistPosNSigmaPrim[j] = dynamic_cast(fOutput->FindObject(Form("hHistPosNSigmaPrim%d",j))); - fHistNegNSigmaPrim[j] = dynamic_cast(fOutput->FindObject(Form("hHistNegNSigmaPrim%d",j))); - fHistPosNSigmaPrimMC[j] = dynamic_cast(fOutput->FindObject(Form("hHistPosNSigmaPrimMC%d",j))); - fHistNegNSigmaPrimMC[j] = dynamic_cast(fOutput->FindObject(Form("hHistNegNSigmaPrimMC%d",j))); - - } - - fNtupleNSigma = dynamic_cast(fOutput->FindObject("fNtupleNSigma")); - fNtupleMC = dynamic_cast(fOutput->FindObject("fNtupleMC")); + printf("Number of Analyzed Events = %f\n",fHistNEvents->GetBinContent(1)); Printf("end of Terminate"); return; diff --git a/PWGLF/SPECTRA/PiKaPr/ITSsa/AliAnalysisTaskSEITSsaSpectra.h b/PWGLF/SPECTRA/PiKaPr/ITSsa/AliAnalysisTaskSEITSsaSpectra.h index fc31f0f975d..b535ff88af4 100644 --- a/PWGLF/SPECTRA/PiKaPr/ITSsa/AliAnalysisTaskSEITSsaSpectra.h +++ b/PWGLF/SPECTRA/PiKaPr/ITSsa/AliAnalysisTaskSEITSsaSpectra.h @@ -17,11 +17,11 @@ class TString; class TTree; class TH1F; +class TF1; class TH2F; class TRandom3; class AliESDEvent; class TNtuple; -class AliESDtrackCuts; class AliITSPIDResponse; #include "AliAnalysisTaskSE.h" @@ -32,6 +32,8 @@ class AliAnalysisTaskSEITSsaSpectra : public AliAnalysisTaskSE { virtual ~AliAnalysisTaskSEITSsaSpectra(); virtual void UserCreateOutputObjects(); + virtual void Init(); + virtual void LocalInit() {Init();} virtual void UserExec(Option_t *); virtual void Terminate(Option_t *); @@ -66,14 +68,16 @@ class AliAnalysisTaskSEITSsaSpectra : public AliAnalysisTaskSE { fLowCentrality=low; fUpCentrality=up; } } - void SetSPDMethodCut(){fSPD=kTRUE;} void SetHImode(){fHImode=kTRUE;} + void SetUseTrackMultiplicityEstimator(){fMultEstimator=0;} + void SetUseTrackletsMultiplicityEstimator(){fMultEstimator=1;} + void SetUseClustersSPD1MultiplicityEstimator(){fMultEstimator=2;} + void SetEtaMax(Double_t maxeta){ fEtaRange=maxeta; } - void SetYear(Int_t year); void SetReadMC(Bool_t flag = kTRUE) {fMC = flag;} void SetFillNtuple(Bool_t fill=kTRUE) {fFillNtuple=fill;} void SetLowEnergypp(Bool_t opt=kTRUE) {fLowEnergypp=opt;} @@ -85,7 +89,7 @@ class AliAnalysisTaskSEITSsaSpectra : public AliAnalysisTaskSE { } Double_t CookdEdx(Double_t *s) const; Double_t Eta2y(Double_t pt, Double_t m, Double_t eta) const; - Bool_t DCAcut(Double_t impactXY, Double_t impactZ, Double_t pt, Bool_t optMC) const; + Bool_t DCAcut(Double_t impactXY, Double_t impactZ, Double_t pt) const; private: AliAnalysisTaskSEITSsaSpectra(const AliAnalysisTaskSEITSsaSpectra &source); @@ -94,9 +98,9 @@ class AliAnalysisTaskSEITSsaSpectra : public AliAnalysisTaskSE { enum {kNbins=22}; AliESDEvent *fESD; //ESD object - AliESDtrackCuts *fesdTrackCutsMult;//cuts for multiplicity TList *fOutput; //! tlist with output + TList *fListCuts; // list of functions storing DCA cut TH1F *fHistNEvents; //! histo with number of events TH1F *fHistMult; //! histo with multiplicity of the events TH1F *fHistCen; //! histo with multiplicity of the events @@ -218,6 +222,9 @@ class AliAnalysisTaskSEITSsaSpectra : public AliAnalysisTaskSE { TH1F *fHistNegNSigmaPrim[3]; //! NSigma histos for 6 species TH1F *fHistNegNSigmaPrimMC[3]; //! NSigma histos for 6 species + TF1* fDCAxyCutFunc; // function with DCAz cut vs. pt + TF1* fDCAzCutFunc; // function with DCAxy cut vs. pt + Double_t fPtBinLimits[kNbins+1]; // limits of Pt Bins AliITSPIDResponse* fITSPIDResponse; //! class with BB parameterizations Int_t fMinSPDPts; // minimum number of SPD Points @@ -233,7 +240,7 @@ class AliAnalysisTaskSEITSsaSpectra : public AliAnalysisTaskSE { Int_t fUpMult; // Multiplicity bin Float_t fLowCentrality;//low Centrality cut Float_t fUpCentrality;//up Centrality cut - Bool_t fSPD;//use spd2 as mulestimator + Int_t fMultEstimator; // multiplicty estimator Bool_t fHImode;//use spd2 as mulestimator Int_t fYear; // Year (2009, 2010) Bool_t fMC; //flag to switch on the MC analysis for the efficiency estimation diff --git a/PWGLF/SPECTRA/PiKaPr/ITSsa/FitDCADistributions.C b/PWGLF/SPECTRA/PiKaPr/ITSsa/FitDCADistributions.C new file mode 100644 index 00000000000..3a4e1285f83 --- /dev/null +++ b/PWGLF/SPECTRA/PiKaPr/ITSsa/FitDCADistributions.C @@ -0,0 +1,303 @@ +#if !defined(__CINT__) || defined(__MAKECINT__) +#include +#include +#include +#include +#include +#include +#include +#include +#endif + +enum species {kPion, kKaon, kProton}; +enum chargesign {kPositive,kNegative}; + +Float_t PerformFit(TH1F* fHistDCA, TObjArray *mc,TString partname,Float_t xyMax); + +void FitDCADistributions(TString period="LHC10d", + TString MCname="LHC10f6aold", + Int_t iSpecies=2, + Int_t cSign=1 + ) +{ + gROOT->SetStyle("Plain"); + gStyle->SetFillColor(0); + gStyle->SetOptStat(0000); + //binning + const Int_t nptbins = 22; + Double_t ptbins[nptbins+1]={0.08,0.10,0.12,0.14,0.16,0.18,0.20,0.25,0.30,0.35,0.40,0.45,0.50,0.55,0.60,0.65,0.70,0.75,0.80,0.85,0.90,0.95,1.0}; + //Final correction + TH1F *hPrimAllDATA=new TH1F("hPrimAllDATA","hPrimAllDATA",nptbins,ptbins); + TH1F *hPrimAllMC=new TH1F("hPrimAllMC","hPrimAllMC",nptbins,ptbins); + TH1F *hPrimAllDATAMC=new TH1F("hPrimAllDATAMC","hPrimAllDATAMC",nptbins,ptbins); + + + TString fnameMC=Form("/home/prino/alice/PID/pp7/%s/AnalysisResults.root",MCname.Data()); + TString lnameMC="clistITSsaMult-1to-1"; + TString lnameCutMC="DCACutMult-1to-1"; + TString fnameDATA=Form("/home/prino/alice/PID/pp7/%s/AnalysisResults.root",period.Data()); + TString lnameDATA="clistITSsaMult-1to-1"; + TString lnameCutDATA="DCACutMult-1to-1"; + + const Int_t firstbin=1; + Int_t rebin=10; + Bool_t optSmooth=kTRUE; + + + TString pCharge="Pos"; + TString pcharge="pos"; + if(cSign==kNegative){ + pCharge="Neg"; + pcharge="neg"; + } + TString species[3]={"Pi","K","P"}; + TString speciesRoot[6]={"#pi^{+}","#pi^{-}","K^{+}","K^{-}","p","#bar{p}"}; + Int_t idPart=iSpecies*2+cSign; + TString partname=Form("%s%s",pCharge.Data(),species[iSpecies].Data()); + printf("%s\n",partname.Data()); + + // MC histograms + TFile *fMC=new TFile(fnameMC,"READ"); + TH1F *fHistDCAMC[nptbins]; + TH1F *fHistDCAPrim[nptbins]; + TH1F *fHistDCAsec[nptbins]; + TH1F *fHistDCAsecSt[nptbins]; + TDirectory *dirFileMC=(TDirectory*)fMC->Get("PWG2SpectraITSsa"); + TList *liMC = (TList*)dirFileMC->Get(lnameMC.Data()); + for(Int_t m=0;mFindObject(Form("fHistDCA%s%i",partname.Data(),m)); + fHistDCAPrim[m] = (TH1F*)liMC->FindObject(Form("fHistMCPrimDCA%s%i",partname.Data(),m)); + fHistDCAsec[m] = (TH1F*)liMC->FindObject(Form("fHistMCSecMatDCA%s%i",partname.Data(),m)); + fHistDCAsecSt[m] = (TH1F*)liMC->FindObject(Form("fHistMCSecStDCA%s%i",partname.Data(),m)); + fHistDCAMC[m]->Rebin(rebin); + fHistDCAPrim[m]->Rebin(rebin); + fHistDCAsec[m]->Rebin(rebin); + fHistDCAsecSt[m]->Rebin(rebin); + } + TList *liCutMC = (TList*)dirFileMC->Get(lnameCutMC.Data()); + Bool_t setDCA=kFALSE; + + Float_t DCAcut=7.; + Double_t xyPMC[3]; //parameters of DCA cut + xyPMC[0]=36.; //MC LHC10d1 + xyPMC[1]=43.9; + xyPMC[2]=1.3; + TF1* fDCAxyCutMC=0x0; + TF1* fDCAzCutMC=0x0; + Float_t nSigmaDCAxyMC=0.; + Float_t nSigmaDCAzMC=0.; + if(liCutMC){ + fDCAxyCutMC=(TF1*)liCutMC->FindObject("fDCAxyCutFunc"); + fDCAzCutMC=(TF1*)liCutMC->FindObject("fDCAzCutFunc"); + nSigmaDCAxyMC=fDCAxyCutMC->GetParameter(3); + nSigmaDCAzMC=fDCAzCutMC->GetParameter(3); + printf("DCA cut parameters from TF1: %f %f %f --- DCA cut = %f\n",fDCAxyCutMC->GetParameter(0),fDCAxyCutMC->GetParameter(1),fDCAxyCutMC->GetParameter(2),nSigmaDCAxyMC); + setDCA=kTRUE; + } + if(!setDCA){ + printf("DCA cut values for MC not found in the output file, use dafault values. DCAcut=%f with default resolution\n",DCAcut); + fDCAxyCutMC = new TF1("fDCAxyCutFunc","[3]*([0]+[1]/TMath::Power(TMath::Abs(x),[2]))",0.05,10.); + for(Int_t ipar=0; ipar<3; ipar++) fDCAxyCutMC->SetParameter(ipar,xyPMC[ipar]); + fDCAxyCutMC->SetParameter(3,DCAcut); + nSigmaDCAxyMC=DCAcut; + } + + //DATA histograms + TFile *fDATA=new TFile(fnameDATA,"READ"); + TDirectory *dirFileDATA=(TDirectory*)fDATA->Get("PWG2SpectraITSsa"); + TList *liDATA = (TList*)dirFileDATA->Get(lnameDATA.Data()); + TH1F *fHistDCADATA[nptbins]; + for(Int_t m=0;mFindObject(Form("fHistDCA%s%i",partname.Data(),m)); + fHistDCADATA[m]->Rebin(rebin); + } + TList *liCutDATA = (TList*)dirFileDATA->Get(lnameCutDATA.Data()); + setDCA=kFALSE; + Float_t nSigmaDCAxyDATA=0.; + Float_t nSigmaDCAzDATA=0.; + Double_t xyP[3]; // parameters of DCA cut + xyP[0]=32.7; //DATA 7 TeV pass2 + xyP[1]=44.8; + xyP[2]=1.3; + TF1* fDCAxyCutDATA=0x0; + TF1* fDCAzCutDATA=0x0; + if(liCutDATA){ + fDCAxyCutDATA=(TF1*)liCutDATA->FindObject("fDCAxyCutFunc"); + fDCAzCutDATA=(TF1*)liCutDATA->FindObject("fDCAzCutFunc"); + nSigmaDCAxyDATA=fDCAxyCutDATA->GetParameter(3); + nSigmaDCAzDATA=fDCAzCutDATA->GetParameter(3); + printf("DCA cut parameters from TF1: %f %f %f --- DCA cut = %f\n",fDCAxyCutDATA->GetParameter(0),fDCAxyCutDATA->GetParameter(1),fDCAxyCutDATA->GetParameter(2),nSigmaDCAxyDATA); + setDCA=kTRUE; + } + if(!setDCA){ + printf("DCA cut values for DATA not found in the output file, use dafault values. DCAcut=%f with default resolution\n",DCAcut); + fDCAxyCutDATA = new TF1("fDCAxyCutFunc","[3]*([0]+[1]/TMath::Power(TMath::Abs(x),[2]))",0.05,10.); + for(Int_t ipar=0; ipar<3; ipar++) fDCAxyCutDATA->SetParameter(ipar,xyP[ipar]); + fDCAxyCutDATA->SetParameter(3,DCAcut); + nSigmaDCAxyDATA=DCAcut; + } + + if(TMath::Abs(nSigmaDCAxyDATA-nSigmaDCAxyMC)<0.01) DCAcut=nSigmaDCAxyDATA; + else{ + printf("ERROR: DCAxy cuts do not match between data and MC\n"); + return; + } + if(TMath::Abs(nSigmaDCAzDATA-nSigmaDCAzMC)<0.01){ + printf("ERROR: DCAz cuts do not match between data and MC\n"); + return; + } + TCanvas *cfitDATA=new TCanvas("cfitDATA","cfitDATA"); + cfitDATA->Divide(6,4); + TCanvas *cfitMC=new TCanvas("cfitMC","cfitMC"); + cfitMC->Divide(6,4); + + for(Int_t ibin=firstbin;ibinEval(TMath::Abs((ptbins[ibin+1]+ptbins[ibin])/2))/10000.; + Double_t xyMaxMC =fDCAxyCutMC->Eval(TMath::Abs((ptbins[ibin+1]+ptbins[ibin])/2))/10000.; + + TObjArray *mcTemplates = 0x0; + if(partname=="PosP"){ + mcTemplates = new TObjArray(3); // MC histograms are put in this array + mcTemplates->Add(fHistDCAPrim[ibin]); + mcTemplates->Add(fHistDCAsecSt[ibin]); + mcTemplates->Add(fHistDCAsec[ibin]); + }else{ + mcTemplates = new TObjArray(2); // MC histograms are put in this array + mcTemplates->Add(fHistDCAPrim[ibin]); + mcTemplates->Add(fHistDCAsecSt[ibin]); + } + if(fHistDCADATA[ibin]->Integral()==0) continue; + + ////////////////// Fit on DATA + cfitDATA->cd(ibin); + gPad->SetLogy(); + Double_t primAllDATA=PerformFit(fHistDCADATA[ibin],mcTemplates,partname,xyMax); + + ////////////////// Fit on MC + cfitMC->cd(ibin); + gPad->SetLogy(); + Double_t primAllMC=PerformFit(fHistDCAMC[ibin],mcTemplates,partname,xyMaxMC); + + if(hPrimAllMC!=0 && hPrimAllDATA!=0){ + hPrimAllDATA->Fill((ptbins[ibin]+ptbins[ibin+1])/2,primAllDATA); + hPrimAllMC->Fill((ptbins[ibin]+ptbins[ibin+1])/2,primAllMC); + hPrimAllDATAMC->Fill((ptbins[ibin]+ptbins[ibin+1])/2,primAllDATA/primAllMC); + } + } + + //Adding MC truth + TH1F* hAll=(TH1F*)liMC->FindObject(Form("hHist%sNSigmaMean%d",pCharge.Data(),iSpecies)); + TH1F* hPrim=(TH1F*)liMC->FindObject(Form("hHist%sNSigmaPrimMean%d",pCharge.Data(),iSpecies)); + TH1F* hPrimRecMC=(TH1F*)liMC->FindObject(Form("fHistPrimMC%sReco%d",pcharge.Data(),iSpecies)); + TH1F* hSecMatRecMC=(TH1F*)liMC->FindObject(Form("fHistSecMatMC%sReco%d",pcharge.Data(),iSpecies)); + TH1F* hSecStrRecMC=(TH1F*)liMC->FindObject(Form("fHistSecStrMC%sReco%d",pcharge.Data(),iSpecies)); + hSecStrRecMC->Add(hSecMatRecMC); + hSecStrRecMC->Add(hPrimRecMC); + hPrim->SetLineColor(8); + hPrim->SetLineWidth(2); + hPrim->SetTitle("Prim/All from MC Truth"); + hPrim->Divide(hAll); + hPrimRecMC->SetLineColor(11); + hPrimRecMC->SetLineWidth(2); + hPrimRecMC->SetTitle("Prim/All from MC Truth Reco"); + hPrimRecMC->Divide(hSecStrRecMC); + + hPrimAllDATAMC->GetXaxis()->SetRangeUser(0.09,0.79); + hPrimAllMC->GetXaxis()->SetRangeUser(0.09,0.79); + hPrimAllDATA->GetXaxis()->SetRangeUser(0.09,0.79); + if(optSmooth) hPrimAllDATAMC->Smooth(1,"R"); + + hPrimAllDATA->SetTitle("Prim/all data"); + hPrimAllDATA->SetLineColor(2); + hPrimAllDATA->SetLineWidth(2); + hPrimAllMC->SetTitle("Prim/all mc"); + hPrimAllMC->SetLineColor(4); + hPrimAllMC->SetLineWidth(2); + hPrimAllDATAMC->SetTitle("Prim/all data/mc"); + hPrimAllDATAMC->SetLineColor(1); + hPrimAllDATAMC->SetLineWidth(2); + + TCanvas *cPrimAll=new TCanvas("cPrimAll","cPrimAll"); + hPrimAllDATA->SetMinimum(0.7); + hPrimAllDATA->SetMaximum(1.1); + hPrimAllDATA->Draw("l"); + hPrimAllMC->Draw("lsame"); + hPrimAllDATAMC->Draw("lsame"); + hPrim->Draw("lsame"); + hPrimRecMC->Draw("lsame"); + gPad->BuildLegend(); + TLatex* tsp=new TLatex(0.4,0.75,speciesRoot[idPart].Data()); + tsp->SetTextFont(43); + tsp->SetTextSize(28); + tsp->SetNDC(); + tsp->Draw(); + cPrimAll->Update(); + + TString fout; + fout=Form("DCACorr%s_%s_%.0fDCA_%s_TFraction.root",period.Data(),MCname.Data(),DCAcut,partname.Data()); + TFile *out=new TFile(fout.Data(),"recreate"); + hPrimAllDATA->Write(); + hPrimAllMC->Write(); + hPrimAllDATAMC->Write(); + out->Close(); + delete out; + +} + +Float_t PerformFit(TH1F* fHistDCA, TObjArray *mc,TString partname,Float_t xyMax){ + + Double_t prim=0,secSt=0,sec=0; + TFractionFitter* fit = new TFractionFitter(fHistDCA, mc); // initialise + fit->Constrain(0,0.0,1.0); // constrain fraction 1 to be between 0 and 1 + fit->Constrain(1,0.0,1.0); // constrain fraction 1 to be between 0 and 1 + if(partname=="PosP")fit->Constrain(2,0.0,1.0); // constrain fraction 1 to be between 0 and 1 + + //fit->SetRangeX(200,1000); // use only the first 15 bins in the fit + Int_t status = fit->Fit(); // perform the fit + cout << "fit status: " << status << endl; + + if (status == 0) { // check on fit status + TH1F* result = (TH1F*) fit->GetPlot(); + TH1F* PrimMCPred=(TH1F*)fit->GetMCPrediction(0); + TH1F* secStMCPred=(TH1F*)fit->GetMCPrediction(1); + TH1F* secMCPred=0x0; + if(partname=="PosP"){ + secMCPred=(TH1F*)fit->GetMCPrediction(2); + secMCPred->SetLineColor(4); + } + //Drawing section + PrimMCPred->SetLineColor(2); + secStMCPred->SetLineColor(6); + fHistDCA->SetTitle("DCA distr - data"); + fHistDCA->DrawCopy("Ep"); + result->SetTitle("Fit result"); + result->DrawCopy("same"); + Double_t value,error; + + fit->GetResult(0,value,error); + PrimMCPred->Scale(fHistDCA->GetSumOfWeights()*value/PrimMCPred->GetSumOfWeights()); + PrimMCPred->SetTitle("Primaries"); + PrimMCPred->DrawCopy("same"); + fit->GetResult(1,value,error); + secStMCPred->Scale(fHistDCA->GetSumOfWeights()*value/secStMCPred->GetSumOfWeights()); + secStMCPred->SetTitle("Sec from strangeness"); + secStMCPred->DrawCopy("same"); + if(partname=="PosP" && secMCPred){ + fit->GetResult(2,value,error); + secMCPred->Scale(fHistDCA->GetSumOfWeights()*value/secMCPred->GetSumOfWeights()); + secMCPred->SetTitle("Sec from material"); + secMCPred->DrawCopy("same"); + } + prim=PrimMCPred->Integral(PrimMCPred->FindBin(-xyMax),PrimMCPred->FindBin(xyMax)); + secSt=secStMCPred->Integral(secStMCPred->FindBin(-xyMax),secStMCPred->FindBin(xyMax)); + if(partname=="PosP")sec=secMCPred->Integral(secMCPred->FindBin(-xyMax),secMCPred->FindBin(xyMax)); + } + else{ + prim=1; + secSt=1; + if(partname=="PosP")sec=1; + } + Printf("Yields: primary=%f material=%f strange=%f\n",prim,sec,secSt); + return prim/(prim+secSt+sec); +} -- 2.43.0