From 8d14dc145334f7037124a5cd3386d1985bd376b1 Mon Sep 17 00:00:00 2001 From: schutz Date: Mon, 18 Dec 2006 11:31:33 +0000 Subject: [PATCH] New TRD QA task --- ESDCheck/AliTRDQATask.cxx | 620 ++++++++++++++++++++++++++++++++ ESDCheck/AliTRDQATask.h | 100 ++++++ ESDCheck/AnalysisCheckLinkDef.h | 2 + ESDCheck/libAnalysisCheck.pkg | 3 +- 4 files changed, 724 insertions(+), 1 deletion(-) create mode 100644 ESDCheck/AliTRDQATask.cxx create mode 100644 ESDCheck/AliTRDQATask.h diff --git a/ESDCheck/AliTRDQATask.cxx b/ESDCheck/AliTRDQATask.cxx new file mode 100644 index 00000000000..6bc7866b17f --- /dev/null +++ b/ESDCheck/AliTRDQATask.cxx @@ -0,0 +1,620 @@ +/************************************************************************** + * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. * + * * + * Author: The ALICE Off-line Project. * + * Contributors are mentioned in the code where appropriate. * + * * + * Permission to use, copy, modify and distribute this software and its * + * documentation strictly for non-commercial purposes is hereby granted * + * without fee, provided that the above copyright notice appears in all * + * copies and that both the copyright notice and this permission notice * + * appear in the supporting documentation. The authors make no claims * + * about the suitability of this software for any purpose. It is * + * provided "as is" without express or implied warranty. * + **************************************************************************/ +//_________________________________________________________________________ +// An analysis task to check the TRD data in simulated data +// +//*-- Sylwester Radomski +////////////////////////////////////////////////////////////////////////////// +// track type codding +// +// tpci = kTPCin +// tpco = kTPCout +// tpcz = kTPCout && !kTRDout +// trdo = kTRDout +// trdr = kTRDref +// trdz = kTRDout && !kTRDref +// + +#include "AliTRDQATask.h" + +#include "TChain.h" +#include "TH1D.h" +#include "TH2D.h" +#include "TFile.h" +#include "TStyle.h" +#include "TGaxis.h" +#include "TCanvas.h" + +#include "AliESD.h" +#include "AliLog.h" + +//______________________________________________________________________________ +AliTRDQATask::AliTRDQATask(const char *name) : + AliAnalysisTask(name,""), + fChain(0), + fESD(0) +{ + // Constructor. + // Input slot #0 works with an Ntuple + DefineInput(0, TChain::Class()); + // Output slot #0 writes into a TH1 container + DefineOutput(0, TObjArray::Class()) ; +} + +//______________________________________________________________________________ +void AliTRDQATask::Init(const Option_t *) +{ + // Initialisation of branch container and histograms + + AliInfo(Form("*** Initialization of %s", GetName())) ; + + // Get input data + fChain = dynamic_cast(GetInputData(0)) ; + if (!fChain) { + AliError(Form("Input 0 for %s not found\n", GetName())); + return ; + } + + if (!fESD) { + // One should first check if the branch address was taken by some other task + char ** address = (char **)GetBranchAddress(0, "ESD"); + if (address) { + fESD = (AliESD *)(*address); + AliInfo("Old ESD found."); + } + if (!fESD) { + fESD = new AliESD(); + SetBranchAddress(0, "ESD", &fESD); + if (fESD) AliInfo("ESD connected."); + } + } + // The output objects will be written to + TDirectory * cdir = gDirectory ; + OpenFile(0, Form("%s.root", GetName()), "RECREATE"); + if (cdir) + cdir->cd() ; + + // create histograms + + fNTracks = new TH1D("ntracks", ";number of all tracks", 500, -0.5, 499.5); + fEventSize = new TH1D("evSize", ";event size (MB)", 100, 0, 5); + + fTrackStatus = new TH1D("trackStatus", ";status bit", 32, -0.5, 31.5); + fKinkIndex = new TH1D("kinkIndex", ";kink index", 41, -20.5, 20.5); + + fParIn = new TH1D("parIn", "Inner Plane", 2, -0.5, 1.5); + fParOut = new TH1D("parOut", "Outer Plane", 2, -0.5, 1.5); + + fXIn = new TH1D("xIn", ";X at the inner plane (cm)", 200, 50, 250); + fXOut = new TH1D("xOut", ";X at the outer plane (cm)", 300, 50, 400); + + const int nNameAlpha = 4; + const char *namesAlpha[nNameAlpha] = {"alphaTPCi", "alphaTPCo", "alphaTRDo", "alphaTRDr"}; + //TH1D *fAlpha[4]; + for(int i=0; iClone(Form("eff_%s_%s", suf[0], suf[1])); + fEffPt[1] = (TH1D*) fPt[0]->Clone(Form("eff_%s_%s", suf[1], suf[3])); + fEffPt[2] = (TH1D*) fPt[0]->Clone(Form("eff_%s_%s", suf[3], suf[4])); + fEffPt[3] = (TH1D*) fPt[0]->Clone(Form("eff_%s_%s", suf[1], suf[4])); + + for(int i=0; i<4; i++) { + fEffPt[i]->Sumw2(); + fEffPt[i]->SetMarkerStyle(20); + fEffPt[i]->SetMinimum(0); + fEffPt[i]->SetMaximum(1.1); + } + + // track features + fClustersTRD[0] = new TH1D("clsTRDo", "TRDo;number of clusters", 130, -0.5, 129.5);; + fClustersTRD[1] = new TH1D("clsTRDr", "TRDr;number of clusters", 130, -0.5, 129.5);; + fClustersTRD[2] = new TH1D("clsTRDz", "TRDz;number of clusters", 130, -0.5, 129.5);; + + // for good refitted tracks only + fTime = new TH1D("time", ";time bin", 25, -0.5, 24.5); + fBudget = new TH1D("budget", ";material budget", 100, 0, 100); + fQuality = new TH1D("quality", ";track quality", 100, 0, 1.1); + fSignal = new TH1D("signal", ";signal", 100, 0, 1e3); + + // dEdX and PID + fTrdSigMom = new TH2D("trdSigMom", ";momentum (GeV/c);signal", 100, 0, 3, 100, 0, 1e3); + fTpcSigMom = new TH2D("tpcSigMom", ";momentum (GeV/c);signal", 100, 0, 3, 100, 0, 200); + + const char *pidName[6] = {"El", "Mu", "Pi", "K", "P", "Ch"}; + for(int i=0; i<6; i++) { + + // TPC + fTpcPID[i] = new TH1D(Form("tpcPid%s",pidName[i]), pidName[i], 100, 0, 1.5); + fTpcPID[i]->GetXaxis()->SetTitle("probability"); + + fTpcSigMomPID[i] = new TH2D(Form("tpcSigMom%s",pidName[i]), "", 100, 0, 3, 100, 0, 200); + fTpcSigMomPID[i]->SetTitle(Form("%s;momentum (GeV/c);signal",pidName[i])); + + // TRD + fTrdPID[i] = new TH1D(Form("trdPid%s",pidName[i]), pidName[i], 100, 0, 1.5); + fTrdPID[i]->GetXaxis()->SetTitle("probability"); + + fTrdSigMomPID[i] = new TH2D(Form("trdSigMom%s",pidName[i]), "", 100, 0, 3, 100, 0, 1e3); + fTrdSigMomPID[i]->SetTitle(Form("%s;momentum (GeV/c);signal",pidName[i])); + } + + + // create output container + fOutputContainer = new TObjArray(80); + + // register histograms to the container + TIter next(gDirectory->GetList()); + TObject *obj; + int counter = 0; + + while (obj = next.Next()) { + if (obj->InheritsFrom("TH1")) fOutputContainer->AddAt(obj, counter++); + } + + AliInfo(Form("Number of histograms = %d", counter)); + + } + +//______________________________________________________________________________ +void AliTRDQATask::Exec(Option_t *) +{ + // Process one event + + Long64_t entry = fChain->GetReadEntry() ; + + // Processing of one event + + if (!fESD) { + AliError("fESD is not connected to the input!") ; + return ; + } + + if ( !((entry-1)%100) ) + AliInfo(Form("%s ----> Processing event # %lld", (dynamic_cast(fChain))->GetFile()->GetName(), entry)) ; + + int nTracks = fESD->GetNumberOfTracks(); + fNTracks->Fill(nTracks); + + // track loop + for(int i=0; iGetTrack(i); + const AliExternalTrackParam *paramOut = track->GetOuterParam(); + const AliExternalTrackParam *paramIn = track->GetInnerParam(); + + fParIn->Fill(!!paramIn); + if (!paramIn) continue; + fXIn->Fill(paramIn->GetX()); + + fParOut->Fill(!!paramOut); + if (!paramOut) continue; + fXOut->Fill(paramOut->GetX()); + + int sector = GetSector(paramOut->GetAlpha()); + if (!CheckSector(sector)) continue; + fSectorTRD->Fill(sector); + + fKinkIndex->Fill(track->GetKinkIndex(0)); + if (track->GetKinkIndex(0)) continue; + + UInt_t u = 1; + UInt_t status = track->GetStatus(); + for(int bit=0; bit<32; bit++) + if (u<Fill(bit); + + const int nbits = 6; + int bit[6] = {0,0,0,0,0,0}; + bit[0] = status & AliESDtrack::kTPCin; + bit[1] = status & AliESDtrack::kTPCout; + bit[2] = (status & AliESDtrack::kTPCout) && !(status & AliESDtrack::kTRDout); + bit[3] = status & AliESDtrack::kTRDout; + bit[4] = status & AliESDtrack::kTRDrefit; + bit[5] = (status & AliESDtrack::kTRDout) && !(status & AliESDtrack::kTRDrefit); + + + // transverse momentum + const double *val = track->GetParameter(); // parameters at the vertex + double pt = 1./TMath::Abs(val[4]); + + for(int b=0; bFill(pt); + fTheta[b]->Fill(val[3]); + fSigmaY[b]->Fill(TMath::Sqrt(paramOut->GetSigmaY2())); + fChi2[b]->Fill(track->GetTRDchi2()/track->GetTRDncls()); + fPlaneYZ[b]->Fill(paramOut->GetY(), paramOut->GetZ()); + } + } + + // sectors + if (bit[1]) { + fAlpha[0]->Fill(paramIn->GetAlpha()); + fAlpha[1]->Fill(paramOut->GetAlpha()); + } + + if (bit[3]) fAlpha[2]->Fill(paramOut->GetAlpha()); + if (bit[4]) fAlpha[3]->Fill(paramOut->GetAlpha()); + + // clusters + for(int b=0; b<3; b++) + if (bit[3+b]) fClustersTRD[b]->Fill(track->GetTRDncls()); + + // refitted only + if (!bit[4]) continue; + + fQuality->Fill(track->GetTRDQuality()); + fBudget->Fill(track->GetTRDBudget()); + fSignal->Fill(track->GetTRDsignal()); + + fTrdSigMom->Fill(track->GetP(), track->GetTRDsignal()); + fTpcSigMom->Fill(track->GetP(), track->GetTPCsignal()); + + // PID only + if (status & AliESDtrack::kTRDpid) { + + for(int l=0; l<6; l++) fTime->Fill(track->GetTRDTimBin(l)); + + // fill pid histograms + double trdr0 = 0, tpcr0 = 0; + int trdBestPid = 5, tpcBestPid = 5; // charged + const double minPidValue = 0.9; + + double pp[5]; + track->GetTPCpid(pp); // ESD inconsequence + + for(int pid=0; pid<5; pid++) { + + trdr0 += track->GetTRDpid(pid); + tpcr0 += pp[pid]; + + fTrdPID[pid]->Fill(track->GetTRDpid(pid)); + fTpcPID[pid]->Fill(pp[pid]); + + if (track->GetTRDpid(pid) > minPidValue) trdBestPid = pid; + if (pp[pid] > minPidValue) tpcBestPid = pid; + } + + fTrdPID[5]->Fill(trdr0); // check unitarity + fTrdSigMomPID[trdBestPid]->Fill(track->GetP(), track->GetTRDsignal()); + + fTpcPID[5]->Fill(tpcr0); // check unitarity + fTpcSigMomPID[tpcBestPid]->Fill(track->GetP(), track->GetTPCsignal()); + } + + } + + CalculateEff(); + PostData(0, fOutputContainer); +} + +//______________________________________________________________________________ +void AliTRDQATask::Terminate(Option_t *) +{ + // Processing when the event loop is ended + AliInfo("TRD QA module"); + + // create efficiency histograms + + CalculateEff(); + PostData(0, fOutputContainer); + + DrawESD() ; + DrawGeoESD() ; + //DrawConvESD() ; + DrawPidESD() ; +} + +//______________________________________________________________________________ +int AliTRDQATask::GetSector(double alpha) +{ + // Gets the sector number + + double size = TMath::DegToRad() * 20.; + int sector = (int)((alpha + TMath::Pi())/size); + return sector; +} + +//______________________________________________________________________________ +int AliTRDQATask::CheckSector(int sector) +{ + // Checks the sector number + const int nSec = 8; + int sec[] = {2,3,5,6,11,12,13,15}; + + for(int i=0; iReset(); + + fEffPt[0]->Add(fPt[1]); + fEffPt[0]->Divide(fPt[0]); + + fEffPt[1]->Add(fPt[3]); + fEffPt[1]->Divide(fPt[1]); + + fEffPt[2]->Add(fPt[4]); + fEffPt[2]->Divide(fPt[3]); + + fEffPt[3]->Add(fPt[4]); + fEffPt[3]->Divide(fPt[1]); +} + +//______________________________________________________________________________ +void AliTRDQATask::DrawESD() +{ + // Makes a few plots + + AliInfo("Plotting....") ; + + gROOT->SetStyle("Plain"); + gStyle->SetPalette(1); + gStyle->SetOptStat(0); + + TGaxis::SetMaxDigits(3); + + gStyle->SetLabelFont(52, "XYZ"); + gStyle->SetTitleFont(62, "XYZ"); + gStyle->SetPadRightMargin(0.02); + + // draw all + + const int nplots = 18; + const int nover[nplots] = {1,1,1,4,1,1,1,1,1,1,2,1,1,3,1,1,1,1}; + const int nnames = 24; + const char *names[nnames] = { + "ntracks", "kinkIndex", "trackStatus", + "ptTPCi", "ptTPCo", "ptTRDo", "ptTRDr", "ptTPCz", "ptTRDz", + "eff_TPCi_TPCo", "eff_TPCo_TRDo", "eff_TRDo_TRDr", "eff_TPCo_TRDr", + "clsTRDo", "clsTRDr", "clsTRDz", + "alphaTPCi", "alphaTPCo", "alphaTRDo", "alphaTRDr", "sectorTRD", + "time", "budget", "signal" + }; + + const int logy[nnames] = { + 1,1,1, + 1,1,1, + 0,0,0,0, + 1,1, + 0,0,0,0,0, + 0,1,1 + }; + + int nhist=0; + for(int i=0; iSetLogy(logy[i]); + + for(int j=0; j(gDirectory->FindObject(names[nhist++])); + if (!hist) continue; + + if (strstr(hist->GetName(), "eff")) { + hist->SetMarkerStyle(20); + hist->SetMinimum(0); + hist->SetMaximum(1.2); + } + + if (!j) hist->Draw(); + else hist->Draw("SAME"); + } + + gPad->Print(Form("%s.gif", names[i])); + } +} + +//______________________________________________________________________________ +void AliTRDQATask::DrawGeoESD() +{ + // Makes a few plots + + AliInfo("Plotting....") ; + + gStyle->SetOptStat(0); + TGaxis::SetMaxDigits(3); + + gStyle->SetLabelFont(52, "XYZ"); + gStyle->SetTitleFont(62, "XYZ"); + + gStyle->SetPadTopMargin(0.06); + gStyle->SetTitleBorderSize(0); + + // draw all + const int nnames = 7; + const char *names[nnames] = { + "xIn", "xOut", + "planeYZTPCo", "planeYZTPCz", "planeYZTRDo", "planeYZTRDr", "planeYZTRDz", + }; + + const char *opt[nnames] = { + "", "", + "colz","colz", "colz", "colz", "colz" + }; + + const int logy[nnames] = { + 1,1, + 0,0,0,0,0 + }; + + for(int i=0; i(gDirectory->FindObject(names[i])); + if (!hist) continue; + + if (i<2) new TCanvas(names[i], names[i], 500, 300); + else new TCanvas(names[i], names[i], 300, 900); + + gPad->SetLogy(logy[i]); + if (strstr(opt[i],"colz")) gPad->SetRightMargin(0.1); + + hist->Draw(opt[i]); + AliInfo(Form("%s\t%d", names[i], hist->GetEntries())); + + gPad->Print(Form("%s.gif", names[i])); + } +} + +//______________________________________________________________________________ +void AliTRDQATask::DrawConvESD() +{ + // Makes a few plots + + AliInfo("Plotting....") ; + + gROOT->SetStyle("Plain"); + gROOT->ForceStyle(); + gStyle->SetPalette(1); + + TGaxis::SetMaxDigits(3); + + gStyle->SetLabelFont(52, "XYZ"); + gStyle->SetTitleFont(62, "XYZ"); + gStyle->SetPadRightMargin(0.02); + + const int nnames = 9; + const int nplots = 5; + const int nover[nplots] = {3,1,1,3,1}; + + const char *names[nnames] = { + "sigmaYTPCo","sigmaYTRDo", "sigmaYTRDr", "sigmaYTPCz", "sigmaYTRDz", + "Chi2TPCo", "Chi2TRDo", "Chi2TRDr", "Chi2TRDz" + }; + + const char *opt[nplots] = { + "", "", "","","", + }; + + const int logy[nplots] = { + 0,0,0,1,1 + }; + + int nhist = 0; + for(int i=0; iSetLogy(logy[i]); + if (strstr(opt[i],"colz")) gPad->SetRightMargin(0.1); + + for(int j=0; j(gDirectory->FindObject(names[nhist++])); + if (!j) hist->Draw(opt[i]); + else hist->Draw("same"); + } + + gPad->Print(Form("%s.eps", names[i])); + } +} + +//______________________________________________________________________________ +void AliTRDQATask::DrawPidESD() +{ + // Makes a few plots + + AliInfo("Plotting....") ; + + gROOT->SetStyle("Plain"); + gROOT->ForceStyle(); + gStyle->SetPalette(1); + gStyle->SetOptStat(0); + + TGaxis::SetMaxDigits(3); + + gStyle->SetLabelFont(52, "XYZ"); + gStyle->SetTitleFont(62, "XYZ"); + + gStyle->SetPadTopMargin(0.05); + gStyle->SetPadRightMargin(0.02); + + // draw all + + const int nnames = 27; + const char *names[nnames] = { + + "signal", "trdSigMom", "tpcSigMom", + + "trdPidEl", "trdPidMu", "trdPidPi", "trdPidK", "trdPidP", "trdPidCh", + "trdSigMomEl", "trdSigMomMu", "trdSigMomPi", "trdSigMomK", "trdSigMomP", "trdSigMomCh", + + "tpcPidEl", "tpcPidMu", "tpcPidPi", "tpcPidK", "tpcPidP", "tpcPidCh", + "tpcSigMomEl", "tpcSigMomMu", "tpcSigMomPi", "tpcSigMomK", "tpcSigMomP", "tpcSigMomCh" + + }; + + const char *opt[nnames] = { + + "", "colz", "colz", + + "", "", "", "", "", "" , + "colz", "colz", "colz", "colz", "colz", "colz", + + "", "", "", "", "", "" , + "colz", "colz", "colz", "colz", "colz", "colz" + }; + + const int logy[nnames] = { + + 0,0,0, + + 1,1,1,1,1,1, + 0,0,0,0,0,0, + + 1,1,1,1,1,1, + 0,0,0,0,0,0 + }; + + for(int i=0; i(gDirectory->FindObject(names[i])); + if (!hist) continue; + + new TCanvas(names[i], names[i], 500, 300); + gPad->SetLogy(logy[i]); + if (strstr(opt[i],"colz")) gPad->SetRightMargin(0.1); + + if (strstr(names[i],"sigMom")) gPad->SetLogz(1); + if (strstr(names[i],"SigMom")) gPad->SetLogz(1); + + hist->Draw(opt[i]); + AliInfo(Form("%s\t%d", names[i], hist->GetEntries())); + + gPad->Print(Form("%s.gif", names[i])); + } + +} diff --git a/ESDCheck/AliTRDQATask.h b/ESDCheck/AliTRDQATask.h new file mode 100644 index 00000000000..42181b6c1dc --- /dev/null +++ b/ESDCheck/AliTRDQATask.h @@ -0,0 +1,100 @@ +#ifndef ALITRDQATASK_H +/* Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. * + * See cxx source for full Copyright notice */ +//______________________________________________________________________________ +// An analysis task to check the TRD data in simulated data +// +//*-- Sylwester Radomski +////////////////////////////////////////////////////////////////////////////// + +#include +#include "AliAnalysisTask.h" + +class AliESD; +class TH1D; +class TH2D; + +class AliTRDQATask : public AliAnalysisTask { + +public: + AliTRDQATask(const char *name); + virtual ~AliTRDQATask() {} + + virtual void Exec(Option_t * opt = ""); + virtual void Init(Option_t * opt = ""); + virtual void Terminate(Option_t * opt = ""); + +private: + + int GetSector(double alpha); + int CheckSector(int sector); + void CalculateEff(); + void DrawESD() ; + void DrawGeoESD() ; + void DrawConvESD() ; + void DrawPidESD() ; + + TTree * fChain; //!pointer to the analyzed TTree or TChain + AliESD * fESD; //! Declaration of leave types + + TObjArray * fOutputContainer; //! output data container + + // options + int fConfSM; + + // Histograms + TH1D *fNTracks; + TH1D *fEventSize; + TH1D *fTrackStatus; + + TH1D *fParIn; + TH1D *fParOut; + TH1D *fKinkIndex; + + // TPC clusters histograms + //TH1D *fTpcNCls; + //TH1D *fTpcFCls; + //TH1D *fTpcRCls; + + // last measurement X plane + TH1D *fXIn; + TH1D *fXOut; + + // sector + TH1D *fAlpha[4]; + TH1D *fSectorTRD; + + //static const int knbits = 5; + + // track parameters + TH1D *fPt[6]; + TH1D *fTheta[6]; + TH1D *fSigmaY[6]; + TH1D *fChi2[6]; + TH2D *fPlaneYZ[6]; + + TH1D *fEffPt[4]; + + // track features + TH1D *fClustersTRD[3]; + + // for good refitted tracks only + TH1D *fTime; + TH1D *fBudget; + TH1D *fQuality; + TH1D *fSignal; + + // PID for TPC and TRD + TH2D *fTrdSigMom; + TH2D *fTpcSigMom; + + TH1D *fTrdPID[6]; + TH2D *fTrdSigMomPID[6]; + + TH1D *fTpcPID[6]; + TH2D *fTpcSigMomPID[6]; + + + ClassDef(AliTRDQATask, 0); // a TRD analysis task +}; +#endif // ALITRDQATASK_H diff --git a/ESDCheck/AnalysisCheckLinkDef.h b/ESDCheck/AnalysisCheckLinkDef.h index 9e893a15841..c84ed6f3d2f 100644 --- a/ESDCheck/AnalysisCheckLinkDef.h +++ b/ESDCheck/AnalysisCheckLinkDef.h @@ -22,4 +22,6 @@ #pragma link C++ class AliFMDQATask+; +#pragma link C++ class AliTRDQATask+; + #endif diff --git a/ESDCheck/libAnalysisCheck.pkg b/ESDCheck/libAnalysisCheck.pkg index 70b6e4f2316..8a8764dbb20 100644 --- a/ESDCheck/libAnalysisCheck.pkg +++ b/ESDCheck/libAnalysisCheck.pkg @@ -6,7 +6,8 @@ SRCS = AliAnalysisGoodies.cxx \ AliHMPIDQATask.cxx \ AliT0QATask.cxx \ AliMUONQATask.cxx \ - AliFMDQATask.cxx + AliFMDQATask.cxx \ + AliTRDQATask.cxx HDRS:= $(SRCS:.cxx=.h) -- 2.43.0