From 11fd38f5d57060a0e9e2943aaf36373324d3bf36 Mon Sep 17 00:00:00 2001 From: wiechula Date: Mon, 5 Dec 2011 23:32:00 +0000 Subject: [PATCH] o AliTRDTenderSupply: updates from Markus Fasel o AliESDv0KineCuts: updates from Linus Feldkamp - initialise KF event internally o AliAnalysisTaskPIDqa: updates from Michael Weber - provide lists of identified particles from V0 decays --- ANALYSIS/AliAnalysisTaskPIDqa.cxx | 283 +++++++++++++++++- ANALYSIS/AliAnalysisTaskPIDqa.h | 13 + ANALYSIS/AliESDv0KineCuts.cxx | 50 +++- ANALYSIS/AliESDv0KineCuts.h | 10 +- .../TenderSupplies/AliTRDTenderSupply.cxx | 256 +++++++++++++--- ANALYSIS/TenderSupplies/AliTRDTenderSupply.h | 28 +- 6 files changed, 572 insertions(+), 68 deletions(-) diff --git a/ANALYSIS/AliAnalysisTaskPIDqa.cxx b/ANALYSIS/AliAnalysisTaskPIDqa.cxx index 7be798d23bc..9f7f5f77723 100644 --- a/ANALYSIS/AliAnalysisTaskPIDqa.cxx +++ b/ANALYSIS/AliAnalysisTaskPIDqa.cxx @@ -39,6 +39,10 @@ #include #include +#include +#include +#include +#include #include "AliAnalysisTaskPIDqa.h" @@ -49,6 +53,11 @@ ClassImp(AliAnalysisTaskPIDqa) AliAnalysisTaskPIDqa::AliAnalysisTaskPIDqa(): AliAnalysisTaskSE(), fPIDResponse(0x0), +fV0cuts(0x0), +fV0electrons(0x0), +fV0pions(0x0), +fV0kaons(0x0), +fV0protons(0x0), fListQA(0x0), fListQAits(0x0), fListQAitsSA(0x0), @@ -59,7 +68,8 @@ fListQAtof(0x0), fListQAemcal(0x0), fListQAhmpid(0x0), fListQAtofhmpid(0x0), -fListQAtpctof(0x0) +fListQAtpctof(0x0), +fListQAV0(0x0) { // // Dummy constructor @@ -70,6 +80,11 @@ fListQAtpctof(0x0) AliAnalysisTaskPIDqa::AliAnalysisTaskPIDqa(const char* name): AliAnalysisTaskSE(name), fPIDResponse(0x0), +fV0cuts(0x0), +fV0electrons(0x0), +fV0pions(0x0), +fV0kaons(0x0), +fV0protons(0x0), fListQA(0x0), fListQAits(0x0), fListQAitsSA(0x0), @@ -80,7 +95,8 @@ fListQAtof(0x0), fListQAemcal(0x0), fListQAhmpid(0x0), fListQAtofhmpid(0x0), -fListQAtpctof(0x0) +fListQAtpctof(0x0), +fListQAV0(0x0) { // // Default constructor @@ -95,6 +111,13 @@ AliAnalysisTaskPIDqa::~AliAnalysisTaskPIDqa() // // Destructor // + + delete fV0cuts; + delete fV0electrons; + delete fV0pions; + delete fV0kaons; + delete fV0protons; + if (!AliAnalysisManager::GetAnalysisManager()->IsProofMode()) delete fListQA; } @@ -116,6 +139,15 @@ void AliAnalysisTaskPIDqa::UserCreateOutputObjects() fPIDResponse=inputHandler->GetPIDResponse(); if (!fPIDResponse) AliError("PIDResponse object was not created"); + // V0 Kine cuts + fV0cuts = new AliESDv0KineCuts; + + // V0 PID Obj arrays + fV0electrons = new TObjArray; + fV0pions = new TObjArray; + fV0kaons = new TObjArray; + fV0protons = new TObjArray; + // fListQA=new TList; fListQA->SetOwner(); @@ -160,6 +192,10 @@ void AliAnalysisTaskPIDqa::UserCreateOutputObjects() fListQAtofhmpid->SetOwner(); fListQAtofhmpid->SetName("TOF_HMPID"); + fListQAV0=new TList; + fListQAV0->SetOwner(); + fListQAV0->SetName("V0decay"); + fListQA->Add(fListQAits); fListQA->Add(fListQAitsSA); fListQA->Add(fListQAitsPureSA); @@ -170,6 +206,7 @@ void AliAnalysisTaskPIDqa::UserCreateOutputObjects() fListQA->Add(fListQAhmpid); fListQA->Add(fListQAtpctof); fListQA->Add(fListQAtofhmpid); + fListQA->Add(fListQAV0); SetupITSqa(); SetupTPCqa(); @@ -179,6 +216,7 @@ void AliAnalysisTaskPIDqa::UserCreateOutputObjects() SetupHMPIDqa(); SetupTPCTOFqa(); SetupTOFHMPIDqa(); + SetupV0qa(); PostData(1,fListQA); } @@ -194,6 +232,8 @@ void AliAnalysisTaskPIDqa::UserExec(Option_t */*option*/) AliVEvent *event=InputEvent(); if (!event||!fPIDResponse) return; + // Start with the V0 task (only possible for ESDs?) + FillV0PIDlist(); FillITSqa(); FillTPCqa(); @@ -206,10 +246,108 @@ void AliAnalysisTaskPIDqa::UserExec(Option_t */*option*/) FillTPCTOFqa(); FillTOFHMPIDqa(); + // Clear the V0 PID arrays + ClearV0PIDlist(); + + + PostData(1,fListQA); } //______________________________________________________________________________ +void AliAnalysisTaskPIDqa::FillV0PIDlist(){ + + // + // Fill the PID object arrays holding the pointers to identified particle tracks + // + + // Dynamic cast to ESD events (DO NOTHING for AOD events) + AliESDEvent *event = dynamic_cast(InputEvent()); + if ( !event ) return; + + if(TString(event->GetBeamType())=="Pb-Pb" || TString(event->GetBeamType())=="A-A"){ + fV0cuts->SetMode(AliESDv0KineCuts::kPurity,AliESDv0KineCuts::kPbPb); + } + else{ + fV0cuts->SetMode(AliESDv0KineCuts::kPurity,AliESDv0KineCuts::kPP); + } + + // V0 selection + // set event + fV0cuts->SetEvent(event); + + // loop over V0 particles + for(Int_t iv0=0; iv0GetNumberOfV0s();iv0++){ + + AliESDv0 *v0 = (AliESDv0 *) event->GetV0(iv0); + + if(!v0) continue; + if(v0->GetOnFlyStatus()) continue; + + // Get the particle selection + Bool_t foundV0 = kFALSE; + Int_t pdgV0, pdgP, pdgN; + + foundV0 = fV0cuts->ProcessV0(v0, pdgV0, pdgP, pdgN); + if(!foundV0) continue; + + Int_t iTrackP = v0->GetPindex(); // positive track + Int_t iTrackN = v0->GetNindex(); // negative track + + // v0 Armenteros plot (QA) + Float_t armVar[2] = {0.0,0.0}; + fV0cuts->Armenteros(v0, armVar); + + TH2 *h=(TH2*)fListQAV0->At(0); + if (!h) continue; + h->Fill(armVar[0],armVar[1]); + + // fill the Object arrays + // positive particles + if( pdgP == -11){ + fV0electrons->Add((AliVTrack*)event->GetTrack(iTrackP)); + } + else if( pdgP == 211){ + fV0pions->Add((AliVTrack*)event->GetTrack(iTrackP)); + } + else if( pdgP == 321){ + fV0kaons->Add((AliVTrack*)event->GetTrack(iTrackP)); + } + else if( pdgP == 2212){ + fV0protons->Add((AliVTrack*)event->GetTrack(iTrackP)); + } + + // negative particles + if( pdgN == 11){ + fV0electrons->Add((AliVTrack*)event->GetTrack(iTrackN)); + } + else if( pdgN == -211){ + fV0pions->Add((AliVTrack*)event->GetTrack(iTrackN)); + } + else if( pdgN == -321){ + fV0kaons->Add((AliVTrack*)event->GetTrack(iTrackN)); + } + else if( pdgN == -2212){ + fV0protons->Add((AliVTrack*)event->GetTrack(iTrackN)); + } + + + } +} +//______________________________________________________________________________ +void AliAnalysisTaskPIDqa::ClearV0PIDlist(){ + + // + // Clear the PID object arrays + // + + fV0electrons->Clear(); + fV0pions->Clear(); + fV0kaons->Clear(); + fV0protons->Clear(); + +} +//______________________________________________________________________________ void AliAnalysisTaskPIDqa::FillITSqa() { // @@ -461,9 +599,6 @@ void AliAnalysisTaskPIDqa::FillEMCALqa() // ULong_t status=track->GetStatus(); // not that nice. status bits not in virtual interface - // TPC refit + ITS refit + - // TOF out + TOFpid + - // kTIME if (!( (status & AliVTrack::kEMCALmatch) == AliVTrack::kEMCALmatch) ) continue; Double_t pt=track->Pt(); @@ -474,8 +609,22 @@ void AliAnalysisTaskPIDqa::FillEMCALqa() Double_t nSigma=fPIDResponse->NumberOfSigmasEMCAL(track, (AliPID::EParticleType)0); h->Fill(pt,nSigma); - //EMCAL signal (E/p vs. pT) - h=(TH2*)fListQAemcal->At(1); + } + + //EMCAL signal (E/p vs. pT) for electrons from V0 + for(Int_t itrack = 0; itrack < fV0electrons->GetEntries(); itrack++){ + AliVTrack *track=(AliVTrack*)fV0electrons->At(itrack); + + // + //basic track cuts + // + ULong_t status=track->GetStatus(); + // not that nice. status bits not in virtual interface + if (!( (status & AliVTrack::kEMCALmatch) == AliVTrack::kEMCALmatch) ) continue; + + Double_t pt=track->Pt(); + + TH2 *h=(TH2*)fListQAemcal->At(1); if (h) { Int_t nMatchClus = track->GetEMCALcluster(); @@ -501,6 +650,90 @@ void AliAnalysisTaskPIDqa::FillEMCALqa() } } } + + //EMCAL signal (E/p vs. pT) for pions from V0 + for(Int_t itrack = 0; itrack < fV0pions->GetEntries(); itrack++){ + AliVTrack *track=(AliVTrack*)fV0pions->At(itrack); + + // + //basic track cuts + // + ULong_t status=track->GetStatus(); + // not that nice. status bits not in virtual interface + if (!( (status & AliVTrack::kEMCALmatch) == AliVTrack::kEMCALmatch) ) continue; + + Double_t pt=track->Pt(); + + TH2 *h=(TH2*)fListQAemcal->At(2); + if (h) { + + Int_t nMatchClus = track->GetEMCALcluster(); + Double_t mom = track->P(); + Double_t eop = -1.; + + if(nMatchClus > -1){ + + AliVCluster *matchedClus = (AliVCluster*)event->GetCaloCluster(nMatchClus); + + if(matchedClus){ + + // matched cluster is EMCAL + if(matchedClus->IsEMCAL()){ + + Double_t fClsE = matchedClus->E(); + eop = fClsE/mom; + + h->Fill(pt,eop); + + } + } + } + } + } + + //EMCAL signal (E/p vs. pT) for protons from V0 + for(Int_t itrack = 0; itrack < fV0protons->GetEntries(); itrack++){ + AliVTrack *track=(AliVTrack*)fV0protons->At(itrack); + + // + //basic track cuts + // + ULong_t status=track->GetStatus(); + // not that nice. status bits not in virtual interface + if (!( (status & AliVTrack::kEMCALmatch) == AliVTrack::kEMCALmatch) ) continue; + + Double_t pt=track->Pt(); + + TH2 *hP=(TH2*)fListQAemcal->At(3); + TH2 *hAP=(TH2*)fListQAemcal->At(4); + if (hP && hAP) { + + Int_t nMatchClus = track->GetEMCALcluster(); + Double_t mom = track->P(); + Int_t charge = track->Charge(); + Double_t eop = -1.; + + if(nMatchClus > -1){ + + AliVCluster *matchedClus = (AliVCluster*)event->GetCaloCluster(nMatchClus); + + if(matchedClus){ + + // matched cluster is EMCAL + if(matchedClus->IsEMCAL()){ + + Double_t fClsE = matchedClus->E(); + eop = fClsE/mom; + + if(charge > 0) hP->Fill(pt,eop); + else if(charge < 0) hAP->Fill(pt,eop); + + } + } + } + } + } + } @@ -856,11 +1089,30 @@ void AliAnalysisTaskPIDqa::SetupEMCALqa() 200,-10,10); fListQAemcal->Add(hNsigmaPt); - TH2F *hSigPt = new TH2F("hSigPt_EMCAL", - "EMCAL signal (E/p) vs. p_{T};p_{T} [GeV]; EMCAL signal (E/p) [arb. units]", + + TH2F *hSigPtEle = new TH2F("hSigPt_EMCAL_Ele", + "EMCAL signal (E/p) vs. p_{T} for electrons;p_{T} [GeV]; EMCAL signal (E/p) [arb. units]", + vX->GetNrows()-1,vX->GetMatrixArray(), + 200,0,2); + fListQAemcal->Add(hSigPtEle); + + TH2F *hSigPtPions = new TH2F("hSigPt_EMCAL_Pions", + "EMCAL signal (E/p) vs. p_{T} for pions;p_{T} [GeV]; EMCAL signal (E/p) [arb. units]", + vX->GetNrows()-1,vX->GetMatrixArray(), + 200,0,2); + fListQAemcal->Add(hSigPtPions); + + TH2F *hSigPtProtons = new TH2F("hSigPt_EMCAL_Protons", + "EMCAL signal (E/p) vs. p_{T} for protons;p_{T} [GeV]; EMCAL signal (E/p) [arb. units]", + vX->GetNrows()-1,vX->GetMatrixArray(), + 200,0,2); + fListQAemcal->Add(hSigPtProtons); + + TH2F *hSigPtAntiProtons = new TH2F("hSigPt_EMCAL_Antiprotons", + "EMCAL signal (E/p) vs. p_{T} for antiprotons;p_{T} [GeV]; EMCAL signal (E/p) [arb. units]", vX->GetNrows()-1,vX->GetMatrixArray(), 200,0,2); - fListQAemcal->Add(hSigPt); + fListQAemcal->Add(hSigPtAntiProtons); delete vX; } @@ -934,6 +1186,17 @@ void AliAnalysisTaskPIDqa::SetupTPCTOFqa() delete vX; } +//______________________________________________________________________________ +void AliAnalysisTaskPIDqa::SetupV0qa() +{ + // + // Create the qa objects for V0 Kine cuts + // + + TH2F *hArmenteros = new TH2F("hArmenteros", "Armenteros plot",200,-1.,1.,200,0.,0.4); + fListQAV0->Add(hArmenteros); + +} //______________________________________________________________________________ TVectorD* AliAnalysisTaskPIDqa::MakeLogBinning(Int_t nbinsX, Double_t xmin, Double_t xmax) diff --git a/ANALYSIS/AliAnalysisTaskPIDqa.h b/ANALYSIS/AliAnalysisTaskPIDqa.h index 92f819adaba..821fbbd6230 100644 --- a/ANALYSIS/AliAnalysisTaskPIDqa.h +++ b/ANALYSIS/AliAnalysisTaskPIDqa.h @@ -23,6 +23,7 @@ class AliPIDResponse; class TList; class AliVEvent; +class AliESDv0KineCuts; class AliAnalysisTaskPIDqa : public AliAnalysisTaskSE { @@ -39,6 +40,13 @@ public: private: AliPIDResponse *fPIDResponse; //! PID response Handler + AliESDv0KineCuts *fV0cuts; //! ESD V0 cuts + + TObjArray *fV0electrons; //! array with pointer to identified particles from V0 decays (electrons) + TObjArray *fV0pions; //! array with pointer to identified particles from V0 decays (pions) + TObjArray *fV0kaons; //! array with pointer to identified particles from V0 decays (kaons) + TObjArray *fV0protons; //! array with pointer to identified particles from V0 decays (ptotons) + TList *fListQA; //! list with all QA histograms TList *fListQAits; //! List with ITS QA histograms TList *fListQAitsSA; //! List with ITS SA QA histograms @@ -50,6 +58,7 @@ private: TList *fListQAhmpid; //! List with EMCAL QA histograms TList *fListQAtofhmpid; //! List with EMCAL QA histograms TList *fListQAtpctof; //! List with combined PID from TPC + TOF + TList *fListQAV0; //! List with V0 kine cuts QA histograms void ExecNewRun(); @@ -63,8 +72,12 @@ private: void SetupHMPIDqa(); void SetupTOFHMPIDqa(); void SetupTPCTOFqa(); + void SetupV0qa(); // + void FillV0PIDlist(); + void ClearV0PIDlist(); + // void FillITSqa(); void FillTPCqa(); void FillTRDqa(); diff --git a/ANALYSIS/AliESDv0KineCuts.cxx b/ANALYSIS/AliESDv0KineCuts.cxx index e9181466a4c..421a31eb960 100644 --- a/ANALYSIS/AliESDv0KineCuts.cxx +++ b/ANALYSIS/AliESDv0KineCuts.cxx @@ -54,6 +54,8 @@ AliESDv0KineCuts::AliESDv0KineCuts() : , fGcutInvMass(0.05) , fK0cutChi2NDF(10) , fLcutChi2NDF(10) + , fUseExternalVertex(kFALSE) + , fDeleteVertex(kFALSE) { // // Default constructor @@ -124,6 +126,8 @@ AliESDv0KineCuts::AliESDv0KineCuts(const AliESDv0KineCuts &ref): , fGcutInvMass(0.05) , fK0cutChi2NDF(10) , fLcutChi2NDF(10) + , fUseExternalVertex(kFALSE) + , fDeleteVertex(kFALSE) { // // Copy operator @@ -156,7 +160,8 @@ void AliESDv0KineCuts::Copy(TObject &ref) const { target.fTPCchi2perCls = fTPCchi2perCls; target.fTPCclsRatio = fTPCclsRatio; target.fNoKinks = fNoKinks; - + target.fUseExternalVertex = fUseExternalVertex; //added december 2nd 2011 + target.fDeleteVertex = fDeleteVertex; //added december 2nd 2011 // default gamma cuts values target.fGcutChi2NDF = fGcutChi2NDF; @@ -949,6 +954,17 @@ void AliESDv0KineCuts::SetEvent(AliESDEvent* const event){ AliErrorClass("Invalid input event pointer"); return; } +if (fUseExternalVertex) return; +else{ + if(fPrimaryVertex && fDeleteVertex){ + delete fPrimaryVertex; + fPrimaryVertex=0x0; + } + fPrimaryVertex = new AliKFVertex(*(event->GetPrimaryVertex())); + fDeleteVertex=kTRUE; + } + + } //____________________________________________________________________ @@ -963,13 +979,43 @@ void AliESDv0KineCuts::SetEvent(AliVEvent* const event){ return; } +if (fUseExternalVertex) return; +else{ + if(fPrimaryVertex && fDeleteVertex){ + delete fPrimaryVertex; + fPrimaryVertex=0x0; + } + fPrimaryVertex = new AliKFVertex(*(event->GetPrimaryVertex())); + fDeleteVertex=kTRUE; } + +} + + +//________________________________________________________________ +void AliESDv0KineCuts::UseExternalVertex(Bool_t use_external){ + // + // Reenable primary Vertex from ESD event + // + if (use_external) fUseExternalVertex =kTRUE; + else fUseExternalVertex =kFALSE; +} + + + + //________________________________________________________________ void AliESDv0KineCuts::SetPrimaryVertex(AliKFVertex* const v){ // // set the primary vertex of the event // - fPrimaryVertex = v; + if(fPrimaryVertex && fDeleteVertex){ + delete fPrimaryVertex; + fPrimaryVertex =0x0; + fDeleteVertex = kFALSE; + } + fUseExternalVertex=kTRUE; + fPrimaryVertex = v; // set primary Vertex if(!fPrimaryVertex){ AliErrorClass("Failed to initialize the primary vertex"); return; diff --git a/ANALYSIS/AliESDv0KineCuts.h b/ANALYSIS/AliESDv0KineCuts.h index 9575d63c9fa..9a3ebbf1399 100644 --- a/ANALYSIS/AliESDv0KineCuts.h +++ b/ANALYSIS/AliESDv0KineCuts.h @@ -52,7 +52,9 @@ class AliESDv0KineCuts : public TObject{ // user can select an operation modes [see .cxx for details] void SetMode(Int_t mode, Int_t type); void SetMode(Int_t mode, const char* type); - + void UseExternalVertex(Bool_t use_external=kTRUE); + AliKFParticle *CreateMotherParticle(const AliVTrack* const pdaughter, const AliVTrack* const ndaughter, Int_t pspec, Int_t nspec) const; + void SetCuts(); // setup cuts for selected fMode and fType, see source file for details // // setter functions for V0 cut values // for default values see the constructor @@ -144,9 +146,6 @@ class AliESDv0KineCuts : public TObject{ void Copy(TObject &ref) const; private: - - AliKFParticle *CreateMotherParticle(const AliVTrack* const pdaughter, const AliVTrack* const ndaughter, Int_t pspec, Int_t nspec) const; - void SetCuts(); // setup cuts for selected fMode and fType, see source file for details Bool_t GammaEffCuts(AliESDv0 * const v0); // set of cuts optimized for high gamma efficiency private: @@ -182,7 +181,8 @@ class AliESDv0KineCuts : public TObject{ Float_t fLcutDCA[2]; // DCA between the daughter tracks [min, max] Float_t fLcutVertexR[2]; // radius of the decay point [min, max] Float_t fLcutInvMass[2]; // invariant mass window - + Bool_t fUseExternalVertex; // Is kTRUE if Vertex is set via SetPrimaryVertex() + Bool_t fDeleteVertex; // Is kTRUE if Vertex has been created in SetEvent() function ClassDef(AliESDv0KineCuts, 0); diff --git a/ANALYSIS/TenderSupplies/AliTRDTenderSupply.cxx b/ANALYSIS/TenderSupplies/AliTRDTenderSupply.cxx index 002d72157b3..a6e8f93cdfd 100644 --- a/ANALYSIS/TenderSupplies/AliTRDTenderSupply.cxx +++ b/ANALYSIS/TenderSupplies/AliTRDTenderSupply.cxx @@ -21,19 +21,24 @@ /////////////////////////////////////////////////////////////////////////////// #include +#include +#include #include #include #include #include +#include #include #include #include +#include #include #include #include #include +#include #include #include #include @@ -41,8 +46,10 @@ #include #include #include +#include #include #include +#include "AliTRDCalChamberStatus.h" #include #include "AliTRDTenderSupply.h" @@ -58,19 +65,25 @@ AliTRDTenderSupply::AliTRDTenderSupply() : fChamberGainNew(NULL), fChamberVdriftOld(NULL), fChamberVdriftNew(NULL), - fRefFilename(""), + fRunByRunCorrection(NULL), fPIDmethod(kNNpid), fNormalizationFactor(1.), fPthreshold(0.8), fNBadChambers(0), + fGeoFile(NULL), fGainCorrection(kTRUE), + fLoadReferences(kFALSE), fLoadReferencesFromCDB(kFALSE), - fHasReferences(kFALSE) + fLoadDeadChambers(kFALSE), + fHasReferences(kFALSE), + fHasNewCalibration(kTRUE), + fDebugMode(kFALSE), + fNameRunByRunCorrection() { // // default ctor // - memset(fBadChamberID, 0, sizeof(Int_t) * kNChambers); + memset(fSlicesForPID, 0, sizeof(UInt_t) * 2); } //_____________________________________________________ @@ -82,18 +95,25 @@ AliTRDTenderSupply::AliTRDTenderSupply(const char *name, const AliTender *tender fChamberGainNew(NULL), fChamberVdriftOld(NULL), fChamberVdriftNew(NULL), - fRefFilename(""), + fRunByRunCorrection(NULL), fPIDmethod(kNNpid), fNormalizationFactor(1.), fPthreshold(0.8), fNBadChambers(0), + fGeoFile(NULL), fGainCorrection(kTRUE), + fLoadReferences(kFALSE), fLoadReferencesFromCDB(kFALSE), - fHasReferences(kFALSE) + fLoadDeadChambers(kFALSE), + fHasReferences(kFALSE), + fHasNewCalibration(kTRUE), + fDebugMode(kFALSE), + fNameRunByRunCorrection() { // // named ctor // + memset(fSlicesForPID, 0, sizeof(UInt_t) * 2); memset(fBadChamberID, 0, sizeof(Int_t) * kNChambers); } @@ -121,8 +141,11 @@ void AliTRDTenderSupply::Init() fTender->GetESDhandler()->SetESDpid(fESDpid); } // Load References - if(!fLoadReferencesFromCDB) LoadReferences(); - fESDpid->GetTRDResponse().SetGainNormalisationFactor(fNormalizationFactor); + if(fLoadReferences && !fLoadReferencesFromCDB) LoadReferences(); + //fESDpid->GetTRDResponse().SetGainNormalisationFactor(fNormalizationFactor); + fESDpid->SetTRDslicesForPID(fSlicesForPID[0], fSlicesForPID[1]); + + if(fNameRunByRunCorrection.Length()) LoadRunByRunCorrection(fNameRunByRunCorrection.Data()); // Set Normalisation Factors if(mgr->GetMCtruthEventHandler()){ @@ -147,9 +170,19 @@ void AliTRDTenderSupply::ProcessEvent() if (fTender->RunChanged()){ AliDebug(0, Form("AliTPCTenderSupply::ProcessEvent - Run Changed (%d)\n",fTender->GetRun())); if (fGainCorrection) SetChamberGain(); - if(!fHasReferences) LoadReferences(); + if(fLoadReferences && !fHasReferences) LoadReferences(); + if(fLoadDeadChambers) LoadDeadChambersFromCDB(); + // Load Geometry + if(AliGeomManager::GetGeometry()){ + AliInfo("Geometry already loaded by other tenders"); + } else { + if(fGeoFile) AliInfo(Form("Load geometry from file %s\n", fGeoFile)); + else AliInfo("Load Geometry from OCDB\n"); + AliGeomManager::LoadGeometry(fGeoFile); + } } + fESD = fTender->GetEvent(); if (!fESD) return; Int_t ntracks=fESD->GetNumberOfTracks(); @@ -157,11 +190,34 @@ void AliTRDTenderSupply::ProcessEvent() // // recalculate PID probabilities // + Int_t detectors[kNPlanes]; for(Int_t itrack = 0; itrack < ntracks; itrack++){ + for(Int_t idet = 0; idet < 5; idet++) detectors[idet] = -1; AliESDtrack *track=fESD->GetTrack(itrack); // Recalculate likelihoods if(!(track->GetStatus() & AliESDtrack::kTRDout)) continue; - if(fGainCorrection) ApplyGainCorrection(track); + AliDebug(2, Form("TRD track found, gain correction: %s, Number of bad chambers: %d\n", fGainCorrection ? "Yes" : "No", fNBadChambers)); + if(GetTRDchamberID(track, detectors)){ + if(fGainCorrection && fHasNewCalibration) ApplyGainCorrection(track, detectors); + if(fNBadChambers) MaskChambers(track, detectors); + } + if(fRunByRunCorrection) ApplyRunByRunCorrection(track); + if(fNormalizationFactor != 1.){ + //printf("Gain Factor: %f\n", fNormalizationFactor); + // Renormalize charge + Double_t qslice = -1; + for(Int_t ily = 0; ily < 6; ily++){ + for(Int_t is = 0; is < track->GetNumberOfTRDslices(); is++){ + qslice = track->GetTRDslice(ily, is); + //printf("Doing layer %d slice %d, value %f\n", ily, is, qslice); + if(qslice >0){ + qslice *= fNormalizationFactor; + //printf("qslice new: %f\n", qslice); + track->SetTRDslice(qslice, ily, is); + } + } + } + } switch(fPIDmethod){ case kNNpid: break; @@ -174,6 +230,32 @@ void AliTRDTenderSupply::ProcessEvent() } } +//_____________________________________________________ +void AliTRDTenderSupply::LoadDeadChambersFromCDB(){ + // + // Load Dead Chambers from the OCDB + // + AliDebug(1, "Loading Dead Chambers from the OCDB"); + AliCDBEntry *en = fTender->GetCDBManager()->Get("TRD/Calib/ChamberStatus",fTender->GetRun()); + if(!en){ + AliError("Dead Chambers not in OCDB"); + return; + } + en->GetId().Print(); + + AliTRDCalChamberStatus* chamberStatus = 0; + if(en){ + chamberStatus = (AliTRDCalChamberStatus*)en->GetObject(); + if(!chamberStatus) AliError("List with the dead chambers not found"); + for(Int_t ichamber = 0; ichamber < 540; ichamber++) { + if(!chamberStatus->IsGood(ichamber)){ + //printf("Chamber not installed %d\n",ichamber); + AddBadChamber(ichamber); + } + } + } +} + //_____________________________________________________ void AliTRDTenderSupply::LoadReferences(){ // @@ -203,19 +285,15 @@ void AliTRDTenderSupply::LoadReferences(){ if(ref){ fESDpid->GetTRDResponse().Load(ref); fHasReferences = kTRUE; - AliInfo("Reference distributions loaded into the PID Response"); + AliDebug(1, "Reference distributions loaded into the PID Response"); } else { AliError("References not found"); } } else { // Backward compatibility mode AliInfo("Loading Reference Distributions from ROOT file"); - if(fRefFilename.Length() != 0){ - fESDpid->GetTRDResponse().Load(fRefFilename.Data()); - fHasReferences = kTRUE; - } else{ - AliError("No file defined"); - } + fESDpid->GetTRDResponse().Load("$TRAIN_ROOT/util/tender/LQ1dRef_v3.root"); + fHasReferences = kTRUE; } } @@ -228,22 +306,26 @@ void AliTRDTenderSupply::SetChamberGain(){ //find previous entry from the UserInfo TTree *tree=((TChain*)fTender->GetInputData(0))->GetTree(); if (!tree) { - AliError("Tree not found in ESDhandler"); + fHasNewCalibration = kFALSE; + AliError("Tree not found in ESDhandler"); return; } TList *userInfo=(TList*)tree->GetUserInfo(); if (!userInfo) { + fHasNewCalibration = kFALSE; AliError("No UserInfo found in tree"); return; } TList *cdbList=(TList*)userInfo->FindObject("cdbList"); if (!cdbList) { + fHasNewCalibration = kFALSE; AliError("No cdbList found in UserInfo"); if (AliLog::GetGlobalLogLevel()>=AliLog::kError) userInfo->Print(); return; } + fHasNewCalibration = kTRUE; TIter nextCDB(cdbList); TObjString *os=0x0; @@ -294,11 +376,54 @@ void AliTRDTenderSupply::SetChamberGain(){ } else AliError("No new drift velocity calibration entry found"); - if(!fChamberGainNew || !fChamberVdriftNew) AliError("No recent calibration found"); + if(!fChamberGainNew || !fChamberVdriftNew){ + AliError("No recent calibration found"); + fHasNewCalibration = kFALSE; + } +} + +//_____________________________________________________ +void AliTRDTenderSupply::LoadRunByRunCorrection(const char *filename){ + // + // Define run by run gain correction for the charge + // + + TDirectory *bkp = gDirectory; + TFile *in = TFile::Open(filename); + bkp->cd(); + fRunByRunCorrection = dynamic_cast(in->Get("TRDchargeCorrection")); + delete in; + if(fRunByRunCorrection ) + AliDebug(2, Form("OADB Container has %d runs\n", fRunByRunCorrection->GetNumberOfEntries())); + /* Temporarily out due to a bug in AliOADBContainer + fRunByRunCorrection = new AliOADBContainer("TRDchargeCorrection"); + Int_t status = fRunByRunCorrection->InitFromFile(filename, "TRDchargeCorrection"); + if(!status) AliDebug(1, Form("Run-dependend gain correction factors loaded from OADB file %s", filename)); + else{ + AliDebug(1, "Failed Loading Run-dependend gain correction factors"); + delete fRunByRunCorrection; + fRunByRunCorrection = NULL; + } + */ +} + +//_____________________________________________________ +Bool_t AliTRDTenderSupply::IsBadChamber(Int_t chamberID){ + // + // Check if the chamber id is in the list of bad chambers + // + Bool_t isBad = kFALSE; + for(UInt_t icam = 0; icam < fNBadChambers; icam++) + if(fBadChamberID[icam] == chamberID){ + isBad = kTRUE; + //printf("cross checking: %i \n",chamberID); + break; + } + return isBad; } //_____________________________________________________ -void AliTRDTenderSupply::ApplyGainCorrection(AliESDtrack * track){ +void AliTRDTenderSupply::ApplyGainCorrection(AliESDtrack * track, const Int_t * const chamberID){ // // Apply new gain factors to the track // @@ -308,32 +433,12 @@ void AliTRDTenderSupply::ApplyGainCorrection(AliESDtrack * track){ } if(!(track->GetStatus() & AliESDtrack::kTRDout)) return; - Double_t p = track->GetOuterParam() ? track->GetOuterParam()->P() : track->P(); - if(p < fPthreshold) return; // Apply low momentum cutoff - Bool_t applyCorrectionVdrift = kFALSE; if(fChamberVdriftOld && fChamberVdriftNew) applyCorrectionVdrift = kTRUE; - Int_t chamberID[kNPlanes]; - for(Int_t il = 0; il < kNPlanes; il++) chamberID[il] = -1; - if(!GetTRDchamberID(track, chamberID)) return; - Int_t nTrackletsPID = 0, nTracklets = track->GetTRDntracklets(); for(Int_t iplane = 0; iplane < kNPlanes; iplane++){ if(chamberID[iplane] < 0) continue; - - // Mask out bad chambers - Bool_t isMasked = kFALSE; - for(UInt_t icam = 0; icam < fNBadChambers; icam++) - if(fBadChamberID[icam] == chamberID[iplane]){ - isMasked = kTRUE; - break; - } - if(isMasked){ - for(Int_t islice = 0; islice < track->GetNumberOfTRDslices(); islice++){ - track->SetTRDslice(0, iplane, islice); - } - continue; - } + if(IsBadChamber(chamberID[iplane])) continue; // Don't apply gain correction for chambers which are in the list of bad chambers // Take old and new gain factor and make ratio Double_t facOld = fChamberGainOld->GetValue(chamberID[iplane]); @@ -346,15 +451,59 @@ void AliTRDTenderSupply::ApplyGainCorrection(AliESDtrack * track){ correction *= vDriftNew/vDriftOld; } AliDebug(2, Form("Applying correction factor %f\n", correction)); - Int_t nSlices = 0; for(Int_t islice = 0; islice < track->GetNumberOfTRDslices(); islice++){ Double_t qslice = track->GetTRDslice(iplane, islice); if(qslice <= 0.) continue; track->SetTRDslice(qslice / correction, iplane, islice); - nSlices++; } - if(nSlices) nTrackletsPID++; } +} + +//_____________________________________________________ +void AliTRDTenderSupply::ApplyRunByRunCorrection(AliESDtrack *const track) { + // + // Equalize charge distribution by applying run-by-run correction (multiplicative) + // + + TVectorD *corrfactor = dynamic_cast(fRunByRunCorrection->GetObject(fTender->GetRun())); + if(!corrfactor) AliDebug(2, "Couldn't derive gain correction factor from OADB"); + else AliDebug(2, Form("Gain factor from OADB %f", (*corrfactor)[0])); + Double_t slice = 0; + for(Int_t ily = 0; ily < kNPlanes; ily++){ + for(Int_t islice = 0; islice < track->GetNumberOfTRDslices(); islice++){ + slice = track->GetTRDslice(ily, islice); + if(slice < 0.001) continue; // do not modify slices which are 0 or negative + slice *= (*corrfactor)[0]; + track->SetTRDslice(slice, ily, islice); + } + } +} + +//_____________________________________________________ +void AliTRDTenderSupply::MaskChambers(AliESDtrack *const track, const Int_t * const chamberID){ + // + // Mask out chambers which are in the list of bad chambers + // Set chamber signal to 0 and reduce the number of tracklets used for PID + // + AliDebug(2, "Masking bad chambers for TRD track"); + Int_t nTrackletsPID = 0, nslice = 0, nTracklets = track->GetTRDntracklets(); + Bool_t badChamber = kFALSE; + //Int_t nbad = 0 ; // Number of bad chambers which contain also a signal + //Int_t nsliceBad = 0; // Number of slices in tracklet in a bad chamber + for(Int_t iplane = 0; iplane < kNPlanes; iplane++){ + badChamber = kFALSE; + nslice = 0; //nsliceBad = 0; + if(IsBadChamber(chamberID[iplane])) badChamber = kTRUE; + for(Int_t islice = 0; islice < track->GetNumberOfTRDslices(); islice++){ + if(badChamber){ + //if(track->GetTRDslice(iplane, islice)) nsliceBad++; + track->SetTRDslice(-1, iplane, islice); + } else if(track->GetTRDslice(iplane, islice) > 0.001) nslice++; + } + //if(nsliceBad) nbad++; + if(nslice > 0) nTrackletsPID++; + } + //if(nbad) track->SetTRDncls(track->GetTRDncls() - 20 * nbad); // subtract mean number of clusters per tracklet for bad tracklets // Use nTrackletsPID to indicate the number of tracklets from good // chambers so they are used for the PID track->SetTRDntracklets(nTrackletsPID | (nTracklets << 3)); @@ -365,14 +514,17 @@ Bool_t AliTRDTenderSupply::GetTRDchamberID(AliESDtrack * const track, Int_t *det // // Calculate TRD chamber ID // + Double_t p = track->GetOuterParam() ? track->GetOuterParam()->P() : track->P(); + if(p < fPthreshold) return kFALSE; // Apply low momentum cutoff + Double_t xLayer[kNPlanes] = {300.2, 312.8, 325.4, 338., 350.6, 363.2}; Double_t etamin[kNStacks] = {0.536, 0.157, -0.145, -0.527,-0.851}; Double_t etamax[kNStacks] = {0.851, 0.527, 0.145, -0.157,-0.536}; for(Int_t ily = 0; ily < kNPlanes; ily++) detectors[ily] = -1; const AliExternalTrackParam *trueparam = NULL; - if(track->GetTPCInnerParam()) trueparam = track->GetTPCInnerParam(); - else if(track->GetOuterParam()) trueparam = track->GetOuterParam(); + if(track->GetOuterParam()) trueparam = track->GetOuterParam(); + else if(track->GetTPCInnerParam()) trueparam = track->GetTPCInnerParam(); else if(track->GetInnerParam()) trueparam = track->GetInnerParam(); if(!trueparam){ AliDebug(2, "No Track Params"); @@ -383,12 +535,24 @@ Bool_t AliTRDTenderSupply::GetTRDchamberID(AliESDtrack * const track, Int_t *det Double_t pos[3]; Int_t nDet = 0; for(Int_t ily = 0; ily < kNPlanes; ily++){ - if(!workparam.PropagateTo(xLayer[ily], fESD->GetMagneticField())) { + if(!AliTrackerBase::PropagateTrackToBxByBz(&workparam, xLayer[ily], 0.139, 100)){ // Assuming the pion mass AliDebug(2, "Propagation failed"); break; } workparam.GetXYZ(pos); Double_t trackAlpha = TMath::ATan2(pos[1], pos[0]); + if(fDebugMode){ + // Compare to simple propagation without magnetic field + AliExternalTrackParam workparam1(*trueparam); // Do calculation on working Copy + Double_t pos1[3]; + if(!workparam1.PropagateTo(xLayer[ily], fESD->GetMagneticField())) { + AliDebug(2, "Propagation failed"); + break; + } + workparam1.GetXYZ(pos1); + Double_t trackAlpha1 = TMath::ATan2(pos1[1], pos1[0]); + AliDebug(2, Form("Alpha: Old %f, New %f, diff %f", trackAlpha1, trackAlpha, trackAlpha-trackAlpha1)); + } if(trackAlpha < 0) trackAlpha = 2 * TMath::Pi() + trackAlpha; Double_t secAlpha = 2 * TMath::Pi() / 18.; diff --git a/ANALYSIS/TenderSupplies/AliTRDTenderSupply.h b/ANALYSIS/TenderSupplies/AliTRDTenderSupply.h index dc380ab2611..82ca5adb88a 100644 --- a/ANALYSIS/TenderSupplies/AliTRDTenderSupply.h +++ b/ANALYSIS/TenderSupplies/AliTRDTenderSupply.h @@ -17,6 +17,7 @@ class AliTRDpidRecalculator; class AliTRDCalDet; class AliESDEvent; +class AliOADBContainer; class AliTRDTenderSupply: public AliTenderSupply { @@ -30,17 +31,22 @@ public: AliTRDTenderSupply(const char *name, const AliTender *tender=NULL); virtual ~AliTRDTenderSupply(); - void SetLoadReferencesFromCDB() { fLoadReferencesFromCDB = kTRUE; } - void SetRefFile(const char* file) {fRefFilename=file;} + void SetRunByRunCorrection(const char *filename) { fNameRunByRunCorrection = filename; } + void SetLoadReferencesFromCDB() { fLoadReferences = kTRUE; fLoadReferencesFromCDB = kTRUE; } + void SetLoadReferencesFromFile() { fLoadReferences = kTRUE; fLoadReferencesFromCDB = kFALSE; } + void SetLoadDeadChambersFromCDB(){ fLoadDeadChambers = kTRUE;} ; void SetPIDmethod(Int_t pidMethod) { fPIDmethod = pidMethod; } void SetNormalizationFactor(Double_t norm) { fNormalizationFactor = norm; } void SetCalibLowpThreshold(Double_t pmin) { fPthreshold = pmin; }; - + void SetGeoFile(const char *filename) { fGeoFile = filename; } + void SetDebugMode() { fDebugMode = kTRUE; } + virtual void Init(); virtual void ProcessEvent(); void SwitchOnGainCorrection() { fGainCorrection = kTRUE; } void SwitchOffGainCorrection() { fGainCorrection = kFALSE; } + void SetSlicesForPID(UInt_t min, UInt_t max) { fSlicesForPID[0] = min; fSlicesForPID[1] = max;} void AddBadChamber(Int_t chamberID){fBadChamberID[fNBadChambers++] = chamberID;}; private: @@ -52,8 +58,13 @@ private: Bool_t GetTRDchamberID(AliESDtrack * const track, Int_t *detectors); void SetChamberGain(); - void ApplyGainCorrection(AliESDtrack *track); + void ApplyGainCorrection(AliESDtrack *track, const Int_t * const detectors); + void ApplyRunByRunCorrection(AliESDtrack *const track); + void MaskChambers(AliESDtrack * const track, const Int_t * const detectors); void LoadReferences(); + void LoadDeadChambersFromCDB(); + void LoadRunByRunCorrection(const char *filename); + Bool_t IsBadChamber(Int_t chamberID); AliESDEvent *fESD; //! the ESD Event AliESDpid *fESDpid; //! ESD PID object @@ -62,16 +73,23 @@ private: AliTRDCalDet *fChamberGainNew; // New TRD Chamber Gain Factor AliTRDCalDet *fChamberVdriftOld; // Old drift velocity calibration AliTRDCalDet *fChamberVdriftNew; // New drift velocity calibration + AliOADBContainer *fRunByRunCorrection; // Run by run gain correction - TString fRefFilename; // path and name to the NNref file Int_t fPIDmethod; // PID method Double_t fNormalizationFactor; // dE/dx Normalization Factor Double_t fPthreshold; // Low Momentum threshold for calibration Int_t fBadChamberID[kNChambers]; // List of Bad Chambers + UInt_t fSlicesForPID[2]; // Select range of slices used in the PID response UInt_t fNBadChambers; // Number of bad chambers + const char *fGeoFile; // File with geometry.root Bool_t fGainCorrection; // Apply gain correction + Bool_t fLoadReferences; // Tender Load references Bool_t fLoadReferencesFromCDB; // Load References from CDB + Bool_t fLoadDeadChambers; // Load dead chambers Bool_t fHasReferences; // has references loaded + Bool_t fHasNewCalibration; // has new calibration + Bool_t fDebugMode; // Run in debug mode + TString fNameRunByRunCorrection; // filename with the run-by-run gain correction AliTRDTenderSupply(const AliTRDTenderSupply&c); AliTRDTenderSupply& operator= (const AliTRDTenderSupply&c); -- 2.39.3