#include "TPDGCode.h"
#include "TF1.h"
#include "TH1F.h"
+#include "TH1D.h"
#include "TH2F.h"
+#include "TProfile.h"
+#include "TProfile2D.h"
#include "TGraph.h"
#include "TGraphErrors.h"
#include <TObjArray.h>
#include <TList.h>
-#include "AliPID.h"
+// #include "AliPID.h"
#include "AliESDEvent.h"
#include "AliESDInputHandler.h"
#include "AliTrackReference.h"
#include "AliAnalysisTask.h"
-#include "AliAnalysisManager.h"
+// #include "AliAnalysisManager.h"
+#include "AliTRDtrackerV1.h"
#include "AliTRDtrackV1.h"
+#include "AliTRDcluster.h"
#include "AliTRDReconstructor.h"
#include "AliCDBManager.h"
-#include "../Cal/AliTRDCalPID.h"
+// #include "../Cal/AliTRDCalPID.h"
+#include "AliTRDpidUtil.h"
#include "AliTRDpidChecker.h"
ClassImp(AliTRDpidChecker)
//________________________________________________________________________
-AliTRDpidChecker::AliTRDpidChecker(const char *name)
- :AliAnalysisTask(name, "")
- ,fObjectContainer(0x0)
- ,fTracks(0x0)
+AliTRDpidChecker::AliTRDpidChecker()
+ :AliTRDrecoTask("PID", "PID Checker")
,fReconstructor(0x0)
-// ,fDebugLevel(1)
-// ,fDebugStream(0x0)
{
//
// Default constructor
fReconstructor = new AliTRDReconstructor();
fReconstructor->SetRecoParam(AliTRDrecoParam::GetLowFluxParam());
-
- DefineInput(0, TObjArray::Class());
- DefineOutput(0, TObjArray::Class());
}
//________________________________________________________________________
AliTRDpidChecker::~AliTRDpidChecker()
{
- if(fObjectContainer){
- //fObjectContainer->Delete();
- //delete fObjectContainer;
- }
-}
-
-
-
-//________________________________________________________________________
-void AliTRDpidChecker::ConnectInputData(Option_t *)
-{
- //
- // Connect input data
- //
-
- fTracks = dynamic_cast<TObjArray*>(GetInputData(0));
+ if(fReconstructor) delete fReconstructor;
}
// Called once
OpenFile(0, "RECREATE");
- fObjectContainer = new TObjArray();
- fObjectContainer->Add(new TH1F("hMom", "momentum distribution", AliTRDCalPID::kNMom, 0.5, 11.5));
+ Int_t xBins = AliPID::kSPECIES*AliTRDCalPID::kNMom;
+ fContainer = new TObjArray();
+ fContainer -> Expand(kGraphNN + 1);
+ const Float_t epsilon = 1./(2*(AliTRDpidUtil::kBins-1)); // get nice histos with bin center at 0 and 1
// histos of the electron probability of all 5 particle species and 11 momenta for the 2-dim LQ method
- const Int_t kBins = 12000; // binning of the histos and eficiency calculation
- for(Int_t iPart = 0; iPart < AliPID::kSPECIES; iPart++){
- for(Int_t iMom = 0; iMom < AliTRDCalPID::kNMom; iMom++){
- fObjectContainer->Add(new TH1F(Form("PID%d_%d_LQ",iPart,iMom), Form("PID distribution for %d_%d", iPart, iMom), kBins, -0.1, 1.1));
- }
- }
+ fContainer->AddAt(
+ new TH2F("PID_LQ", "",
+ xBins, -0.5, xBins - 0.5,
+ AliTRDpidUtil::kBins, 0.-epsilon, 1.+epsilon)
+ ,kLQlikelihood);
- // histos of the electron probability of all 5 particle species and 11 momenta for the neural network method
- for(Int_t iPart = 0; iPart < AliPID::kSPECIES; iPart++){
- for(Int_t iMom = 0; iMom < AliTRDCalPID::kNMom; iMom++){
- fObjectContainer->Add(new TH1F(Form("PID%d_%d_NN",iPart,iMom), Form("PID distribution for %d_%d", iPart, iMom), kBins, -0.1, 1.1));
- }
- }
- // frame and TGraph of the pion efficiencies
- fObjectContainer -> Add(new TH2F("hFrame", "", 10, 0.4, 12., 10, 0.0005, 0.1));
- fObjectContainer -> Add(new TGraph(AliTRDCalPID::kNMom));
- fObjectContainer -> Add(new TGraphErrors(AliTRDCalPID::kNMom));
- fObjectContainer -> Add(new TGraph(AliTRDCalPID::kNMom));
- fObjectContainer -> Add(new TGraphErrors(AliTRDCalPID::kNMom));
+ // histos of the electron probability of all 5 particle species and 11 momenta for the neural network method
+ fContainer->AddAt(
+ new TH2F("PID_NN", "",
+ xBins, -0.5, xBins - 0.5,
+ AliTRDpidUtil::kBins, 0.-epsilon, 1.+epsilon)
+ ,kNNlikelihood);
+
+ // histos of the dE/dx distribution for all 5 particle species and 11 momenta
+ fContainer->AddAt(
+ new TH2F("dEdx", "",
+ xBins, -0.5, xBins - 0.5,
+ 200, 0, 10000)
+ ,kdEdx);
+
+ // histos of the pulse height distribution for all 5 particle species and 11 momenta
+ fContainer->AddAt(
+ new TProfile2D("PH", "",
+ xBins, -0.5, xBins - 0.5,
+ AliTRDtrackerV1::GetNTimeBins(), -0.5, AliTRDtrackerV1::GetNTimeBins() - 0.5)
+ ,kPH);
+
+
+ // momentum distributions - absolute and in momentum bins
+ fContainer->AddAt(new TH1F("hMom", "momentum distribution", 100, 0., 12.),kMomentum);
+ fContainer->AddAt(new TH1F("hMomBin", "momentum distribution in momentum bins", AliTRDCalPID::kNMom, 0.5, 11.5),kMomentumBin);
+
+
+ // TGraph of the pion efficiencies
+
+ TGraphErrors *gEffisLQ = 0x0;
+ TGraphErrors *gEffisNN = 0x0;
+
+ fContainer->AddAt(gEffisLQ = new TGraphErrors(), kGraphLQ);
+ gEffisLQ->SetLineColor(kBlue);
+ gEffisLQ->SetMarkerColor(kBlue);
+ gEffisLQ->SetMarkerStyle(29);
+
+ fContainer -> AddAt(gEffisNN = new TGraphErrors(),kGraphNN);
+ gEffisNN->SetLineColor(kRed);
+ gEffisNN->SetMarkerColor(kRed);
+ gEffisNN->SetMarkerStyle(29);
+
+ gEffisLQ -> SetNameTitle("gEffisLQErr", "Efficiencies and Errors of the 2-dim LQ method");
+ gEffisNN -> SetNameTitle("gEffisNNErr", "Efficiencies and Errors of the NN method");
+
}
//________________________________________________________________________
// AliTracker::SetFieldMap(field, kTRUE);
// }
- TH1F *hMom = (TH1F*)fObjectContainer->UncheckedAt(0);
- TH1F *hPIDLQ[AliPID::kSPECIES][AliTRDCalPID::kNMom];
- TH1F *hPIDNN[AliPID::kSPECIES][AliTRDCalPID::kNMom];
+ TH2F *hPIDLQ;
+ TH2F *hPIDNN;
+ TH2F *hdEdx;
+ TProfile2D *hPH;
- for(Int_t iPart = 0; iPart < AliPID::kSPECIES; iPart++){
- for(Int_t iMom = 0; iMom < AliTRDCalPID::kNMom; iMom++){
- hPIDLQ[iPart][iMom] = (TH1F*)fObjectContainer->At(iPart*AliTRDCalPID::kNMom+iMom+1);
- hPIDNN[iPart][iMom] = (TH1F*)fObjectContainer->At(iPart*AliTRDCalPID::kNMom+iMom+1+AliPID::kSPECIES*AliTRDCalPID::kNMom);
- }
- }
-
+ hPIDLQ = (TH2F*)fContainer->At(kLQlikelihood);
+ hPIDNN = (TH2F*)fContainer->At(kNNlikelihood);
+ hdEdx = (TH2F*)fContainer->At(kdEdx);
+ hPH = (TProfile2D*)fContainer->At(kPH);
+
+ TH1F *hMom = (TH1F*)fContainer->At(kMomentum);
+ TH1F *hMomBin = (TH1F*)fContainer->At(kMomentumBin);
+
Int_t labelsacc[10000];
memset(labelsacc, 0, sizeof(Int_t) * 10000);
Float_t mom;
ULong_t status;
Int_t nTRD = 0;
-
+ Float_t *fdEdx;
AliTRDtrackInfo *track = 0x0;
- AliTRDtrackV1 *TRDtrack = 0x0;
+ AliTRDtrackV1 *TRDtrack = 0x0;
AliTrackReference *ref = 0x0;
AliExternalTrackParam *esd = 0x0;
+
+ AliTRDseedV1 *TRDtracklet[AliTRDgeometry::kNlayer];
+ for(Int_t iChamb = 0; iChamb < AliTRDgeometry::kNlayer; iChamb++)
+ TRDtracklet[iChamb] = 0x0;
+
+ AliTRDcluster *TRDcluster = 0x0;
+
+ AliTRDpidUtil *util = new AliTRDpidUtil();
for(Int_t itrk=0; itrk<fTracks->GetEntriesFast(); itrk++){
track = (AliTRDtrackInfo*)fTracks->UncheckedAt(itrk);
if(!track->HasESDtrack()) continue;
if(!(TRDtrack = track->GetTRDtrack())) continue;
//&&(track->GetNumberOfClustersRefit()
- // use only tracks taht hit 6 chambers
- if(!(TRDtrack->GetNumberOfTracklets() == AliTRDCalPID::kNPlane)) continue;
+
+ // use only tracks that hit 6 chambers
+ if(!(TRDtrack->GetNumberOfTracklets() == AliTRDgeometry::kNlayer)) continue;
ref = track->GetTrackRef(0);
esd = track->GetOuterParam();
labelsacc[nTRD] = track->GetLabel();
nTRD++;
- // fill momentum histo to have the momentum distribution
- hMom -> Fill(mom);
-
// set the 11 momentum bins
Int_t iMomBin = -1;
- if(mom < .7) iMomBin = 0;
- else if(mom < .9) iMomBin = 1;
- else if(mom < 1.25) iMomBin = 2;
- else if(mom < 1.75) iMomBin = 3;
- else if(mom < 2.5) iMomBin = 4;
- else if(mom < 3.5) iMomBin = 5;
- else if(mom < 4.5) iMomBin = 6;
- else if(mom < 5.5) iMomBin = 7;
- else if(mom < 7.0) iMomBin = 8;
- else if(mom < 9.0) iMomBin = 9;
- else iMomBin = 10;
+ iMomBin = util -> GetMomentumBin(mom);
+ if(fDebugLevel>=4) Printf("MomBin[%d] MomTot[%f]", iMomBin, mom);
+
+
+ // fill momentum histo to have the momentum distribution
+ hMom -> Fill(mom);
+ hMomBin -> Fill(iMomBin);
+
// set the reconstructor
TRDtrack -> SetReconstructor(fReconstructor);
- // calculate the probabilities and fill histograms for electrons using 2-dim LQ
+ // if no monte carlo data available -> use TRDpid
+ if(!HasMCdata()){
+ fReconstructor -> SetOption("nn");
+ TRDtrack -> CookPID();
+ if(TRDtrack -> GetPID(0) > TRDtrack -> GetPID(1) + TRDtrack -> GetPID(2) + TRDtrack -> GetPID(3) + TRDtrack -> GetPID(4)){
+ track -> SetPDG(kElectron);
+ }
+ else if(TRDtrack -> GetPID(4) > TRDtrack -> GetPID(2) && TRDtrack -> GetPID(4) > TRDtrack -> GetPID(3) && TRDtrack -> GetPID(4) > TRDtrack -> GetPID(1)){
+ track -> SetPDG(kProton);
+ }
+ else if(TRDtrack -> GetPID(3) > TRDtrack -> GetPID(1) && TRDtrack -> GetPID(3) > TRDtrack -> GetPID(2)){
+ track -> SetPDG(kKPlus);
+ }
+ else if(TRDtrack -> GetPID(1) > TRDtrack -> GetPID(2)){
+ track -> SetPDG(kMuonPlus);
+ }
+ else{
+ track -> SetPDG(kPiPlus);
+ }
+ }
+
+
+ // calculate the probabilities for electron probability using 2-dim LQ, the deposited charge per chamber and the pulse height spectra and fill histograms
fReconstructor -> SetOption("!nn");
TRDtrack -> CookPID();
+ if(fDebugLevel>=4) Printf("PIDmethod[%d] Slices[%d] PDG[%d] LQLike[%f]", fReconstructor->GetPIDMethod(), fReconstructor->GetNdEdxSlices(), track->GetPDG(), TRDtrack -> GetPID(0));
+
+
+ Float_t SumdEdx[AliTRDgeometry::kNlayer];
+ for(Int_t iChamb = 0; iChamb < AliTRDgeometry::kNlayer; iChamb++){
+ TRDtracklet[iChamb] = TRDtrack -> GetTracklet(iChamb);
+ SumdEdx[iChamb] = 0.;
+ fdEdx = TRDtracklet[iChamb] -> GetdEdx();
+ SumdEdx[iChamb] += fdEdx[0] + fdEdx[1] + fdEdx[2];
+ }
+
switch(track->GetPDG()){
case kElectron:
case kPositron:
- hPIDLQ[AliPID::kElectron][iMomBin] -> Fill(TRDtrack -> GetPID(0));
+ hPIDLQ -> Fill(AliPID::kElectron * AliTRDCalPID::kNMom + iMomBin, TRDtrack -> GetPID(0));
+ for(Int_t iChamb = 0; iChamb < AliTRDgeometry::kNlayer; iChamb++){
+ hdEdx -> Fill(AliPID::kElectron * AliTRDCalPID::kNMom + iMomBin, SumdEdx[iChamb]);
+ for(Int_t iClus = 0; iClus < AliTRDtrackerV1::GetNTimeBins(); iClus++){
+ if(!(TRDcluster = (AliTRDcluster*)TRDtracklet[iChamb] -> GetClusters(iClus)))
+ continue;
+ hPH -> Fill(AliPID::kElectron * AliTRDCalPID::kNMom + iMomBin, TRDcluster -> GetLocalTimeBin(), TRDtracklet[iChamb] -> GetdQdl(iClus));
+ }
+ }
break;
case kMuonPlus:
case kMuonMinus:
- hPIDLQ[AliPID::kMuon][iMomBin] -> Fill(TRDtrack -> GetPID(0));
+ hPIDLQ -> Fill(AliPID::kMuon * AliTRDCalPID::kNMom + iMomBin, TRDtrack -> GetPID(0));
+ for(Int_t iChamb = 0; iChamb < AliTRDgeometry::kNlayer; iChamb++){
+ hdEdx -> Fill(AliPID::kMuon * AliTRDCalPID::kNMom + iMomBin, SumdEdx[iChamb]);
+ for(Int_t iClus = 0; iClus < AliTRDtrackerV1::GetNTimeBins(); iClus++){
+ if(!(TRDcluster = (AliTRDcluster*)TRDtracklet[iChamb] -> GetClusters(iClus)))
+ continue;
+ hPH -> Fill(AliPID::kMuon * AliTRDCalPID::kNMom + iMomBin, TRDcluster -> GetLocalTimeBin(), TRDtracklet[iChamb] -> GetdQdl(iClus));
+ }
+ }
break;
case kPiPlus:
case kPiMinus:
- hPIDLQ[AliPID::kPion][iMomBin] -> Fill(TRDtrack -> GetPID(0));
+ hPIDLQ -> Fill(AliPID::kPion * AliTRDCalPID::kNMom + iMomBin, TRDtrack -> GetPID(0));
+ for(Int_t iChamb = 0; iChamb < AliTRDgeometry::kNlayer; iChamb++){
+ hdEdx -> Fill(AliPID::kPion * AliTRDCalPID::kNMom + iMomBin, SumdEdx[iChamb]);
+ for(Int_t iClus = 0; iClus < AliTRDtrackerV1::GetNTimeBins(); iClus++){
+ if(!(TRDcluster = (AliTRDcluster*)TRDtracklet[iChamb] -> GetClusters(iClus)))
+ continue;
+ hPH -> Fill(AliPID::kPion * AliTRDCalPID::kNMom + iMomBin, TRDcluster -> GetLocalTimeBin(), TRDtracklet[iChamb] -> GetdQdl(iClus));
+ }
+ }
break;
case kKPlus:
case kKMinus:
- hPIDLQ[AliPID::kKaon][iMomBin] -> Fill(TRDtrack -> GetPID(0));
+ hPIDLQ -> Fill(AliPID::kKaon * AliTRDCalPID::kNMom + iMomBin, TRDtrack -> GetPID(0));
+ for(Int_t iChamb = 0; iChamb < AliTRDgeometry::kNlayer; iChamb++){
+ hdEdx -> Fill(AliPID::kKaon * AliTRDCalPID::kNMom + iMomBin, SumdEdx[iChamb]);
+ for(Int_t iClus = 0; iClus < AliTRDtrackerV1::GetNTimeBins(); iClus++){
+ if(!(TRDcluster = (AliTRDcluster*)TRDtracklet[iChamb] -> GetClusters(iClus)))
+ continue;
+ hPH -> Fill(AliPID::kKaon * AliTRDCalPID::kNMom + iMomBin, TRDcluster -> GetLocalTimeBin(), TRDtracklet[iChamb] -> GetdQdl(iClus));
+ }
+ }
break;
case kProton:
case kProtonBar:
- hPIDLQ[AliPID::kMuon][iMomBin] -> Fill(TRDtrack -> GetPID(0));
+ hPIDLQ -> Fill(AliPID::kProton * AliTRDCalPID::kNMom + iMomBin, TRDtrack -> GetPID(0));
+ for(Int_t iChamb = 0; iChamb < AliTRDgeometry::kNlayer; iChamb++){
+ hdEdx -> Fill(AliPID::kProton * AliTRDCalPID::kNMom + iMomBin, SumdEdx[iChamb]);
+ for(Int_t iClus = 0; iClus < AliTRDtrackerV1::GetNTimeBins(); iClus++){
+ if(!(TRDcluster = (AliTRDcluster*)TRDtracklet[iChamb] -> GetClusters(iClus)))
+ continue;
+ hPH -> Fill(AliPID::kProton * AliTRDCalPID::kNMom + iMomBin, TRDcluster -> GetLocalTimeBin(), TRDtracklet[iChamb] -> GetdQdl(iClus));
+ }
+ }
break;
}
// calculate the probabilities and fill histograms for electrons using NN
fReconstructor -> SetOption("nn");
TRDtrack->CookPID();
+
+
+ if(fDebugLevel>=4) Printf("PIDmethod[%d] Slices[%d] PDG[%d] NNLike[%f]", fReconstructor->GetPIDMethod(), fReconstructor->GetNdEdxSlices(), track->GetPDG(), TRDtrack -> GetPID(0));
+
+
switch(track->GetPDG()){
case kElectron:
case kPositron:
- hPIDNN[AliPID::kElectron][iMomBin] -> Fill(TRDtrack -> GetPID(0));
+ hPIDNN -> Fill(AliPID::kElectron * AliTRDCalPID::kNMom + iMomBin, TRDtrack -> GetPID(0));
break;
case kMuonPlus:
case kMuonMinus:
- hPIDNN[AliPID::kMuon][iMomBin] -> Fill(TRDtrack -> GetPID(0));
+ hPIDNN -> Fill(AliPID::kMuon * AliTRDCalPID::kNMom + iMomBin, TRDtrack -> GetPID(0));
break;
case kPiPlus:
case kPiMinus:
- hPIDNN[AliPID::kPion][iMomBin] -> Fill(TRDtrack -> GetPID(0));
+ hPIDNN -> Fill(AliPID::kPion * AliTRDCalPID::kNMom + iMomBin, TRDtrack -> GetPID(0));
break;
case kKPlus:
case kKMinus:
- hPIDNN[AliPID::kKaon][iMomBin] -> Fill(TRDtrack -> GetPID(0));
+ hPIDNN -> Fill(AliPID::kKaon * AliTRDCalPID::kNMom + iMomBin, TRDtrack -> GetPID(0));
break;
case kProton:
case kProtonBar:
- hPIDNN[AliPID::kMuon][iMomBin] -> Fill(TRDtrack -> GetPID(0));
+ hPIDNN -> Fill(AliPID::kProton * AliTRDCalPID::kNMom + iMomBin, TRDtrack -> GetPID(0));
break;
}
+
}
- PostData(0, fObjectContainer);
+ util -> Delete();
+ PostData(0, fContainer);
}
+//________________________________________________________
+void AliTRDpidChecker::GetRefFigure(Int_t ifig, Int_t &first, Int_t &last, Option_t *opt)
+{
+ opt = "pl";
+ switch(ifig){
+ case 0:
+ first = (Int_t)kGraphStart; last = first+2;
+ break;
+ default:
+ first = (Int_t)kGraphStart; last = first+2;
+ break;
+ }
+}
+
+
//________________________________________________________________________
-void AliTRDpidChecker::Terminate(Option_t *)
+Bool_t AliTRDpidChecker::PostProcess()
{
// Draw result to the screen
// Called once at the end of the query
- fObjectContainer = dynamic_cast<TObjArray*>(GetOutputData(0));
- if (!fObjectContainer) {
+ if (!fContainer) {
Printf("ERROR: list not available");
- return;
+ return kFALSE;
}
+// return kTRUE; // testing protection
+
-
// container for the pion efficiencies and the errors
Double_t PionEffiLQ[AliTRDCalPID::kNMom],
- PionEffiNN[AliTRDCalPID::kNMom],
PionEffiErrorLQ[AliTRDCalPID::kNMom],
- PionEffiErrorNN[AliTRDCalPID::kNMom];
+ EleEffiLQ[AliTRDCalPID::kNMom],
+ ThresholdLQ[AliTRDCalPID::kNMom];
+
+ Double_t PionEffiNN[AliTRDCalPID::kNMom],
+ PionEffiErrorNN[AliTRDCalPID::kNMom],
+ EleEffiNN[AliTRDCalPID::kNMom],
+ ThresholdNN[AliTRDCalPID::kNMom];
+ Float_t mom = 0.;
+
+ TH1D *Histo1=0x0, *Histo2=0x0;
+
+ TH2F *hPIDLQ=0x0, *hPIDNN=0x0;
+ hPIDLQ = (TH2F*)fContainer->At(kLQlikelihood);
+ hPIDNN = (TH2F*)fContainer->At(kNNlikelihood);
// calculate the pion efficiencies and the errors for 90% electron efficiency (2-dim LQ)
for(Int_t iMom = 0; iMom < AliTRDCalPID::kNMom; iMom++){
- PionEffiLQ[iMom] = GetPionEfficiency(iMom+1,iMom+23);
- PionEffiErrorLQ[iMom] = GetError(iMom+1,iMom+23);
- Printf("Pion Efficiency for 2-dim LQ is : %f +/- %f\n\n", PionEffiLQ[iMom], PionEffiErrorLQ[iMom]);
- }
-
- // calculate the pion efficiencies and the errors for 90% electron efficiency (NN)
- for(Int_t iMom = 0; iMom < AliTRDCalPID::kNMom; iMom++){
- PionEffiNN[iMom] = GetPionEfficiency(iMom+56,iMom+78);
- PionEffiErrorNN[iMom] = GetError(iMom+56,iMom+78);
- Printf("Pion Efficiency for NN is : %f +/- %f\n\n", PionEffiNN[iMom], PionEffiErrorNN[iMom]);
- }
-
- // create TGraph to plot the pion efficiencies
- TGraph *gEffisLQ=0x0, *gEffisNN=0x0;
- TGraphErrors *gEffisErrLQ=0x0, *gEffisErrNN=0x0;
+ AliTRDpidUtil *util = new AliTRDpidUtil();
+ mom = AliTRDCalPID::GetMomentum(iMom);
- gEffisLQ = (TGraph*)fObjectContainer->At(112);
- gEffisErrLQ = (TGraphErrors*)fObjectContainer->At(113);
- gEffisNN = (TGraph*)fObjectContainer->At(114);
- gEffisErrNN = (TGraphErrors*)fObjectContainer->At(115);
+ Histo1 = hPIDLQ -> ProjectionY("LQ_ele",AliTRDCalPID::kNMom*AliPID::kElectron+iMom+1,AliTRDCalPID::kNMom*AliPID::kElectron+iMom+1);
+ Histo2 = hPIDLQ -> ProjectionY("LQ_pio",AliTRDCalPID::kNMom*AliPID::kPion+iMom+1,AliTRDCalPID::kNMom*AliPID::kPion+iMom+1);
- for(Int_t iBin = 0; iBin < AliTRDCalPID::kNMom; iBin++){
- Float_t momentum = AliTRDCalPID::GetMomentum(iBin);
- gEffisLQ->SetPoint(iBin, momentum, PionEffiLQ[iBin]);
- gEffisErrLQ->SetPoint(iBin, momentum, PionEffiLQ[iBin]);
- gEffisErrLQ->SetPointError(iBin, 0., PionEffiErrorLQ[iBin]);
+ util -> CalculatePionEffi(Histo1, Histo2);
- gEffisNN->SetPoint(iBin, momentum, PionEffiNN[iBin]);
- gEffisErrNN->SetPoint(iBin, momentum, PionEffiNN[iBin]);
- gEffisErrNN->SetPointError(iBin, 0., PionEffiErrorNN[iBin]);
- }
+ PionEffiLQ[iMom] = util -> GetPionEfficiency();
+ PionEffiErrorLQ[iMom] = util -> GetError();
+ EleEffiLQ[iMom] = util -> GetCalcElectronEfficiency();
+ ThresholdLQ[iMom] = util -> GetThreshold();
- gEffisLQ -> SetNameTitle("gEffisLQ", "Efficiencies of the 2-dim LQ method");
- gEffisErrLQ -> SetNameTitle("gEffisLQErr", "Efficiencies and Errors of the 2-dim LQ method");
- gEffisNN -> SetNameTitle("gEffisNN", "Efficiencies of the NN method");
- gEffisErrNN -> SetNameTitle("gEffisNNErr", "Efficiencies and Errors of the NN method");
-}
+ if(fDebugLevel>=1) Printf("Pion Efficiency for 2-dim LQ is : %f +/- %f\n\n", PionEffiLQ[iMom], PionEffiErrorLQ[iMom]);
+ util -> Delete();
+ }
+
-//________________________________________________________________________
-Double_t AliTRDpidChecker::GetPionEfficiency(Int_t Index1, Int_t Index2){
+ // calculate the pion efficiencies and the errors for 90% electron efficiency (NN)
+ for(Int_t iMom = 0; iMom < AliTRDCalPID::kNMom; iMom++){
- Float_t Multi = 0.9; // electron efficiency
- Int_t abin, bbin;
- Double_t SumElecs, SumPions; // integrated sum of elecs and pions in the histos
- Double_t aBinSum, bBinSum; // content of a single bin in the histos
-
- TH1F *Histo1 = (TH1F*)fObjectContainer->At(Index1); // electron histo
- TH1F *Histo2 = (TH1F*)fObjectContainer->At(Index2); // pion histo
+ AliTRDpidUtil *util = new AliTRDpidUtil();
+ mom = AliTRDCalPID::GetMomentum(iMom);
+ Histo1 = hPIDNN -> ProjectionY("NN_ele",AliTRDCalPID::kNMom*AliPID::kElectron+iMom+1,AliTRDCalPID::kNMom*AliPID::kElectron+iMom+1);
+ Histo2 = hPIDNN -> ProjectionY("NN_pio",AliTRDCalPID::kNMom*AliPID::kPion+iMom+1,AliTRDCalPID::kNMom*AliPID::kPion+iMom+1);
- SumElecs = 0.;
- if(!(Histo1 -> GetEntries() > 0 && Histo2 -> GetEntries() > 0)){
- Printf("Warning: Histo momentum intervall %d has no Entries!", Index1-1);
- return -1.;
- }
- Histo1 -> Scale(1./Histo1->GetEntries());
- Histo2 -> Scale(1./Histo2->GetEntries());
+ util -> CalculatePionEffi(Histo1, Histo2);
+ PionEffiNN[iMom] = util -> GetPionEfficiency();
+ PionEffiErrorNN[iMom] = util -> GetError();
+ EleEffiNN[iMom] = util -> GetCalcElectronEfficiency();
+ ThresholdNN[iMom] = util -> GetThreshold();
- // calculate threshold for pion efficiency
- for (abin = (Histo1 -> GetNbinsX()); abin >= 0; abin--){
- aBinSum = 0;
- aBinSum = Histo1 -> GetBinContent(abin);
- if(!(aBinSum == 0)){
- SumElecs = SumElecs + aBinSum;
- }
- if (SumElecs >= Multi){
- break;
- }
- }
+ if(fDebugLevel>=1) Printf("Pion Efficiency for NN is : %f +/- %f\n\n", PionEffiNN[iMom], PionEffiErrorNN[iMom]);
-
- // calculate pion efficiency
- SumPions = 0.;
- for (bbin = (Histo2 -> GetNbinsX()); bbin >= abin; bbin--){
- bBinSum = 0;
- bBinSum = Histo2 -> GetBinContent(bbin);
- if(!(bBinSum == 0)){
- SumPions = SumPions + bBinSum;
- }
+ util -> Delete();
}
-
- // print the electron efficiency and its cuts
- Printf("Cut for momentum intervall %d and electron efficiency of %f for: 0.%d", Index1-1, SumElecs, abin-1000);
- Printf("(%.0f electrons and %.0f pions)",Histo1 -> GetEntries(), Histo2 -> GetEntries());
-
- // return the pion efficiency
- return SumPions;
+ // create TGraph to plot the pion efficiencies
+ TGraphErrors *gEffisLQ=0x0, *gEffisNN=0x0;
+ gEffisLQ = (TGraphErrors*)fContainer->At(kGraphLQ);
+ gEffisNN = (TGraphErrors*)fContainer->At(kGraphNN);
-}
-
-//________________________________________________________________________
-Double_t AliTRDpidChecker::GetError(Int_t Index1, Int_t Index2){
+ for(Int_t iBin = 0; iBin < AliTRDCalPID::kNMom; iBin++){
-
- const Int_t kBins = 12000; // binning of the histos and eficiency calculation
- const Float_t Multi = 0.9; // electron efficiency
- Int_t abinE, bbinE, cbinE;
- Double_t SumElecsE[kBins], SumPionsE[kBins]; // array of the integrated sum in each bin
- Double_t aBinSumE, bBinSumE; // content of a single bin
- Double_t EleEffi, PioEffi; // electron and pion efficiency
- Bool_t bFirst = 1; // checks if threshold is crossed for the first time
- Double_t fError = 0.; // error of the pion efficiency
-
-
- TH1F *Histo1 = (TH1F*)fObjectContainer->At(Index1); // electron histo
- TH1F *Histo2 = (TH1F*)fObjectContainer->At(Index2); // pion histo
-
- if(!(Histo1 -> GetEntries() > 0 && Histo2 -> GetEntries() > 0)){
- Printf("Warning: Histo momentum intervall %d has no Entries!", Index1-1);
- return -1.;
- }
+ Float_t momentum = AliTRDCalPID::GetMomentum(iBin);
+ gEffisLQ->SetPoint(iBin, momentum, PionEffiLQ[iBin]);
+ gEffisLQ->SetPointError(iBin, 0., PionEffiErrorLQ[iBin]);
- for(Int_t iBin = 0; iBin < kBins; iBin++){
- SumElecsE[iBin] = 0.;
- SumPionsE[iBin] = 0.;
+ gEffisNN->SetPoint(iBin, momentum, PionEffiNN[iBin]);
+ gEffisNN->SetPointError(iBin, 0., PionEffiErrorNN[iBin]);
}
- EleEffi = 0.;
- PioEffi = 0.;
- cbinE = -1;
+ fNRefFigures = 1;
+ return kTRUE; // testing protection
+}
- // calculate electron efficiency of each bin
- for (abinE = (Histo1 -> GetNbinsX())-2; abinE >= 0; abinE--){
- aBinSumE = 0;
- aBinSumE = Histo1 -> GetBinContent(abinE);
-
- SumElecsE[abinE] = SumElecsE[abinE+1] + aBinSumE;
- if((SumElecsE[abinE] >= Multi) && (bFirst == 1)){
- bFirst = 0;
- cbinE = abinE;
- EleEffi = (SumElecsE[cbinE]);
- }
- }
-
- // calculate pion efficiency of each bin
- for (bbinE = (Histo2 -> GetNbinsX())-2; bbinE >= abinE; bbinE--){
- bBinSumE = 0;
- bBinSumE = Histo2 -> GetBinContent(bbinE);
+//________________________________________________________________________
+void AliTRDpidChecker::Terminate(Option_t *)
+{
+ // Draw result to the screen
+ // Called once at the end of the query
- SumPionsE[bbinE] = SumPionsE[bbinE+1] + bBinSumE;
- if(bbinE == cbinE){
- PioEffi = (SumPionsE[cbinE]);
- }
+ fContainer = dynamic_cast<TObjArray*>(GetOutputData(0));
+ if (!fContainer) {
+ Printf("ERROR: list not available");
+ return;
}
-
-
- // pion efficiency vs electron efficiency
- TGraph *gEffis = new TGraph(kBins, SumElecsE, SumPionsE);
+}
- // use fit function to get derivate of the TGraph for error calculation
- TF1 *f1 = new TF1("f1","[0]*x*x+[1]*x+[2]", Multi-.05, Multi+.05);
- gEffis -> Fit("f1","Q","",Multi-.05, Multi+.05);
- Printf("Derivative at %4.2f : %f", Multi, f1 -> Derivative(Multi));
- // return the error of the pion efficiency
- fError = sqrt(((1/Histo2 -> GetEntries())*PioEffi*(1-PioEffi))+((f1 -> Derivative(Multi))*(f1 -> Derivative(Multi))*(1/Histo1 -> GetEntries())*EleEffi*(1-EleEffi)));
- return fError;
-}
git://git.uio.no / u / mrichter / AliRoot.git / blobdiff