From f58a4769819ca460dccefd23e8e9e47bf1397d3f Mon Sep 17 00:00:00 2001 From: hdalsgaa Date: Tue, 19 May 2009 21:03:30 +0000 Subject: [PATCH] Analysis with correction for double hits (work in progress) and analysis independent of OCDB and geometry --- FMD/analysis/AliFMDAnaParameters.cxx | 68 +++++++++++++- FMD/analysis/AliFMDAnaParameters.h | 3 +- FMD/analysis/AliFMDAnalysisTaskDensity.cxx | 93 +++++++++++++++---- FMD/analysis/AliFMDAnalysisTaskDensity.h | 5 +- FMD/analysis/AliFMDAnalysisTaskDndeta.cxx | 10 +- .../AliFMDAnalysisTaskGenerateBackground.cxx | 55 ++++++----- .../AliFMDAnalysisTaskGenerateBackground.h | 1 + FMD/analysis/AliFMDAnalysisTaskSharing.cxx | 13 ++- FMD/analysis/FMDanalysisLinkDef.h | 2 + FMD/analysis/RunAliEnFMDAnalysisSE.C | 4 +- FMD/analysis/RunLocalFMDAnalysisSE.C | 5 + FMD/libFMDanalysis.pkg | 3 +- 12 files changed, 210 insertions(+), 52 deletions(-) diff --git a/FMD/analysis/AliFMDAnaParameters.cxx b/FMD/analysis/AliFMDAnaParameters.cxx index 98f423a81c9..89ee4fc03eb 100644 --- a/FMD/analysis/AliFMDAnaParameters.cxx +++ b/FMD/analysis/AliFMDAnaParameters.cxx @@ -24,12 +24,14 @@ #include "AliFMDAnaParameters.h" // ALIFMDPARAMETERS_H #include // ALICDBMANAGER_H #include // ALICDBMANAGER_H +#include "AliFMDRing.h" #include #include #include #include #include #include +#include //==================================================================== ClassImp(AliFMDAnaParameters) @@ -238,7 +240,71 @@ TH2F* AliFMDAnaParameters::GetBackgroundCorrection(Int_t det, return fBackground->GetBgCorrection(det,ring,vtxbin); } - +//_____________________________________________________________________ +Float_t AliFMDAnaParameters::GetPhiFromSector(UShort_t det, Char_t ring, UShort_t sec) +{ + Int_t nsec = (ring == 'I' ? 20 : 40); + Float_t basephi = 0; + if(det == 1) + basephi = 1.72787594; + if(det == 2 && ring == 'I') + basephi = 0.15707963; + if(det == 2 && ring == 'O') + basephi = 0.078539818; + if(det == 3 && ring == 'I') + basephi = 2.984513044; + if(det == 3 && ring == 'O') + basephi = 3.06305289; + + Float_t step = 2*TMath::Pi() / nsec; + Float_t phi = 0; + if(det == 3) + phi = basephi - sec*step; + else + phi = basephi + sec*step; + + if(phi < 0) + phi = phi +2*TMath::Pi(); + if(phi > 2*TMath::Pi() ) + phi = phi - 2*TMath::Pi(); + + return phi; +} +//_____________________________________________________________________ +Float_t AliFMDAnaParameters::GetEtaFromStrip(UShort_t det, Char_t ring, UShort_t sec, UShort_t strip, Float_t zvtx) +{ + AliFMDRing fmdring(ring); + fmdring.Init(); + Float_t rad = fmdring.GetMaxR()-fmdring.GetMinR(); + Float_t segment = rad / fmdring.GetNStrips(); + Float_t r = fmdring.GetMinR() + segment*strip; + Float_t z = 0; + Int_t hybrid = sec / 2; + + if(det == 1) { + if(!(hybrid%2)) z = 320.266; else z = 319.766; + } + if(det == 2 && ring == 'I' ) { + if(!(hybrid%2)) z = 83.666; else z = 83.166; + } + if(det == 2 && ring == 'O' ) { + if(!(hybrid%2)) z = 74.966; else z = 75.466; + } + if(det == 3 && ring == 'I' ) { + if(!(hybrid%2)) z = -63.066; else z = -62.566; + } + if(det == 3 && ring == 'O' ) { + if(!(hybrid%2)) z = -74.966; else z = -75.466; + } + + //std::cout< #include "TAxis.h" #include "TH2F.h" +#include "TF1.h" #include "AliFMDAnalysisTaskDensity.h" #include "AliAnalysisManager.h" #include "AliESDFMD.h" @@ -59,6 +60,13 @@ AliFMDAnalysisTaskDensity::AliFMDAnalysisTaskDensity(const char* name, Bool_t SE DefineInput (1, AliESDVertex::Class()); DefineOutput(0, TList::Class()); } + + fFuncPos = new TF1("funcPos","pol1",0,6); + fFuncPos->SetParameters(0.99925,0.00298301); + fFuncNeg = new TF1("funcNeg","pol1",-6,0); + fFuncNeg->SetParameters(0.987583,-0.0101022); + + } //_____________________________________________________________________ void AliFMDAnalysisTaskDensity::CreateOutputObjects() @@ -118,6 +126,10 @@ void AliFMDAnalysisTaskDensity::ConnectInputData(Option_t */*option*/) fVertex = (AliESDVertex*)GetInputData(1); } } + + + + //_____________________________________________________________________ void AliFMDAnalysisTaskDensity::Exec(Option_t */*option*/) { @@ -140,7 +152,9 @@ void AliFMDAnalysisTaskDensity::Exec(Option_t */*option*/) Double_t vertexBinDouble = (vertex[2] + pars->GetVtxCutZ()) / delta; Int_t vtxbin = (Int_t)vertexBinDouble; + + fVertexString.SetString(Form("%d",vtxbin)); //Reset everything for(UShort_t det=1;det<=3;det++) { @@ -171,27 +185,34 @@ void AliFMDAnalysisTaskDensity::Exec(Option_t */*option*/) for(UShort_t sec =0; sec < nsec; sec++) { for(UShort_t strip = 0; strip < nstr; strip++) { Float_t mult = fESD->Multiplicity(det,ring,sec,strip); - Float_t eta = fESD->Eta(det,ring,sec,strip); + //Float_t eta = fESD->Eta(det,ring,sec,strip); if(mult == 0 || mult == AliESDFMD::kInvalidMult) continue; //Particle number cut goes here... - Double_t x,y,z; - geo->Detector2XYZ(det,ring,sec,strip,x,y,z); - 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[2]); - Float_t etacalc = -1*TMath::Log(TMath::Tan(0.5*theta)); - eta = etacalc; - + //Double_t x,y,z; + //geo->Detector2XYZ(det,ring,sec,strip,x,y,z); + // Float_t phi = TMath::ATan2(y,x); + // if(phi<0) + // phi = phi+2*TMath::Pi(); + + Float_t phi = pars->GetPhiFromSector(det,ring,sec); + Float_t eta = pars->GetEtaFromStrip(det,ring,sec,strip,vertex[2]); + //std::cout<GetMPV(det,ring,eta); Float_t s = pars->GetSigma(det,ring,eta); - AliFMDParameters* recopars = AliFMDParameters::Instance(); - - Float_t mult_cut = m-2*s; //0.15;//0.2;//m-3*s;// 0.2;//0.01;//m-2*s;//0.2; + //AliFMDParameters* recopars = AliFMDParameters::Instance(); - mult_cut = (4*recopars->GetPedestalWidth(det,ring,sec,strip))/(recopars->GetPulseGain(det,ring,sec,strip)*recopars->GetDACPerMIP()); + Float_t mult_cut = 0.15;//m-2*s;//0.15;//0.2;//m-3*s;// 0.2;//0.01;//m-2*s;//0.2; + if(ring == 'I') + mult_cut = 0.10; + //mult_cut = (5*recopars->GetPedestalWidth(det,ring,sec,strip))/(recopars->GetPulseGain(det,ring,sec,strip)*recopars->GetDACPerMIP()); Float_t nParticles = 0; if(fESD->GetUniqueID() == kTRUE) { @@ -232,12 +253,21 @@ void AliFMDAnalysisTaskDensity::Exec(Option_t */*option*/) Float_t correction = GetAcceptanceCorrection(ring,strip); + + //std::cout<<"before "<Eval(eta); + else + correction = correction / fFuncPos->Eval(eta); + + // std::cout< 0) hMult->Fill(eta,phi,nParticles); //if(det == 1 && ring =='I' && nParticles >0) - // std::cout< 0) + // std::cout< 2*TMath::Pi() ) + phi = phi - 2*TMath::Pi(); + + return phi; +} + // //EOF diff --git a/FMD/analysis/AliFMDAnalysisTaskDensity.h b/FMD/analysis/AliFMDAnalysisTaskDensity.h index 33dcf7ffea9..e79a178cf51 100644 --- a/FMD/analysis/AliFMDAnalysisTaskDensity.h +++ b/FMD/analysis/AliFMDAnalysisTaskDensity.h @@ -15,7 +15,7 @@ class AliESDEvent; class TChain; class AliAODEvent; - +class TF1; class AliFMDAnalysisTaskDensity : public AliAnalysisTask @@ -48,6 +48,7 @@ class AliFMDAnalysisTaskDensity : public AliAnalysisTask void SetInputESDFMD(AliESDFMD* esdfmd) {fESD = esdfmd;} void SetInputVertex(AliESDVertex* vertex) {fVertex = vertex;} Bool_t GetEventStatus() { return fStatus; } + Float_t GetPhiFromSector(UShort_t det, Char_t ring, UShort_t sec); private: Int_t fDebug; // Debug flag @@ -58,6 +59,8 @@ class AliFMDAnalysisTaskDensity : public AliAnalysisTask AliESDVertex* fVertex; Bool_t fStandalone; Bool_t fStatus; + TF1* fFuncPos; + TF1* fFuncNeg; ClassDef(AliFMDAnalysisTaskDensity, 0); // Analysis task for FMD analysis }; diff --git a/FMD/analysis/AliFMDAnalysisTaskDndeta.cxx b/FMD/analysis/AliFMDAnalysisTaskDndeta.cxx index 52693340610..b7a1693db54 100644 --- a/FMD/analysis/AliFMDAnalysisTaskDndeta.cxx +++ b/FMD/analysis/AliFMDAnalysisTaskDndeta.cxx @@ -271,20 +271,20 @@ void AliFMDAnalysisTaskDndeta::ProcessPrimary() { Float_t thisStripTrack = fLastTrackByStrip.operator()(det,ring,sec,strip); if(particle->Charge() != 0 && i != thisStripTrack ) { Double_t x,y,z; - AliFMDGeometry* fmdgeo = AliFMDGeometry::Instance(); + /*AliFMDGeometry* fmdgeo = AliFMDGeometry::Instance(); fmdgeo->Detector2XYZ(det,ring,sec,strip,x,y,z); 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)); + Float_t theta = TMath::ATan2(r,z-vertex.At(2));*/ + Float_t eta = pars->GetEtaFromStrip(det,ring,sec,strip,vertex.At(2));//-1*TMath::Log(TMath::Tan(0.5*theta)); TH1F* hHits = (TH1F*)fOutputList->FindObject(Form("hHits_FMD%d%c_vtxbin%d",det,ring,vertexBin)); hHits->Fill(eta); Float_t nstrips = (ring =='O' ? 256 : 512); - // if(det == 1 && ring == 'I') - // std::cout<<"hit in FMD 1I "<0) diff --git a/FMD/analysis/AliFMDAnalysisTaskGenerateBackground.cxx b/FMD/analysis/AliFMDAnalysisTaskGenerateBackground.cxx index 6005da00788..e8d081db274 100644 --- a/FMD/analysis/AliFMDAnalysisTaskGenerateBackground.cxx +++ b/FMD/analysis/AliFMDAnalysisTaskGenerateBackground.cxx @@ -25,32 +25,18 @@ ClassImp(AliFMDAnalysisTaskGenerateBackground) //_____________________________________________________________________ AliFMDAnalysisTaskGenerateBackground::AliFMDAnalysisTaskGenerateBackground(): -AliAnalysisTaskSE(), - fListOfHits(), - fListOfPrimaries(), - fListOfCorrection(), - fVertexBins(), - fLastTrackByStrip(), - fZvtxCut(10), - fNvtxBins(10), - fNbinsEta(100), - fBackground(0) +AliAnalysisTaskSE() { // Default constructor } //_____________________________________________________________________ AliFMDAnalysisTaskGenerateBackground::AliFMDAnalysisTaskGenerateBackground(const char* name): AliAnalysisTaskSE(name), - fListOfHits(), - fListOfPrimaries(), - fListOfCorrection(), - fVertexBins(), - fLastTrackByStrip(), fZvtxCut(10), fNvtxBins(10), - fNbinsEta(100), - fBackground(0) + fNbinsEta(100) { + DefineOutput(1, TList::Class()); DefineOutput(2, TList::Class()); DefineOutput(3, TH1F::Class()); @@ -92,8 +78,14 @@ void AliFMDAnalysisTaskGenerateBackground::UserCreateOutputObjects() } } } + TH1F* doubleHits = new TH1F("DoubleHits", "DoubleHits", + fNbinsEta, -6,6); + TH1F* allHits = new TH1F("allHits", "allHits", + fNbinsEta, -6,6); fVertexBins.SetName("VertexBins"); fVertexBins.GetXaxis()->Set(fNvtxBins,-1*fZvtxCut,fZvtxCut); + fListOfHits.Add(allHits); + fListOfHits.Add(doubleHits); } //_____________________________________________________________________ void AliFMDAnalysisTaskGenerateBackground::Init() @@ -110,8 +102,9 @@ void AliFMDAnalysisTaskGenerateBackground::UserExec(Option_t */*option*/) { fLastTrackByStrip.Reset(-1); + fHitsByStrip.Reset(0); AliMCEvent* mcevent = MCEvent(); - //AliFMDAnaParameters* pars = AliFMDAnaParameters::Instance(); + AliFMDAnaParameters* pars = AliFMDAnaParameters::Instance(); AliMCParticle* particle = 0; AliStack* stack = mcevent->Stack(); @@ -168,6 +161,22 @@ void AliFMDAnalysisTaskGenerateBackground::UserExec(Option_t */*option*/) 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.operator()(det,ring,sec,strip) +=1; + TH1F* allHits = (TH1F*)fListOfHits.FindObject("allHits"); + TH1F* doubleHits = (TH1F*)fListOfHits.FindObject("DoubleHits"); + + if(fHitsByStrip.operator()(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.operator()(det,ring,sec,strip) = (Float_t)i; if(strip >0) @@ -179,9 +188,7 @@ void AliFMDAnalysisTaskGenerateBackground::UserExec(Option_t */*option*/) } - - - + PostData(1, &fListOfHits); PostData(2, &fListOfPrimaries); PostData(3, &fVertexBins); @@ -189,8 +196,14 @@ void AliFMDAnalysisTaskGenerateBackground::UserExec(Option_t */*option*/) //_____________________________________________________________________ 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() { diff --git a/FMD/analysis/AliFMDAnalysisTaskGenerateBackground.h b/FMD/analysis/AliFMDAnalysisTaskGenerateBackground.h index 6b0156054ad..2b2b442915b 100644 --- a/FMD/analysis/AliFMDAnalysisTaskGenerateBackground.h +++ b/FMD/analysis/AliFMDAnalysisTaskGenerateBackground.h @@ -44,6 +44,7 @@ class AliFMDAnalysisTaskGenerateBackground : public AliAnalysisTaskSE TList fListOfCorrection; TH1F fVertexBins; AliFMDFloatMap fLastTrackByStrip; + AliFMDFloatMap fHitsByStrip; Float_t fZvtxCut; Int_t fNvtxBins; Int_t fNbinsEta; diff --git a/FMD/analysis/AliFMDAnalysisTaskSharing.cxx b/FMD/analysis/AliFMDAnalysisTaskSharing.cxx index 0d34cddfbf6..901e64eb435 100644 --- a/FMD/analysis/AliFMDAnalysisTaskSharing.cxx +++ b/FMD/analysis/AliFMDAnalysisTaskSharing.cxx @@ -158,7 +158,8 @@ void AliFMDAnalysisTaskSharing::Exec(Option_t */*option*/) if(mult == AliESDFMD::kInvalidMult || mult == 0) continue; - Double_t eta = EtaFromStrip(det,ring,sec,strip,vertex[2]);//fmd->Eta(det,ring,sec,strip); + //Double_t eta = EtaFromStrip(det,ring,sec,strip,vertex[2]);//fmd->Eta(det,ring,sec,strip); + Double_t eta = fmd->Eta(det,ring,sec,strip); //std::cout<Eta(det,ring,sec,strip)<Fill(mult); @@ -212,12 +213,16 @@ Float_t AliFMDAnalysisTaskSharing::GetMultiplicityOfStrip(Float_t mult, Float_t merged_energy = 0; //Float_t nParticles = 0; Float_t cutLow = 0.15; - AliFMDParameters* recopars = AliFMDParameters::Instance(); - cutLow = (4*recopars->GetPedestalWidth(det,ring,sec,strip))/(recopars->GetPulseGain(det,ring,sec,strip)*recopars->GetDACPerMIP()); + if(ring == 'I') + cutLow = 0.1; + //cutLow = 0; + //AliFMDParameters* recopars = AliFMDParameters::Instance(); + //cutLow = (5*recopars->GetPedestalWidth(det,ring,sec,strip))/(recopars->GetPulseGain(det,ring,sec,strip)*recopars->GetDACPerMIP()); - Float_t cutHigh = pars->GetMPV(det,ring,eta);// - pars->GetSigma(det,ring,eta); + + Float_t cutHigh = pars->GetMPV(det,ring,eta) - 3*pars->GetSigma(det,ring,eta); // Float_t cutPart = pars->GetMPV(det,ring,eta) - 5*pars->GetSigma(det,ring,eta); Float_t Etotal = mult; diff --git a/FMD/analysis/FMDanalysisLinkDef.h b/FMD/analysis/FMDanalysisLinkDef.h index f27998511c9..981718c2c6f 100644 --- a/FMD/analysis/FMDanalysisLinkDef.h +++ b/FMD/analysis/FMDanalysisLinkDef.h @@ -39,6 +39,8 @@ #pragma link C++ class AliFMDAnalysisTaskDensity+; #pragma link C++ class AliFMDAnalysisTaskBackgroundCorrection+; #pragma link C++ class AliFMDAnalysisTaskCollector+; +#pragma link C++ class AliFMDAnalysisTaskGenerateBackground+; + #else # error Not for compilation #endif diff --git a/FMD/analysis/RunAliEnFMDAnalysisSE.C b/FMD/analysis/RunAliEnFMDAnalysisSE.C index 2b2ffa0fd56..8df236437d0 100644 --- a/FMD/analysis/RunAliEnFMDAnalysisSE.C +++ b/FMD/analysis/RunAliEnFMDAnalysisSE.C @@ -47,7 +47,7 @@ void RunAliEnFMDAnalysisSE(const Char_t* collectionName="collection.xml", const AliCDBManager* cdb = AliCDBManager::Instance(); cdb->SetDefaultStorage(cdbPath); - cdb->SetSpecificStorage("FMD/*","local://$ALICE_ROOT"); + cdb->SetRun(0); AliFMDAnaParameters* pars = AliFMDAnaParameters::Instance(); @@ -63,7 +63,7 @@ void RunAliEnFMDAnalysisSE(const Char_t* collectionName="collection.xml", const timer.Start(); if (mgr->InitAnalysis()) { mgr->PrintStatus(); - mgr->StartAnalysis("local",chain, 2000); + mgr->StartAnalysis("local",chain, 5000); } timer.Stop(); timer.Print(); diff --git a/FMD/analysis/RunLocalFMDAnalysisSE.C b/FMD/analysis/RunLocalFMDAnalysisSE.C index 0914202ab3c..c44b3661443 100644 --- a/FMD/analysis/RunLocalFMDAnalysisSE.C +++ b/FMD/analysis/RunLocalFMDAnalysisSE.C @@ -59,8 +59,13 @@ void RunLocalFMDAnalysisSE(const Char_t* filename= "AliESDs.root", const Char_t* AliFMDAnaParameters* pars = AliFMDAnaParameters::Instance(); pars->Init(); + //if (AliGeomManager::GetGeometry() == NULL) + // AliGeomManager::LoadGeometry("geometry.root"); + if (AliGeomManager::GetGeometry() == NULL) AliGeomManager::LoadGeometry(); + AliFMDParameters* recopars = AliFMDParameters::Instance(); + recopars->Init(); AliFMDGeometry* geo = AliFMDGeometry::Instance(); geo->Init(); diff --git a/FMD/libFMDanalysis.pkg b/FMD/libFMDanalysis.pkg index 74308f8c9eb..3bb69ce3346 100644 --- a/FMD/libFMDanalysis.pkg +++ b/FMD/libFMDanalysis.pkg @@ -10,7 +10,8 @@ SRCS = analysis/AliFMDAnaParameters.cxx \ analysis/AliFMDAnalysisTaskBackgroundCorrection.cxx \ analysis/AliFMDAnalysisTaskCollector.cxx \ analysis/AliFMDAnalysisTaskDndeta.cxx \ - analysis/AliFMDAnalysisTaskSE.cxx + analysis/AliFMDAnalysisTaskSE.cxx \ + analysis/AliFMDAnalysisTaskGenerateBackground.cxx HDRS = $(SRCS:.cxx=.h) -- 2.39.3