From 4b8bdb60a48b52d84dcd6d8c5af0116688524b1f Mon Sep 17 00:00:00 2001 From: hdalsgaa Date: Wed, 18 Aug 2010 10:43:05 +0000 Subject: [PATCH] Scripts to make FMD analysis correction objects --- PWG2/FORWARD/analysis/AliFMDAnaParameters.cxx | 2 + .../AliFMDAnalysisTaskBFCorrelation.cxx | 394 ++++++++++++++++-- .../AliFMDAnalysisTaskBFCorrelation.h | 37 +- .../AliFMDAnalysisTaskGenerateBackground.cxx | 332 --------------- .../AliFMDAnalysisTaskGenerateBackground.h | 66 --- .../AliFMDAnalysisTaskGenerateCorrection.cxx | 10 +- PWG2/FORWARD/analysis/AliFMDDndeta.cxx | 8 +- PWG2/FORWARD/analysis/SubmitFMDCorrections.C | 47 +++ .../analysis/SubmitFMDEnergyDistributions.C | 49 +++ .../analysis/SubmitSharingEffCorrection.C | 55 +++ 10 files changed, 549 insertions(+), 451 deletions(-) delete mode 100644 PWG2/FORWARD/analysis/AliFMDAnalysisTaskGenerateBackground.cxx delete mode 100644 PWG2/FORWARD/analysis/AliFMDAnalysisTaskGenerateBackground.h create mode 100644 PWG2/FORWARD/analysis/SubmitFMDCorrections.C create mode 100644 PWG2/FORWARD/analysis/SubmitFMDEnergyDistributions.C create mode 100644 PWG2/FORWARD/analysis/SubmitSharingEffCorrection.C diff --git a/PWG2/FORWARD/analysis/AliFMDAnaParameters.cxx b/PWG2/FORWARD/analysis/AliFMDAnaParameters.cxx index 7c43376f2bf..e20ec78d81c 100644 --- a/PWG2/FORWARD/analysis/AliFMDAnaParameters.cxx +++ b/PWG2/FORWARD/analysis/AliFMDAnaParameters.cxx @@ -99,6 +99,8 @@ AliFMDAnaParameters::AliFMDAnaParameters() : // Default constructor fPhysicsSelection = new AliPhysicsSelection; fPhysicsSelection->SetAnalyzeMC(kTRUE); //For the background correction. This is reset in Init if relevant + // fPhysicsSelection->SetUseBXNumbers(kFALSE); + AliBackgroundSelection* backgroundSelection = new AliBackgroundSelection("bg","bg"); diff --git a/PWG2/FORWARD/analysis/AliFMDAnalysisTaskBFCorrelation.cxx b/PWG2/FORWARD/analysis/AliFMDAnalysisTaskBFCorrelation.cxx index 790321e0eeb..c6a8490236f 100644 --- a/PWG2/FORWARD/analysis/AliFMDAnalysisTaskBFCorrelation.cxx +++ b/PWG2/FORWARD/analysis/AliFMDAnalysisTaskBFCorrelation.cxx @@ -1,4 +1,3 @@ - #include #include #include @@ -6,8 +5,11 @@ #include #include #include +#include +#include #include "TH1F.h" #include "TH2F.h" +#include "TProfile.h" #include "AliFMDAnalysisTaskBFCorrelation.h" #include "AliAnalysisManager.h" #include "AliESDFMD.h" @@ -27,50 +29,188 @@ //#include "TDatabasePDG.h" //#include "TParticlePDG.h" #include "AliESDInputHandler.h" -ClassImp(AliFMDAnalysisTaskBFCorrelation) +ClassImp(AliFMDAnalysisTaskBFCorrelation) AliFMDAnalysisTaskBFCorrelation::AliFMDAnalysisTaskBFCorrelation() : fDebug(0), fOutputList(0), fInputList(0), + fInternalList(0), fVertexString(0x0), - fStandalone(kTRUE) + fStandalone(kTRUE), + fEvent(0), + fnBinsX(200), + fXmin(-6), + fXmax(6), + fnBinsY(20), + fYmin(0), + fYmax(2 * TMath::Pi()), + c(0), + debug0(0), + debug1(0) { // Default constructor DefineInput (0, TList::Class()); DefineOutput(0, TList::Class()); } + //_____________________________________________________________________ AliFMDAnalysisTaskBFCorrelation::AliFMDAnalysisTaskBFCorrelation(const char* name, Bool_t SE): AliAnalysisTask(name,name), - fDebug(0), - fOutputList(0), - fInputList(0), - fVertexString(0x0), - fStandalone(kTRUE) + fDebug(0), + fOutputList(0), + fInputList(0), + fInternalList(0), + fVertexString(0x0), + fStandalone(kTRUE), + fEvent(0), + fnBinsX(200), + fXmin(-6), + fXmax(6), + fnBinsY(20), + fYmin(0), + fYmax(2 * TMath::Pi()), + c(0), + debug0(0), + debug1(0) { fStandalone = SE; if(fStandalone) { DefineInput (0, TList::Class()); DefineInput(1, TObjString::Class()); DefineOutput(0, TList::Class()); - } } + //_____________________________________________________________________ void AliFMDAnalysisTaskBFCorrelation::CreateOutputObjects() { - //AliFMDAnaParameters* pars = AliFMDAnaParameters::Instance(); - + // Setup the list for storing results, if it does not exist + if(!fOutputList) { fOutputList = new TList(); fOutputList->SetName("BackgroundCorrected"); } + + // Setup the list for temporary storage during the analysis + + if(!fInternalList) { + fInternalList = new TList(); + fInternalList->SetName("InternalBFList"); + } + + // Set up histograms for analysis and storage. 4 different binnings - TH1F* test = new TH1F("test","test",10,0,10); - fOutputList->Add(test); + for (Int_t i = 1; i <= 4; i++) { + + // Temporary histograms for storing hist per eta summed over phi + + TH1D *hESDMult = new TH1D(Form("hESDMult_binning%d", i), + Form("Multiplicity per event vs eta ESD (%d bins)", 30*i), + 18*i, -9, 9); + TH1D *hMCMult = new TH1D(Form("hMCMult_binning%d", i), + Form("Multiplicity per event vs eta MC-truth (%d bins)", 30*i), + 18*i, -9, 9); + fInternalList->Add(hESDMult); + fInternalList->Add(hMCMult); + + + // Histograms for storing the 5 values to calculate b (ESD & MC) + + TH1D *pnFnB = new TH1D(Form("pnFnB_binning%d", i), + Form("Correlation ESD (%d bins)", i*9), + 9*i, 0, 9); + TH1D *pnF2 = new TH1D(Form("pnF2_binning%d", i), + Form("Fwd^2 ESD (%d bins)", i*9), + 9*i, 0, 9); + TH1D *pnB2 = new TH1D(Form("pnB2_binning%d", i), + Form("Bwd^2 ESD (%d bins)", i*9), + 9*i, 0, 9); + TH1D *pnF = new TH1D(Form("pnF_binning%d", i), + Form("Fwd ESD (%d bins)", i*9), + 9*i, 0, 9); + TH1D *pnB = new TH1D(Form("pnB_binning%d", i), + Form("Bwd ESD (%d bins)", i*9), + 9*i, 0, 9); + + TH1D *pnFnB_MC = new TH1D(Form("pnFnB_MC_binning%d", i), + Form("Correlation MC (%d bins)", i*9), + 9*i, 0, 9); + TH1D *pnF2_MC = new TH1D(Form("pnF2_MC_binning%d", i), + Form("Fwd^2 MC (%d bins)", i*9), + 9*i, 0, 9); + TH1D *pnB2_MC = new TH1D(Form("pnB2_MC_binning%d", i), + Form("Bwd^2 MC (%d bins)", i*9), + 9*i, 0, 9); + TH1D *pnF_MC = new TH1D(Form("pnF_MC_binning%d", i), + Form("Fwd MC (%d bins)", i*9), + 9*i, 0, 9); + TH1D *pnB_MC = new TH1D(Form("pnB_MC_binning%d", i), + Form("Backward MC (%d bins)" ,i*9), + 9*i, 0, 9); + fOutputList->Add(pnFnB); + fOutputList->Add(pnF2); + fOutputList->Add(pnB2); + fOutputList->Add(pnF); + fOutputList->Add(pnB); + + fOutputList->Add(pnFnB_MC); + fOutputList->Add(pnF2_MC); + fOutputList->Add(pnB2_MC); + fOutputList->Add(pnF_MC); + fOutputList->Add(pnB_MC); + + // Debugging histograms : "dNdEta" + + // TH1D *hESDMultvsEta = new TH1D(Form("hESDMultvsEta_binning%d" ,i), + // Form("Mult vs Eta (%d bins, ESD)", 18*i), + // 18*i, -9, 9); + TH1D *hESDMultvsEta = new TH1D(Form("hESDMultvsEta_binning%d" ,i), + Form("Mult vs Eta (%d bins, ESD)", i*18), + 200, -6, 6); + TH1D *hMCMultvsEta = new TH1D(Form("hMCMultvsEta_binning%d", i), + Form("Mult vs Eta (%d bins, MC)", 18*i), + 200, -6, 6); + //TH1D *hMCMultvsEta = new TH1D(Form("hMCMultvsEta_binning%d", i), + // Form("Mult vs Eta (%d bins, MC)", 18*i), + // 18*i, -9, 9); + fOutputList->Add(hESDMultvsEta); + fOutputList->Add(hMCMultvsEta); + + TH1D *hESDMultvsEtaC = new TH1D(Form("hESDMultvsEtaC_binning%d" ,i), + Form("Mult vs Eta C (%d bins, ESD)", 18*i), + 18*i, -9, 9); + TH1D *hMCMultvsEtaC = new TH1D(Form("hMCMultvsEtaC_binning%d", i), + Form("Mult vs Eta C (%d bins, MC)", 18*i), + 18*i, -9, 9); + fOutputList->Add(hESDMultvsEtaC); + fOutputList->Add(hMCMultvsEtaC); + + // Debugging histograms : Different distributions + + TH2D *hESDBwdDist = new TH2D(Form("hESDBwdDist_binning%d", i), + Form("Distribution of Bwd values (%d bins)", i*9), + 9*i, 0, 9, 20, 0, 20); + + TH2D *hESDBwd2Dist = new TH2D(Form("hESDBwd2Dist_binning%d", i), + Form("Distribution of Bwd^2 values (%d bins)", i*9), + 9*i, 0, 9, 400, 0, 400); + + TH2D *hMCBwdDist = new TH2D(Form("hMCBwdDist_binning%d", i), + Form("Distribution of Bwd values (%d bins)", i*9), + 9*i, 0, 9, 20, 0, 20); + + TH2D *hMCBwd2Dist = new TH2D(Form("hMCBwd2Dist_binning%d", i), + Form("Distribution of Bwd^2 values (%d bins)", i*9), + 9*i, 0, 9, 400, 0, 400); + fOutputList->Add(hESDBwdDist); + fOutputList->Add(hESDBwd2Dist); + fOutputList->Add(hMCBwdDist); + fOutputList->Add(hMCBwd2Dist); + } } + //_____________________________________________________________________ void AliFMDAnalysisTaskBFCorrelation::ConnectInputData(Option_t */*option*/) { @@ -79,45 +219,187 @@ void AliFMDAnalysisTaskBFCorrelation::ConnectInputData(Option_t */*option*/) fVertexString = (TObjString*)GetInputData(1); } } + //_____________________________________________________________________ void AliFMDAnalysisTaskBFCorrelation::Exec(Option_t */*option*/) { + fEvent++; + if (fEvent % 1000 == 0) + std::cout << "Event # " << fEvent << std::endl; + AliFMDAnaParameters* pars = AliFMDAnaParameters::Instance(); fVertexString = (TObjString*)fInputList->At(0); - Int_t vtxbin = fVertexString->GetString().Atoi(); - //for(UShort_t det=1;det<=3;det++) { - // Int_t nRings = (det==1 ? 1 : 2); - // for (UShort_t ir = 0; ir < nRings; ir++) { - // Char_t ringChar = (ir == 0 ? 'I' : 'O'); - //TH2F* hMultTotal = (TH2F*)fOutputList->FindObject(Form("dNdeta_FMD%d%c_vtxbin%d",det,ringChar,vtxbin)); - // TH2F* hMultTotalTrVtx = (TH2F*)fOutputList->FindObject(Form("dNdetaTrVtx_FMD%d%c_vtxbin%d",det,ringChar,vtxbin)); - + if(vtxbin !=4) return; + SetBounds(); + CountESDHits(); + CalculateValues("ESD"); - TH1F* test = (TH1F*)fOutputList->FindObject("test"); - test->Fill(vtxbin); if(pars->GetProcessPrimary()) ProcessPrimary(); if(fStandalone) { PostData(0, fOutputList); } +} +//_____________________________________________________________________ +void AliFMDAnalysisTaskBFCorrelation::CountESDHits() { + + TH2F *dNdEtadPhiHist = (TH2F*)fInputList->FindObject("dNdetadphiHistogramSPDTrVtx"); + + for (Int_t i = 1; i<=4; i++) { + + TH1D *hESDMult = (TH1D*)fInternalList->FindObject(Form("hESDMult_binning%d", i)); + hESDMult->Reset(); + } + + Int_t etamax = dNdEtadPhiHist->GetNbinsX(); + Int_t phimax = dNdEtadPhiHist->GetNbinsY(); + + TH1D *hprod = dNdEtadPhiHist->ProjectionX(); + + TH1D *hESDMultvsEta = (TH1D*)fOutputList->FindObject(Form("hESDMultvsEta_binning%d" ,4)); + hESDMultvsEta->Add(hprod); + + for (Int_t etabin = 1; etabin <= etamax; etabin++) { + Float_t val = 0; + Float_t eta = dNdEtadPhiHist->GetXaxis()->GetBinCenter(etabin); + + + for (Int_t i = 1; i <= 4; i++) { + TH1D *hESDMult = (TH1D*)fInternalList->FindObject(Form("hESDMult_binning%d", i)); + hESDMult->Fill(eta, hprod->GetBinContent(etabin)); + } + } + + //hESDMult->Add(hprod); } //_____________________________________________________________________ -void AliFMDAnalysisTaskBFCorrelation::Terminate(Option_t */*option*/) { +void AliFMDAnalysisTaskBFCorrelation::CalculateValues(TString type) { + + type.ToUpper(); + + const char *stype = type.Data(); + + TString sinput; + TString soutput; + TString soutputC; + TString snFnB; + TString snF2; + TString snB2; + TString snF; + TString snB; + + if (type == "ESD") { + snFnB = "pnFnB"; + snF2 = "pnF2"; + snB2 = "pnB2"; + snF = "pnF"; + snB = "pnB"; + } else { + snFnB = "pnFnB_MC"; + snF2 = "pnF2_MC"; + snB2 = "pnB2_MC"; + snF = "pnF_MC"; + snB = "pnB_MC"; + } + + for (Int_t i = 1; i <= 4; i++) { + + sinput = Form("h%sMult_binning%d", stype, i); + soutput = Form("h%sMultvsEta_binning%d", stype, i); + soutputC = Form("h%sMultvsEtaC_binning%d", stype, i); + + TH1D *input = (TH1D*)fInternalList->FindObject(sinput); + TH1D *output = (TH1D*)fOutputList->FindObject(soutput); + TH1D *outputC = (TH1D*)fOutputList->FindObject(soutputC); + + //output->Add(input); + + for (Int_t bin = 1; bin <= input->GetNbinsX()/2; bin++) { + + Double_t eta = input->GetBinCenter(bin); + Double_t bwd = input->GetBinContent(bin); + Double_t bwd2 = TMath::Power(bwd, 2); + Double_t fwd = input->GetBinContent(input->FindBin(TMath::Abs(eta))); + Double_t fwd2 = TMath::Power(fwd, 2); + Double_t fwdbwd = fwd*bwd; + + snFnB += Form("_binning%d", i); + snF2 += Form("_binning%d", i); + snB2 += Form("_binning%d", i); + snF += Form("_binning%d", i); + snB += Form("_binning%d", i); + + TH1D *nFnB = static_cast(fOutputList->FindObject(snFnB)); + TH1D *nF2 = static_cast(fOutputList->FindObject(snF2)); + TH1D *nB2 = static_cast(fOutputList->FindObject(snB2)); + TH1D *nF = static_cast(fOutputList->FindObject(snF)); + TH1D *nB = static_cast(fOutputList->FindObject(snB)); + + nFnB->Fill(TMath::Abs(eta), fwdbwd); + nF2->Fill(TMath::Abs(eta), fwd2); + nB2->Fill(TMath::Abs(eta), bwd2); + nF->Fill(TMath::Abs(eta), fwd); + nB->Fill(TMath::Abs(eta), bwd); + + outputC->Fill(TMath::Abs(eta), fwd); + outputC->Fill(eta, bwd); + + if (type == "ESD") { + + snFnB.Remove(5,9); + snF2.Remove(4,9); + snB2.Remove(4,9); + snF.Remove(3,9); + snB.Remove(3,9); + TH2D *hESDBwdDist = (TH2D*)fOutputList->FindObject(Form("hESDBwdDist_binning%d", i)); + TH2D *hESDBwd2Dist = (TH2D*)fOutputList->FindObject(Form("hESDBwd2Dist_binning%d", i)); + hESDBwdDist->Fill(TMath::Abs(eta), bwd); + hESDBwd2Dist->Fill(TMath::Abs(eta), bwd2); + + } else { + + snFnB.Remove(8,9); + snF2.Remove(7,9); + snB2.Remove(7,9); + snF.Remove(6,9); + snB.Remove(6,9); + TH2D *hMCBwdDist = (TH2D*)fOutputList->FindObject(Form("hMCBwdDist_binning%d", i)); + TH2D *hMCBwd2Dist = (TH2D*)fOutputList->FindObject(Form("hMCBwd2Dist_binning%d", i)); + hMCBwdDist->Fill(TMath::Abs(eta), bwd); + hMCBwd2Dist->Fill(TMath::Abs(eta), bwd2); + } + } + } +} + +//_____________________________________________________________________ +void AliFMDAnalysisTaskBFCorrelation::SetBounds() { + + TH2D *hTemp = (TH2D*)fInputList->FindObject("multTrVtx_FMD1I_vtxbin0"); + fnBinsX = hTemp->GetNbinsX(); + fnBinsY = hTemp->GetNbinsY(); + fXmin = hTemp->GetXaxis()->GetXmin(); + fXmax = hTemp->GetXaxis()->GetXmax(); + fYmin = hTemp->GetYaxis()->GetXmin(); + fYmax = hTemp->GetYaxis()->GetXmax(); +} + +//_____________________________________________________________________ +void AliFMDAnalysisTaskBFCorrelation::Terminate(Option_t */*option*/) { } //_____________________________________________________________________ void AliFMDAnalysisTaskBFCorrelation::ProcessPrimary() { - /* + AliMCEventHandler* eventHandler = dynamic_cast (AliAnalysisManager::GetAnalysisManager()->GetMCtruthEventHandler()); AliMCEvent* mcEvent = eventHandler->MCEvent(); if(!mcEvent) return; - AliFMDAnaParameters* pars = AliFMDAnaParameters::Instance(); @@ -128,39 +410,61 @@ void AliFMDAnalysisTaskBFCorrelation::ProcessPrimary() { AliHeader* header = mcEvent->Header(); AliGenEventHeader* genHeader = header->GenEventHeader(); + AliGenPythiaEventHeader* pythiaGenHeader = dynamic_cast(genHeader); + + if (!pythiaGenHeader) { + std::cout<<" no pythia header!"<ProcessType(); + + if(pythiaType==92||pythiaType==93){ + std::cout<<"single diffractive"<PrimaryVertex(vertex); if(TMath::Abs(vertex.At(2)) > pars->GetVtxCutZ()) return; - Double_t delta = 2*pars->GetVtxCutZ()/pars->GetNvtxBins(); - Double_t vertexBinDouble = (vertex.At(2) + pars->GetVtxCutZ()) / delta; - Int_t vertexBin = (Int_t)vertexBinDouble; + //Double_t delta = 2*pars->GetVtxCutZ()/pars->GetNvtxBins(); + //Double_t vertexBinDouble = (vertex.At(2) + pars->GetVtxCutZ()) / delta; + //Int_t vertexBin = (Int_t)vertexBinDouble; - Bool_t firstTrack = kTRUE; + //Bool_t firstTrack = kTRUE; // we loop over the primaries only unless we need the hits (diagnostics running slowly) Int_t nTracks = stack->GetNprimary(); - if(pars->GetProcessHits()) - nTracks = stack->GetNtrack(); + // if(pars->GetProcessHits()) + // nTracks = stack->GetNtrack(); + TH1D *hMCMultvsEta = (TH1D*)fOutputList->FindObject(Form("hMCMultvsEta_binning%d" ,4)); - for(Int_t i = 0 ;iGetTrack(i); - if(!particle) + for (Int_t i = 1; i <= 4; i++) { + + TH1D *hMCMult = (TH1D*)fInternalList->FindObject(Form("hMCMult_binning%d", i)); + hMCMult->Reset(); + + for(Int_t ii = 0 ;iiGetTrack(ii); + if(!particle) continue; - - if(stack->IsPhysicalPrimary(i) && particle->Charge() != 0) { - hPrimary->Fill(particle->Eta()); - - TH1F* hPrimVtxBin = (TH1F*)fOutputList->FindObject(Form("primmult_vtxbin%d",vertexBin)); - hPrimVtxBin->Fill(particle->Eta()); - if(firstTrack) { - fNMCevents.Fill(vertexBin); - firstTrack = kFALSE; + if(stack->IsPhysicalPrimary(ii) && particle->Charge() != 0) { + if (i == 1) + hPrimary->Fill(particle->Eta()); + if(i==4) + hMCMultvsEta->Fill(particle->Eta()); + hMCMult->Fill(particle->Eta()); } } } - */ + CalculateValues("MC"); } //_____________________________________________________________________ // diff --git a/PWG2/FORWARD/analysis/AliFMDAnalysisTaskBFCorrelation.h b/PWG2/FORWARD/analysis/AliFMDAnalysisTaskBFCorrelation.h index e429d9eaea6..3ffc86cce5c 100644 --- a/PWG2/FORWARD/analysis/AliFMDAnalysisTaskBFCorrelation.h +++ b/PWG2/FORWARD/analysis/AliFMDAnalysisTaskBFCorrelation.h @@ -10,8 +10,10 @@ #include "TObjString.h" #include "TArrayI.h" #include "TH1I.h" +#include "TH2.h" #include "AliMCEvent.h" #include "AliFMDFloatMap.h" +#include "TCanvas.h" /** * @ingroup FMD_ana @@ -26,9 +28,20 @@ class AliFMDAnalysisTaskBFCorrelation : public AliAnalysisTask fDebug(o.fDebug), fOutputList(0), fInputList(0), + fInternalList(0), fVertexString(o.fVertexString), - fStandalone(o.fStandalone) - {} + fStandalone(o.fStandalone), + fEvent(0), + fnBinsX(200), + fXmin(-6), + fXmax(6), + fnBinsY(20), + fYmin(2), + fYmax(2 * TMath::Pi()), + c(0), + debug0(0), + debug1(0) + {} AliFMDAnalysisTaskBFCorrelation& operator=(const AliFMDAnalysisTaskBFCorrelation&) { return *this; } // Implementation of interface methods @@ -42,17 +55,37 @@ class AliFMDAnalysisTaskBFCorrelation : public AliAnalysisTask void SetInputList(TList* inputList) {fInputList = inputList;} void SetInputVertex(TObjString* vtxString) {fVertexString = vtxString;} void SetOutputList(TList* outputList) {fOutputList = outputList;} + void CountESDHits(); + void CalculateValues(TString type); + void SetBounds(); + void MergeEvent(TH2D *hMultTrVtxFull); void ProcessPrimary(); TList* GetOutputList() {return fOutputList;} private: + Int_t fDebug; // Debug flag TList* fOutputList; TList* fInputList; + TList* fInternalList; TObjString* fVertexString; Bool_t fStandalone; + + Int_t fEvent; + Int_t fnBinsX; + Float_t fXmin; + Float_t fXmax; + Int_t fnBinsY; + Float_t fYmin; + Float_t fYmax; + + TCanvas *c; + + Double_t debug0; + Double_t debug1; + ClassDef(AliFMDAnalysisTaskBFCorrelation, 0); // Analysis task for FMD analysis }; diff --git a/PWG2/FORWARD/analysis/AliFMDAnalysisTaskGenerateBackground.cxx b/PWG2/FORWARD/analysis/AliFMDAnalysisTaskGenerateBackground.cxx deleted file mode 100644 index 35e617888d9..00000000000 --- a/PWG2/FORWARD/analysis/AliFMDAnalysisTaskGenerateBackground.cxx +++ /dev/null @@ -1,332 +0,0 @@ -#include "AliFMDAnalysisTaskGenerateBackground.h" -#include "AliESDEvent.h" -#include "iostream" -#include "AliESDFMD.h" -#include "TH2F.h" -#include "AliTrackReference.h" -#include "AliStack.h" -#include "AliFMDAnaParameters.h" -#include "AliFMDStripIndex.h" -#include "AliStack.h" -#include "AliMCParticle.h" -#include "AliMCEvent.h" -//#include "AliFMDGeometry.h" -#include "TArray.h" -#include "AliGenEventHeader.h" -#include "AliHeader.h" -#include "AliFMDAnaCalibBackgroundCorrection.h" -//#include "AliCDBManager.h" -//#include "AliCDBId.h" -//#include "AliCDBMetaData.h" -#include "TSystem.h" -#include "TROOT.h" -#include "TAxis.h" -ClassImp(AliFMDAnalysisTaskGenerateBackground) - -//_____________________________________________________________________ -AliFMDAnalysisTaskGenerateBackground::AliFMDAnalysisTaskGenerateBackground(): -AliAnalysisTaskSE(), - fListOfHits(), - fListOfPrimaries(), - fListOfCorrection(), - fVertexBins(), - fLastTrackByStrip(0), - fHitsByStrip(0), - fZvtxCut(10), - fNvtxBins(10), - fNbinsEta(200), - fBackground(0) -{ - // Default constructor -} -//_____________________________________________________________________ -AliFMDAnalysisTaskGenerateBackground::AliFMDAnalysisTaskGenerateBackground(const char* name): - AliAnalysisTaskSE(name), - fListOfHits(), - fListOfPrimaries(), - fListOfCorrection(), - fVertexBins(), - fLastTrackByStrip(0), - fHitsByStrip(0), - fZvtxCut(10), - fNvtxBins(10), - fNbinsEta(200), - fBackground(0) -{ - - DefineOutput(1, TList::Class()); - DefineOutput(2, TList::Class()); - DefineOutput(3, TH1F::Class()); - DefineOutput(4, TList::Class()); -} -//_____________________________________________________________________ -void AliFMDAnalysisTaskGenerateBackground::UserCreateOutputObjects() -{ -// Create the output containers -// - - std::cout<<"Creating output objects"<Sumw2(); - fListOfPrimaries.Add(hPrimary); - } - } - - - for(Int_t det =1; det<=3;det++) { - Int_t nRings = (det==1 ? 1 : 2); - for(Int_t ring = 0;ringSumw2(); - allHits->Sumw2(); - fListOfHits.Add(allHits); - fListOfHits.Add(doubleHits); - - for(Int_t v=0; vSumw2(); - fListOfHits.Add(hHits); - - } - } - } - - fVertexBins.SetName("VertexBins"); - fVertexBins.GetXaxis()->Set(fNvtxBins,-1*fZvtxCut,fZvtxCut); - -} -//_____________________________________________________________________ -void AliFMDAnalysisTaskGenerateBackground::Init() -{ - fLastTrackByStrip.Reset(-1); - - -} -//_____________________________________________________________________ -void AliFMDAnalysisTaskGenerateBackground::UserExec(Option_t */*option*/) -{ - - fLastTrackByStrip.Reset(-1); - fHitsByStrip.Reset(0); - AliMCEvent* mcevent = MCEvent(); - AliFMDAnaParameters* pars = AliFMDAnaParameters::Instance(); - - AliMCParticle* particle = 0; - AliStack* stack = mcevent->Stack(); - - UShort_t det,sec,strip; - Char_t ring; - - Int_t nTracks = mcevent->GetNumberOfTracks(); - AliHeader* header = mcevent->Header(); - AliGenEventHeader* genHeader = header->GenEventHeader(); - - TArrayF vertex; - genHeader->PrimaryVertex(vertex); - - if(TMath::Abs(vertex.At(2)) > fZvtxCut) - return; - for(Int_t i = 0 ;iGetTrack(i); - - if(!particle) - continue; - - Double_t delta = 2*fZvtxCut/fNvtxBins; - Double_t vertexBinDouble = (vertex.At(2) + fZvtxCut) / delta; - Int_t vertexBin = (Int_t)vertexBinDouble; - - if(stack->IsPhysicalPrimary(i) && particle->Charge() != 0) { - - - TH2F* hPrimaryInner = (TH2F*)fListOfPrimaries.FindObject( Form("hPrimary_FMD_%c_vtx%d",'I',vertexBin)); - TH2F* hPrimaryOuter = (TH2F*)fListOfPrimaries.FindObject( Form("hPrimary_FMD_%c_vtx%d",'O',vertexBin)); - hPrimaryInner->Fill(particle->Eta(),particle->Phi()); - hPrimaryOuter->Fill(particle->Eta(),particle->Phi()); - } - - for(Int_t j=0; jGetNumberOfTrackReferences();j++) { - - AliTrackReference* ref = particle->GetTrackReference(j); - - if(ref->DetectorId() != AliTrackReference::kFMD) - continue; - AliFMDStripIndex::Unpack(ref->UserId(),det,ring,sec,strip); - Float_t thisStripTrack = fLastTrackByStrip(det,ring,sec,strip); - if(particle->Charge() != 0 && i != thisStripTrack ) { - //Double_t x,y,z; - //AliFMDGeometry* fmdgeo = AliFMDGeometry::Instance(); - //fmdgeo->Detector2XYZ(det,ring,sec,strip,x,y,z); - Float_t phi = pars->GetPhiFromSector(det,ring,sec); - Float_t eta = pars->GetEtaFromStrip(det,ring,sec,strip,vertex.At(2)); - //Float_t phi = TMath::ATan2(y,x); - //if(phi<0) phi = phi+2*TMath::Pi(); - // Float_t r = TMath::Sqrt(TMath::Power(x,2)+TMath::Power(y,2)); - //Float_t theta = TMath::ATan2(r,z-vertex.At(2)); - //Float_t eta = -1*TMath::Log(TMath::Tan(0.5*theta)); - TH2F* hHits = (TH2F*)fListOfHits.FindObject(Form("hHits_FMD%d%c_vtx%d", det,ring,vertexBin)); - hHits->Fill(eta,phi); - Float_t nstrips = (ring =='O' ? 256 : 512); - fHitsByStrip(det,ring,sec,strip) +=1; - TH1F* allHits = (TH1F*)fListOfHits.FindObject(Form("allHits_FMD%d%c",det,ring)); - TH1F* doubleHits = (TH1F*)fListOfHits.FindObject(Form("DoubleHits_FMD%d%c",det,ring)); - - if(fHitsByStrip(det,ring,sec,strip) == 1) - allHits->Fill(eta); - - doubleHits->Fill(eta); - /*if(fHitsByStrip.operator()(det,ring,sec,strip) == 2){ - TH1F* doubleHits = (TH1F*)fListOfHits.FindObject("DoubleHits"); - doubleHits->Fill(eta,2); - }*/ - //if(fHitsByStrip.operator()(det,ring,sec,strip) > 1){ - // doubleHits->Fill(eta); - // } - - fLastTrackByStrip(det,ring,sec,strip) = (Float_t)i; - if(strip >0) - fLastTrackByStrip(det,ring,sec,strip-1) = (Float_t)i; - if(strip < (nstrips - 1)) - fLastTrackByStrip(det,ring,sec,strip+1) = (Float_t)i; - } - } - - } - - - PostData(1, &fListOfHits); - PostData(2, &fListOfPrimaries); - PostData(3, &fVertexBins); -} -//_____________________________________________________________________ -void AliFMDAnalysisTaskGenerateBackground::Terminate(Option_t */*option*/) -{ - /* TH1F* allHits = (TH1F*)fListOfHits.FindObject("allHits"); - TH1F* doubleHits = (TH1F*)fListOfHits.FindObject("DoubleHits"); - - doubleHits->Divide(allHits); - GenerateCorrection(); - PostData(1, &fListOfHits); - PostData(4, &fListOfCorrection);*/ - -} -//_____________________________________________________________________ -void AliFMDAnalysisTaskGenerateBackground::GenerateCorrection() { - - fBackground = new AliFMDAnaCalibBackgroundCorrection(); - - for(Int_t det= 1; det <=3; det++) { - Int_t nRings = (det==1 ? 1 : 2); - - for(Int_t iring = 0; iringSetDoubleHitCorrection(det,ring,doubleHits); - doubleHits->Divide(allHits); - for(Int_t vertexBin=0;vertexBinClone(Form("FMD%d%c_vtxbin_%d_correction",det,ring,vertexBin)); - hCorrection->Divide(hPrimary); - hCorrection->SetTitle(hCorrection->GetName()); - fListOfCorrection.Add(hCorrection); - fBackground->SetBgCorrection(det,ring,vertexBin,hCorrection); - - - } - - } - } - TAxis refAxis(fNvtxBins,-1*fZvtxCut,fZvtxCut); - fBackground->SetRefAxis(&refAxis); - -} -//_____________________________________________________________________ -void AliFMDAnalysisTaskGenerateBackground::ReadFromFile(const Char_t* filename, Bool_t storeInOCDB, Int_t /*runNo*/) { - - TFile infile(filename); - TH1F* hVertex = (TH1F*)infile.Get("VertexBins"); - fZvtxCut = hVertex->GetXaxis()->GetXmax(); - fNvtxBins = hVertex->GetXaxis()->GetNbins(); - fVertexBins.SetName("VertexBins"); - fVertexBins.GetXaxis()->Set(fNvtxBins,-1*fZvtxCut,fZvtxCut); - - TList* listOfHits = (TList*)infile.Get("Hits"); - TList* listOfPrim = (TList*)infile.Get("Primaries"); - - for(Int_t det =1; det<=3;det++) - { - Int_t nRings = (det==1 ? 1 : 2); - for(Int_t ring = 0;ringFindObject(Form("allHits_FMD%d%c",det,ringChar)); - TH1F* doubleHits = (TH1F*)listOfHits->FindObject(Form("DoubleHits_FMD%d%c",det,ringChar)); - fListOfHits.Add(allHits); - fListOfHits.Add(doubleHits); - for(Int_t v=0; vFindObject(Form("hHits_FMD%d%c_vtx%d", det,ringChar,v)); - fListOfHits.Add(hHits); - } - } - } - for(Int_t iring = 0; iring<2;iring++) { - Char_t ringChar = (iring == 0 ? 'I' : 'O'); - for(Int_t v=0; vFindObject( Form("hPrimary_FMD_%c_vtx%d",ringChar,v)); - fListOfPrimaries.Add(hPrimary); - - } - } - GenerateCorrection(); - - TFile fout("backgroundFromFile.root","recreate"); - fListOfHits.Write(); - fListOfPrimaries.Write(); - fListOfCorrection.Write(); - fVertexBins.Write(); - - fout.Close(); - - if(storeInOCDB) { - TFile fcalib("$ALICE_ROOT/FMD/Correction/Background/background.root","RECREATE"); - fBackground->Write(AliFMDAnaParameters::GetBackgroundID()); - fcalib.Close(); - /* AliCDBManager* cdb = AliCDBManager::Instance(); - cdb->SetDefaultStorage("local://$ALICE_ROOT/OCDB"); - AliCDBId id(AliFMDAnaParameters::GetBackgroundPath(),runNo,999999999); - - AliCDBMetaData* meta = new AliCDBMetaData; - meta->SetResponsible(gSystem->GetUserInfo()->fRealName.Data()); - meta->SetAliRootVersion(gROOT->GetVersion()); - meta->SetBeamPeriod(1); - meta->SetComment("Background Correction for FMD"); - meta->SetProperty("key1", fBackground ); - cdb->Put(fBackground, id, meta); - */ - } - -} -//_____________________________________________________________________ -// -// EOF -// diff --git a/PWG2/FORWARD/analysis/AliFMDAnalysisTaskGenerateBackground.h b/PWG2/FORWARD/analysis/AliFMDAnalysisTaskGenerateBackground.h deleted file mode 100644 index 4de2bbd03cb..00000000000 --- a/PWG2/FORWARD/analysis/AliFMDAnalysisTaskGenerateBackground.h +++ /dev/null @@ -1,66 +0,0 @@ -#ifndef ALIFMDANALYSISTASKGENERATEBACKGROUND_H -#define ALIFMDANALYSISTASKGENERATEBACKGROUND_H - -#include "AliAnalysisTaskSE.h" -#include "TList.h" -#include "AliFMDFloatMap.h" -#include "TH1F.h" - -/** - * Make a background distribution from simulated data - * @ingroup FMD_ana - * - * - */ - -class AliFMDAnaCalibBackgroundCorrection; - -class AliFMDAnalysisTaskGenerateBackground : public AliAnalysisTaskSE -{ - public: - AliFMDAnalysisTaskGenerateBackground(); - AliFMDAnalysisTaskGenerateBackground(const char* name); - ~AliFMDAnalysisTaskGenerateBackground() {;} - AliFMDAnalysisTaskGenerateBackground(const AliFMDAnalysisTaskGenerateBackground& o) : AliAnalysisTaskSE(), - fListOfHits(), - fListOfPrimaries(), - fListOfCorrection(), - fVertexBins(o.fVertexBins), - fLastTrackByStrip(o.fLastTrackByStrip), - fHitsByStrip(o.fHitsByStrip), - fZvtxCut(o.fZvtxCut), - fNvtxBins(o.fNvtxBins), - fNbinsEta(o.fNbinsEta), - fBackground(o.fBackground) - {} - AliFMDAnalysisTaskGenerateBackground& operator=(const AliFMDAnalysisTaskGenerateBackground&) { return *this; } - - virtual void Init(); - virtual void UserCreateOutputObjects(); - virtual void UserExec(Option_t* /*option*/); - void Terminate(Option_t */*option*/); - void SetZvtxCut(Float_t vtxcut) {fZvtxCut = vtxcut;} - void SetNvtxBins(Int_t nvtxbins) {fNvtxBins = nvtxbins;} - void SetNbinsEta(Int_t netabins) {fNbinsEta = netabins;} - void ReadFromFile(const Char_t* filename = "background.root", Bool_t storeInOCDB = kFALSE, Int_t runNo=0); - private: - - void GenerateCorrection(); - - TList fListOfHits; - TList fListOfPrimaries; - TList fListOfCorrection; - TH1F fVertexBins; - AliFMDFloatMap fLastTrackByStrip; - AliFMDFloatMap fHitsByStrip; - Float_t fZvtxCut; - Int_t fNvtxBins; - Int_t fNbinsEta; - AliFMDAnaCalibBackgroundCorrection* fBackground; - ClassDef(AliFMDAnalysisTaskGenerateBackground, 1); - -}; -#endif -// Local Variables: -// mode: C++ -// End: diff --git a/PWG2/FORWARD/analysis/AliFMDAnalysisTaskGenerateCorrection.cxx b/PWG2/FORWARD/analysis/AliFMDAnalysisTaskGenerateCorrection.cxx index a04038c11a5..ab49f67a02f 100644 --- a/PWG2/FORWARD/analysis/AliFMDAnalysisTaskGenerateCorrection.cxx +++ b/PWG2/FORWARD/analysis/AliFMDAnalysisTaskGenerateCorrection.cxx @@ -277,10 +277,10 @@ void AliFMDAnalysisTaskGenerateCorrection::UserExec(Option_t */*option*/) if(!vtxStatus) vtxFound = kFALSE; - if(TMath::Abs(vertex.At(2)) > fZvtxCut) { - vtxFound = kFALSE; + //if(TMath::Abs(vertex.At(2)) > fZvtxCut) { + // vtxFound = kFALSE; - } + //} Bool_t isTriggered = pars->IsEventTriggered(AliFMDAnaParameters::kMB1); Bool_t isTriggeredNSD = pars->IsEventTriggered(AliFMDAnaParameters::kNSD); if(vtxFound && isTriggered) hEventsSelected->Fill(vertexBin); @@ -297,6 +297,10 @@ void AliFMDAnalysisTaskGenerateCorrection::UserExec(Option_t */*option*/) // if(!vtxFound || !isTriggered) return; + if(TMath::Abs(vertex.At(2)) > fZvtxCut) { + return; + } + for(Int_t i = 0 ;iGetTrack(i); diff --git a/PWG2/FORWARD/analysis/AliFMDDndeta.cxx b/PWG2/FORWARD/analysis/AliFMDDndeta.cxx index 8c41549cd7f..462165f00d9 100644 --- a/PWG2/FORWARD/analysis/AliFMDDndeta.cxx +++ b/PWG2/FORWARD/analysis/AliFMDDndeta.cxx @@ -376,9 +376,9 @@ void AliFMDDndeta::GenerateMult(Analysis what) { nNonZero++; } Int_t nBinsOld = fNbinsToCut; - if(det == 1 && ringChar =='I') { - fNbinsToCut = 0; - } + //if(det == 1 && ringChar =='I') { + // fNbinsToCut = 0; + // } TH1F* hRingMult = (TH1F*)fMultList[what]->FindObject(Form("hRingMult_FMD%d%c_%s",det,ringChar,fAnalysisNames[what])); for(Int_t i=1; i<=hRingMult->GetNbinsX(); i++) { @@ -1013,6 +1013,8 @@ void AliFMDDndeta::DrawDndeta(Analysis what, Int_t rebin, Bool_t realdata) { else if(pars->GetEnergy() == AliFMDAnaParameters::k900) fpyt = TFile::Open("/home/canute/ALICE/FMDanalysis/FirstAnalysis/pythia_study/pythiahists900.root","READ"); + + if(realdata ) { if(fpyt) { hPythiaMB = (TH1F*)fpyt->Get("hPythiaMB"); diff --git a/PWG2/FORWARD/analysis/SubmitFMDCorrections.C b/PWG2/FORWARD/analysis/SubmitFMDCorrections.C new file mode 100644 index 00000000000..66e43dc9330 --- /dev/null +++ b/PWG2/FORWARD/analysis/SubmitFMDCorrections.C @@ -0,0 +1,47 @@ +void SubmitFMDCorrections(const Char_t* filename, Bool_t store, Float_t energy, Int_t trigger, Float_t mag, Int_t collsystem) { + + gSystem->Load("libANALYSIS"); + gSystem->Load("libANALYSISalice"); + gSystem->Load("libPWG0base"); + gSystem->Load("libPWG2forward"); + + AliFMDAnaParameters* pars = AliFMDAnaParameters::Instance(); + + if(energy == 900) + pars->SetEnergy(AliFMDAnaParameters::k900); + else if(energy == 7000) + pars->SetEnergy(AliFMDAnaParameters::k7000); + else if(energy == 10000) + pars->SetEnergy(AliFMDAnaParameters::k10000); + else if(energy == 14000) + pars->SetEnergy(AliFMDAnaParameters::k14000); + else if(energy == 5500) + pars->SetEnergy(AliFMDAnaParameters::k5500); + + if(trigger == 0) + pars->SetTriggerDefinition(AliFMDAnaParameters::kMB1); + else if(trigger == 1) + pars->SetTriggerDefinition(AliFMDAnaParameters::kMB2); + + if(mag==0) + pars->SetMagField(AliFMDAnaParameters::k0G); + else if(mag==1) + pars->SetMagField(AliFMDAnaParameters::k5G); + + if(collsystem == 0) + pars->SetCollisionSystem(AliFMDAnaParameters::kPP); + else if(collsystem == 1) + pars->SetCollisionSystem(AliFMDAnaParameters::kPbPb); + + pars->PrintStatus(); + + std::cout<<"creating background object"<Load("libANALYSIS"); + gSystem->Load("libANALYSISalice"); + gSystem->Load("libPWG0base"); + gSystem->Load("libPWG2forward"); + + AliFMDAnaParameters* pars = AliFMDAnaParameters::Instance(); + + if(energy == 900) + pars->SetEnergy(AliFMDAnaParameters::k900); + else if(energy == 7000) + pars->SetEnergy(AliFMDAnaParameters::k7000); + else if(energy == 10000) + pars->SetEnergy(AliFMDAnaParameters::k10000); + else if(energy == 14000) + pars->SetEnergy(AliFMDAnaParameters::k14000); + else if(energy == 5500) + pars->SetEnergy(AliFMDAnaParameters::k5500); + + if(trigger == 0) + pars->SetTriggerDefinition(AliFMDAnaParameters::kMB1); + else if(trigger == 1) + pars->SetTriggerDefinition(AliFMDAnaParameters::kMB2); + + if(mag==0) + pars->SetMagField(AliFMDAnaParameters::k0G); + else if(mag==1) + pars->SetMagField(AliFMDAnaParameters::k5G); + + if(collsystem == 0) + pars->SetCollisionSystem(AliFMDAnaParameters::kPP); + else if(collsystem == 1) + pars->SetCollisionSystem(AliFMDAnaParameters::kPbPb); + + pars->SetRealData(realdata); + + pars->PrintStatus(); + pars->Init(kTRUE,AliFMDAnaParameters::kBackgroundCorrection); + std::cout<<"creating energy distribution object"<Load("libANALYSIS"); + gSystem->Load("libANALYSISalice"); + gSystem->Load("libPWG2forward"); + + gStyle->SetTextFont(132); + gStyle->SetLabelFont(132,"X"); + gStyle->SetLabelFont(132,"Y"); + gStyle->SetLabelFont(132,"Z"); + gStyle->SetTitleFont(132,"X"); + gStyle->SetTitleFont(132,"Y"); + gStyle->SetTitleFont(132,"Z"); + + AliFMDAnaParameters* pars = AliFMDAnaParameters::Instance(); + + if(energy == 900) + pars->SetEnergy(AliFMDAnaParameters::k900); + else if(energy == 7000) + pars->SetEnergy(AliFMDAnaParameters::k7000); + else if(energy == 10000) + pars->SetEnergy(AliFMDAnaParameters::k10000); + else if(energy == 14000) + pars->SetEnergy(AliFMDAnaParameters::k14000); + + if(trigger == 0) + pars->SetTriggerDefinition(AliFMDAnaParameters::kMB1); + else if(trigger == 1) + pars->SetTriggerDefinition(AliFMDAnaParameters::kMB2); + + if(mag==0) + pars->SetMagField(AliFMDAnaParameters::k0G); + else if(mag==1) + pars->SetMagField(AliFMDAnaParameters::k5G); + + if(collsystem == 0) + pars->SetCollisionSystem(AliFMDAnaParameters::kPP); + else if(collsystem == 1) + pars->SetCollisionSystem(AliFMDAnaParameters::kPbPb); + + pars->PrintStatus(); + + std::cout<<"creating sharing efficiency object"<