[u/mrichter/AliRoot.git] / TRD / qaRec / AliTRDpidChecker.cxx

#include "TPDGCode.h"
#include "TF1.h"
#include "TH1F.h"
#include "TH2F.h"
#include "TGraph.h"
#include "TGraphErrors.h"

#include <TClonesArray.h>
#include <TObjArray.h>
#include <TList.h>

#include "AliPID.h"
#include "AliESDEvent.h"
#include "AliESDInputHandler.h"
#include "AliTrackReference.h"

#include "AliTRDrecoTask.h"
//#include "AliAnalysisManager.h"

#include "AliTRDtrackV1.h"
#include "AliTRDReconstructor.h"
#include "AliCDBManager.h"
#include "../Cal/AliTRDCalPID.h"


#include "AliTRDpidChecker.h"
#include "AliTRDtrackInfo/AliTRDtrackInfo.h"

// calculate pion efficiency at 90% electron efficiency for 11 momentum bins
// this task should be used with simulated data only

ClassImp(AliTRDpidChecker)

//________________________________________________________________________
AliTRDpidChecker::AliTRDpidChecker() 
  :AliTRDrecoTask("PID", "PID Checker")
  ,fReconstructor(0x0)
{
  //
  // Default constructor
  //

  fReconstructor = new AliTRDReconstructor();
  fReconstructor->SetRecoParam(AliTRDrecoParam::GetLowFluxParam());
}


//________________________________________________________________________
AliTRDpidChecker::~AliTRDpidChecker() 
{
  if(fReconstructor) delete fReconstructor;
}


//________________________________________________________________________
void AliTRDpidChecker::CreateOutputObjects()
{
  // Create histograms
  // Called once

  OpenFile(0, "RECREATE");
  fContainer = new TObjArray();
  fContainer->Add(new TH1F("hMom", "momentum distribution", AliTRDCalPID::kNMom, 0.5, 11.5));


  // 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++){
      fContainer->Add(new TH1F(Form("PID%d_%d_LQ",iPart,iMom), Form("PID distribution for %d_%d", iPart, iMom), kBins, -0.1, 1.1));
    }
  }

  // histos of the electron probability of all 5 particle species and 11 momenta for the neural network method
  const Float_t epsilon = 1.E-3; 
  for(Int_t iPart = 0; iPart < AliPID::kSPECIES; iPart++){
    for(Int_t iMom = 0; iMom < AliTRDCalPID::kNMom; iMom++){
      fContainer->Add(new TH1F(Form("PID%d_%d_NN",iPart,iMom), Form("PID distribution for %d_%d", iPart, iMom), kBins, 0.-epsilon, 1.+epsilon));
    }
  }

  // histos of the dE/dx distribution for all 5 particle species and 11 momenta 
  for(Int_t iPart = 0; iPart < AliPID::kSPECIES; iPart++){
    for(Int_t iMom = 0;  iMom < AliTRDCalPID::kNMom; iMom++){
      fContainer->Add(new TH1F(Form("dEdx%d_%d",iPart,iMom), Form("dEdx distribution for %d_%d", iPart, iMom), 200, 0, 10000));
    }
  }

  // frame and TGraph of the pion efficiencies
  //fContainer -> Add(new TH2F("hFrame", "", 10, 0.4, 12., 10, 0.0005, 0.1));
  fContainer -> Add(new TGraph(AliTRDCalPID::kNMom));
  fContainer -> Add(new TGraphErrors(AliTRDCalPID::kNMom));
  fContainer -> Add(new TGraph(AliTRDCalPID::kNMom));
  fContainer -> Add(new TGraphErrors(AliTRDCalPID::kNMom));
}

//________________________________________________________________________
void AliTRDpidChecker::Exec(Option_t *) 
{
  // Main loop
  // Called for each event


//   if(!AliTracker::GetFieldMap()){
//     AliMagFMaps* field = new AliMagFMaps("Maps","Maps", 2, 1., 10., AliMagFMaps::k5kG);
//     AliTracker::SetFieldMap(field, kTRUE);
//   }

  TH1F *hMom = (TH1F*)fContainer->UncheckedAt(0);	
  TH1F *hPIDLQ[AliPID::kSPECIES][AliTRDCalPID::kNMom];
  TH1F *hPIDNN[AliPID::kSPECIES][AliTRDCalPID::kNMom];
  TH1F *hdEdx[AliPID::kSPECIES][AliTRDCalPID::kNMom];

  for(Int_t iPart = 0; iPart < AliPID::kSPECIES; iPart++){
    for(Int_t iMom = 0; iMom < AliTRDCalPID::kNMom; iMom++){
      hPIDLQ[iPart][iMom] = (TH1F*)fContainer->At(iPart*AliTRDCalPID::kNMom+iMom+1);
      hPIDNN[iPart][iMom] = (TH1F*)fContainer->At(iPart*AliTRDCalPID::kNMom+iMom+1+AliPID::kSPECIES*AliTRDCalPID::kNMom);
      hdEdx[iPart][iMom]  = (TH1F*)fContainer->At(iPart*AliTRDCalPID::kNMom+iMom+1+AliPID::kSPECIES*AliTRDCalPID::kNMom*2);
    }
  }
	
  
  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;
  AliTrackReference     *ref = 0x0;
  AliExternalTrackParam *esd = 0x0;

  AliTRDseedV1 *TRDtracklet[6];
  for(Int_t iChamb = 0; iChamb < 6; iChamb++)
    TRDtracklet[iChamb] = 0x0;

  for(Int_t itrk=0; itrk<fTracks->GetEntriesFast(); itrk++){
    track = (AliTRDtrackInfo*)fTracks->UncheckedAt(itrk);
    if(!track->HasESDtrack()) continue;
    status = track->GetStatus();
    if(!(status&AliESDtrack::kTPCout)) continue;

    if(!(TRDtrack = track->GetTRDtrack())) continue; 
    //&&(track->GetNumberOfClustersRefit()

    // use only tracks that hit 6 chambers
    if(!(TRDtrack->GetNumberOfTracklets() == AliTRDCalPID::kNPlane)) continue;
     
    ref = track->GetTrackRef(0);
    esd = track->GetOuterParam();
    mom = ref ? ref->P(): esd->P();

    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;

    // set the reconstructor
    TRDtrack -> SetReconstructor(fReconstructor);

    
    // calculate the probabilities for electron probability using 2-dim LQ and the deposited charge per chamber and fill histograms
    fReconstructor -> SetOption("!nn");
    TRDtrack -> CookPID();
    Float_t SumdEdx[AliTRDCalPID::kNPlane];
    for(Int_t iChamb = 0; iChamb < AliTRDCalPID::kNPlane; 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));
      for(Int_t iChamb = 0; iChamb < AliTRDCalPID::kNPlane; iChamb++){
        hdEdx[AliPID::kElectron][iMomBin] -> Fill(SumdEdx[iChamb]);
      }
      break;
    case kMuonPlus:
    case kMuonMinus:
      hPIDLQ[AliPID::kMuon][iMomBin] -> Fill(TRDtrack -> GetPID(0));
      for(Int_t iChamb = 0; iChamb < AliTRDCalPID::kNPlane; iChamb++){
        hdEdx[AliPID::kMuon][iMomBin] -> Fill(SumdEdx[iChamb]);
      }
      break;
    case kPiPlus:
    case kPiMinus:
      hPIDLQ[AliPID::kPion][iMomBin] -> Fill(TRDtrack -> GetPID(0));
      for(Int_t iChamb = 0; iChamb < AliTRDCalPID::kNPlane; iChamb++){
        hdEdx[AliPID::kPion][iMomBin] -> Fill(SumdEdx[iChamb]);
      }
      break;
    case kKPlus:
    case kKMinus:
      hPIDLQ[AliPID::kKaon][iMomBin] -> Fill(TRDtrack -> GetPID(0));
      for(Int_t iChamb = 0; iChamb < AliTRDCalPID::kNPlane; iChamb++){
        hdEdx[AliPID::kKaon][iMomBin] -> Fill(SumdEdx[iChamb]);
      }
      break;
    case kProton:
    case kProtonBar:
      hPIDLQ[AliPID::kProton][iMomBin] -> Fill(TRDtrack -> GetPID(0));
      for(Int_t iChamb = 0; iChamb < AliTRDCalPID::kNPlane; iChamb++){
        hdEdx[AliPID::kProton][iMomBin] -> Fill(SumdEdx[iChamb]);
      }
      break;
    }


    // calculate the probabilities and fill histograms for electrons using NN
    fReconstructor -> SetOption("nn");
    TRDtrack->CookPID();
    switch(track->GetPDG()){
    case kElectron:
    case kPositron:
      hPIDNN[AliPID::kElectron][iMomBin] -> Fill(TRDtrack -> GetPID(0));
      break;
    case kMuonPlus:
    case kMuonMinus:
      hPIDNN[AliPID::kMuon][iMomBin] -> Fill(TRDtrack -> GetPID(0));
      break;
    case kPiPlus:
    case kPiMinus:
      hPIDNN[AliPID::kPion][iMomBin] -> Fill(TRDtrack -> GetPID(0));
      break;
    case kKPlus:
    case kKMinus:
      hPIDNN[AliPID::kKaon][iMomBin] -> Fill(TRDtrack -> GetPID(0));
      break;
    case kProton:
    case kProtonBar:
      hPIDNN[AliPID::kProton][iMomBin] -> Fill(TRDtrack -> GetPID(0));
      break;
    }
  }

  PostData(0, fContainer);
}

//________________________________________________________________________
void AliTRDpidChecker::Terminate(Option_t *) 
{
  // Draw result to the screen
  // Called once at the end of the query

  fContainer = dynamic_cast<TObjArray*>(GetOutputData(0));
  if (!fContainer) {
    Printf("ERROR: list not available");
    return;
  }

  
  // normalize the dE/dx histos
  const Int_t kNSpectra = AliPID::kSPECIES*AliTRDCalPID::kNMom; 
  TH1F *hdEdx[kNSpectra];
  for(Int_t iHist = 0; iHist < kNSpectra; iHist++){
    hdEdx[iHist] = (TH1F*)fContainer->At(111+iHist);
    if(hdEdx[iHist] -> GetEntries() > 0)
      hdEdx[iHist] -> Scale(1./hdEdx[iHist] -> GetEntries());
    else continue;
  }


  // container for the pion efficiencies and the errors
  Double_t  PionEffiLQ[AliTRDCalPID::kNMom], 
            PionEffiNN[AliTRDCalPID::kNMom],
            PionEffiErrorLQ[AliTRDCalPID::kNMom], 
            PionEffiErrorNN[AliTRDCalPID::kNMom];


  // 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;

  gEffisLQ = (TGraph*)fContainer->At(166);
  gEffisErrLQ = (TGraphErrors*)fContainer->At(167);
  gEffisNN = (TGraph*)fContainer->At(168);
  gEffisErrNN = (TGraphErrors*)fContainer->At(169);

  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]);

    gEffisNN->SetPoint(iBin, momentum, PionEffiNN[iBin]);
    gEffisErrNN->SetPoint(iBin, momentum, PionEffiNN[iBin]);
    gEffisErrNN->SetPointError(iBin, 0., PionEffiErrorNN[iBin]);
  }

  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");
}


//________________________________________________________________________
Double_t  AliTRDpidChecker::GetPionEfficiency(Int_t Index1, Int_t Index2){

  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*)fContainer->At(Index1);  // electron histo
  TH1F *Histo2 = (TH1F*)fContainer->At(Index2);  // pion histo


  SumElecs = 0.;
  if(!(Histo1 -> GetEntries() > 0 && Histo2 -> GetEntries() > 0)){
    Printf("AliTRDpidChecker::GetPionEfficiency : Warning: Histo momentum intervall %d has no Entries!", Index1-1);
    return -1.;
  }
  Histo1 -> Scale(1./Histo1->GetEntries());
  Histo2 -> Scale(1./Histo2->GetEntries());


  // 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;
    }
  }

    
  // calculate pion efficiency
  SumPions = 0.;
  for (bbin = (Histo2 -> GetNbinsX()); bbin >= abin; bbin--){	
    bBinSum = 0;
    bBinSum = Histo2 -> GetBinContent(bbin);
    if(!(bBinSum == 0)){
      SumPions = SumPions + bBinSum;
    }
  }
  
  
  // 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;

} 
  

//________________________________________________________________________
Double_t  AliTRDpidChecker::GetError(Int_t Index1, Int_t Index2){

  
  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*)fContainer->At(Index1);  // electron histo
  TH1F *Histo2 = (TH1F*)fContainer->At(Index2);  // pion histo

  if(Index1 > 10)              // print the correct momentum index for neural networks
    Index1 = Index1 - 55;
  if(!(Histo1 -> GetEntries() > 0 && Histo2 -> GetEntries() > 0)){
    Printf("Warning: Histo momentum intervall %d has no Entries!", Index1-1);
    return -1.;
  }

  for(Int_t iBin = 0; iBin < kBins; iBin++){
    SumElecsE[iBin] = 0.;
    SumPionsE[iBin] = 0.;
  }

  EleEffi = 0.;
  PioEffi = 0.;
  cbinE = -1;


  // 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);

    SumPionsE[bbinE] = SumPionsE[bbinE+1] + bBinSumE;
    if(bbinE == cbinE){
      PioEffi = (SumPionsE[cbinE]);
    }
  }
  

  // 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
  if(((1.-PioEffi) < 0) || ((1.-EleEffi) < 0)){
    Printf("AliTRDpidChecker::GetError : Warning: EleEffi or PioEffi > 1. Error can not be calculated. Please increase statistics or check your simulation!");
    return -1;
  }
  fError = sqrt(((1/Histo2 -> GetEntries())*PioEffi*(1-PioEffi))+((f1 -> Derivative(Multi))*(f1 -> Derivative(Multi))*(1/Histo1 -> GetEntries())*EleEffi*(1-EleEffi)));
  return fError;
}
Commit	Line	Data
773d3f8c	1	#include "TPDGCode.h"
	2	#include "TF1.h"
	3	#include "TH1F.h"
	4	#include "TH2F.h"
	5	#include "TGraph.h"
	6	#include "TGraphErrors.h"
	7
	8	#include <TClonesArray.h>
	9	#include <TObjArray.h>
	10	#include <TList.h>
	11
	12	#include "AliPID.h"
	13	#include "AliESDEvent.h"
	14	#include "AliESDInputHandler.h"
	15	#include "AliTrackReference.h"
	16
3d86166d	17	#include "AliTRDrecoTask.h"
3d86166d	18	//#include "AliAnalysisManager.h"
773d3f8c	19
	20	#include "AliTRDtrackV1.h"
	21	#include "AliTRDReconstructor.h"
	22	#include "AliCDBManager.h"
	23	#include "../Cal/AliTRDCalPID.h"
	24
	25
	26	#include "AliTRDpidChecker.h"
	27	#include "AliTRDtrackInfo/AliTRDtrackInfo.h"
	28
	29	// calculate pion efficiency at 90% electron efficiency for 11 momentum bins
	30	// this task should be used with simulated data only
	31
	32	ClassImp(AliTRDpidChecker)
	33
	34	//________________________________________________________________________
3d86166d	35	AliTRDpidChecker::AliTRDpidChecker()
3d86166d	36	:AliTRDrecoTask("PID", "PID Checker")
773d3f8c	37	,fReconstructor(0x0)
773d3f8c	38	{
	39	//
	40	// Default constructor
	41	//
	42
	43	fReconstructor = new AliTRDReconstructor();
	44	fReconstructor->SetRecoParam(AliTRDrecoParam::GetLowFluxParam());
773d3f8c	45	}
	46
	47
	48	//________________________________________________________________________
	49	AliTRDpidChecker::~AliTRDpidChecker()
	50	{
3d86166d	51	if(fReconstructor) delete fReconstructor;
773d3f8c	52	}
	53
	54
	55	//________________________________________________________________________
	56	void AliTRDpidChecker::CreateOutputObjects()
	57	{
	58	// Create histograms
	59	// Called once
	60
	61	OpenFile(0, "RECREATE");
3d86166d	62	fContainer = new TObjArray();
3d86166d	63	fContainer->Add(new TH1F("hMom", "momentum distribution", AliTRDCalPID::kNMom, 0.5, 11.5));
773d3f8c	64
	65
	66	// histos of the electron probability of all 5 particle species and 11 momenta for the 2-dim LQ method
	67	const Int_t kBins = 12000; // binning of the histos and eficiency calculation
	68	for(Int_t iPart = 0; iPart < AliPID::kSPECIES; iPart++){
	69	for(Int_t iMom = 0; iMom < AliTRDCalPID::kNMom; iMom++){
3d86166d	70	fContainer->Add(new TH1F(Form("PID%d_%d_LQ",iPart,iMom), Form("PID distribution for %d_%d", iPart, iMom), kBins, -0.1, 1.1));
773d3f8c	71	}
	72	}
	73
1a6d7c9a	74	// histos of the electron probability of all 5 particle species and 11 momenta for the neural network method
1a6d7c9a	75	const Float_t epsilon = 1.E-3;
773d3f8c	76	for(Int_t iPart = 0; iPart < AliPID::kSPECIES; iPart++){
773d3f8c	77	for(Int_t iMom = 0; iMom < AliTRDCalPID::kNMom; iMom++){
3d86166d	78	fContainer->Add(new TH1F(Form("PID%d_%d_NN",iPart,iMom), Form("PID distribution for %d_%d", iPart, iMom), kBins, 0.-epsilon, 1.+epsilon));
1a6d7c9a	79	}
	80	}
	81
	82	// histos of the dE/dx distribution for all 5 particle species and 11 momenta
	83	for(Int_t iPart = 0; iPart < AliPID::kSPECIES; iPart++){
	84	for(Int_t iMom = 0; iMom < AliTRDCalPID::kNMom; iMom++){
3d86166d	85	fContainer->Add(new TH1F(Form("dEdx%d_%d",iPart,iMom), Form("dEdx distribution for %d_%d", iPart, iMom), 200, 0, 10000));
773d3f8c	86	}
	87	}
	88
	89	// frame and TGraph of the pion efficiencies
3d86166d	90	//fContainer -> Add(new TH2F("hFrame", "", 10, 0.4, 12., 10, 0.0005, 0.1));
	91	fContainer -> Add(new TGraph(AliTRDCalPID::kNMom));
	92	fContainer -> Add(new TGraphErrors(AliTRDCalPID::kNMom));
	93	fContainer -> Add(new TGraph(AliTRDCalPID::kNMom));
	94	fContainer -> Add(new TGraphErrors(AliTRDCalPID::kNMom));
773d3f8c	95	}
	96
	97	//________________________________________________________________________
	98	void AliTRDpidChecker::Exec(Option_t *)
	99	{
	100	// Main loop
	101	// Called for each event
	102
	103
	104	// if(!AliTracker::GetFieldMap()){
	105	// AliMagFMaps* field = new AliMagFMaps("Maps","Maps", 2, 1., 10., AliMagFMaps::k5kG);
	106	// AliTracker::SetFieldMap(field, kTRUE);
	107	// }
	108
3d86166d	109	TH1F hMom = (TH1F)fContainer->UncheckedAt(0);
773d3f8c	110	TH1F *hPIDLQ[AliPID::kSPECIES][AliTRDCalPID::kNMom];
773d3f8c	111	TH1F *hPIDNN[AliPID::kSPECIES][AliTRDCalPID::kNMom];
1a6d7c9a	112	TH1F *hdEdx[AliPID::kSPECIES][AliTRDCalPID::kNMom];
773d3f8c	113
	114	for(Int_t iPart = 0; iPart < AliPID::kSPECIES; iPart++){
	115	for(Int_t iMom = 0; iMom < AliTRDCalPID::kNMom; iMom++){
3d86166d	116	hPIDLQ[iPart][iMom] = (TH1F)fContainer->At(iPartAliTRDCalPID::kNMom+iMom+1);
	117	hPIDNN[iPart][iMom] = (TH1F)fContainer->At(iPartAliTRDCalPID::kNMom+iMom+1+AliPID::kSPECIES*AliTRDCalPID::kNMom);
	118	hdEdx[iPart][iMom] = (TH1F)fContainer->At(iPartAliTRDCalPID::kNMom+iMom+1+AliPID::kSPECIESAliTRDCalPID::kNMom2);
773d3f8c	119	}
	120	}
	121
1a6d7c9a	122
773d3f8c	123	Int_t labelsacc[10000];
	124	memset(labelsacc, 0, sizeof(Int_t) * 10000);
	125
	126	Float_t mom;
	127	ULong_t status;
	128	Int_t nTRD = 0;
1a6d7c9a	129	Float_t *fdEdx;
773d3f8c	130
	131	AliTRDtrackInfo *track = 0x0;
	132	AliTRDtrackV1 *TRDtrack = 0x0;
	133	AliTrackReference *ref = 0x0;
	134	AliExternalTrackParam *esd = 0x0;
1a6d7c9a	135
	136	AliTRDseedV1 *TRDtracklet[6];
	137	for(Int_t iChamb = 0; iChamb < 6; iChamb++)
	138	TRDtracklet[iChamb] = 0x0;
	139
773d3f8c	140	for(Int_t itrk=0; itrk<fTracks->GetEntriesFast(); itrk++){
	141	track = (AliTRDtrackInfo*)fTracks->UncheckedAt(itrk);
	142	if(!track->HasESDtrack()) continue;
	143	status = track->GetStatus();
	144	if(!(status&AliESDtrack::kTPCout)) continue;
	145
	146	if(!(TRDtrack = track->GetTRDtrack())) continue;
	147	//&&(track->GetNumberOfClustersRefit()
1a6d7c9a	148
1a6d7c9a	149	// use only tracks that hit 6 chambers
773d3f8c	150	if(!(TRDtrack->GetNumberOfTracklets() == AliTRDCalPID::kNPlane)) continue;
	151
	152	ref = track->GetTrackRef(0);
	153	esd = track->GetOuterParam();
	154	mom = ref ? ref->P(): esd->P();
	155
	156	labelsacc[nTRD] = track->GetLabel();
	157	nTRD++;
	158
	159	// fill momentum histo to have the momentum distribution
	160	hMom -> Fill(mom);
	161
	162
	163	// set the 11 momentum bins
	164	Int_t iMomBin = -1;
	165	if(mom < .7) iMomBin = 0;
	166	else if(mom < .9) iMomBin = 1;
	167	else if(mom < 1.25) iMomBin = 2;
	168	else if(mom < 1.75) iMomBin = 3;
	169	else if(mom < 2.5) iMomBin = 4;
	170	else if(mom < 3.5) iMomBin = 5;
	171	else if(mom < 4.5) iMomBin = 6;
	172	else if(mom < 5.5) iMomBin = 7;
	173	else if(mom < 7.0) iMomBin = 8;
	174	else if(mom < 9.0) iMomBin = 9;
	175	else iMomBin = 10;
	176
	177	// set the reconstructor
	178	TRDtrack -> SetReconstructor(fReconstructor);
	179
	180
1a6d7c9a	181	// calculate the probabilities for electron probability using 2-dim LQ and the deposited charge per chamber and fill histograms
773d3f8c	182	fReconstructor -> SetOption("!nn");
773d3f8c	183	TRDtrack -> CookPID();
1a6d7c9a	184	Float_t SumdEdx[AliTRDCalPID::kNPlane];
	185	for(Int_t iChamb = 0; iChamb < AliTRDCalPID::kNPlane; iChamb++){
	186	TRDtracklet[iChamb] = TRDtrack -> GetTracklet(iChamb);
	187	SumdEdx[iChamb] = 0.;
	188	fdEdx = TRDtracklet[iChamb] -> GetdEdx();
	189	SumdEdx[iChamb] += fdEdx[0] + fdEdx[1] + fdEdx[2];
	190	}
	191
773d3f8c	192
	193	switch(track->GetPDG()){
	194	case kElectron:
	195	case kPositron:
	196	hPIDLQ[AliPID::kElectron][iMomBin] -> Fill(TRDtrack -> GetPID(0));
1a6d7c9a	197	for(Int_t iChamb = 0; iChamb < AliTRDCalPID::kNPlane; iChamb++){
	198	hdEdx[AliPID::kElectron][iMomBin] -> Fill(SumdEdx[iChamb]);
	199	}
773d3f8c	200	break;
	201	case kMuonPlus:
	202	case kMuonMinus:
	203	hPIDLQ[AliPID::kMuon][iMomBin] -> Fill(TRDtrack -> GetPID(0));
1a6d7c9a	204	for(Int_t iChamb = 0; iChamb < AliTRDCalPID::kNPlane; iChamb++){
	205	hdEdx[AliPID::kMuon][iMomBin] -> Fill(SumdEdx[iChamb]);
	206	}
773d3f8c	207	break;
	208	case kPiPlus:
	209	case kPiMinus:
	210	hPIDLQ[AliPID::kPion][iMomBin] -> Fill(TRDtrack -> GetPID(0));
1a6d7c9a	211	for(Int_t iChamb = 0; iChamb < AliTRDCalPID::kNPlane; iChamb++){
	212	hdEdx[AliPID::kPion][iMomBin] -> Fill(SumdEdx[iChamb]);
	213	}
773d3f8c	214	break;
	215	case kKPlus:
	216	case kKMinus:
	217	hPIDLQ[AliPID::kKaon][iMomBin] -> Fill(TRDtrack -> GetPID(0));
1a6d7c9a	218	for(Int_t iChamb = 0; iChamb < AliTRDCalPID::kNPlane; iChamb++){
	219	hdEdx[AliPID::kKaon][iMomBin] -> Fill(SumdEdx[iChamb]);
	220	}
773d3f8c	221	break;
	222	case kProton:
	223	case kProtonBar:
1a6d7c9a	224	hPIDLQ[AliPID::kProton][iMomBin] -> Fill(TRDtrack -> GetPID(0));
	225	for(Int_t iChamb = 0; iChamb < AliTRDCalPID::kNPlane; iChamb++){
	226	hdEdx[AliPID::kProton][iMomBin] -> Fill(SumdEdx[iChamb]);
	227	}
773d3f8c	228	break;
	229	}
	230
	231
	232	// calculate the probabilities and fill histograms for electrons using NN
	233	fReconstructor -> SetOption("nn");
	234	TRDtrack->CookPID();
	235	switch(track->GetPDG()){
	236	case kElectron:
	237	case kPositron:
	238	hPIDNN[AliPID::kElectron][iMomBin] -> Fill(TRDtrack -> GetPID(0));
	239	break;
	240	case kMuonPlus:
	241	case kMuonMinus:
	242	hPIDNN[AliPID::kMuon][iMomBin] -> Fill(TRDtrack -> GetPID(0));
	243	break;
	244	case kPiPlus:
	245	case kPiMinus:
	246	hPIDNN[AliPID::kPion][iMomBin] -> Fill(TRDtrack -> GetPID(0));
	247	break;
	248	case kKPlus:
	249	case kKMinus:
	250	hPIDNN[AliPID::kKaon][iMomBin] -> Fill(TRDtrack -> GetPID(0));
	251	break;
	252	case kProton:
	253	case kProtonBar:
1a6d7c9a	254	hPIDNN[AliPID::kProton][iMomBin] -> Fill(TRDtrack -> GetPID(0));
773d3f8c	255	break;
	256	}
	257	}
	258
3d86166d	259	PostData(0, fContainer);
773d3f8c	260	}
	261
	262	//________________________________________________________________________
	263	void AliTRDpidChecker::Terminate(Option_t *)
	264	{
	265	// Draw result to the screen
	266	// Called once at the end of the query
	267
3d86166d	268	fContainer = dynamic_cast<TObjArray*>(GetOutputData(0));
3d86166d	269	if (!fContainer) {
773d3f8c	270	Printf("ERROR: list not available");
	271	return;
	272	}
	273
	274
1a6d7c9a	275	// normalize the dE/dx histos
	276	const Int_t kNSpectra = AliPID::kSPECIES*AliTRDCalPID::kNMom;
	277	TH1F *hdEdx[kNSpectra];
	278	for(Int_t iHist = 0; iHist < kNSpectra; iHist++){
3d86166d	279	hdEdx[iHist] = (TH1F*)fContainer->At(111+iHist);
1a6d7c9a	280	if(hdEdx[iHist] -> GetEntries() > 0)
	281	hdEdx[iHist] -> Scale(1./hdEdx[iHist] -> GetEntries());
	282	else continue;
	283	}
	284
	285
773d3f8c	286	// container for the pion efficiencies and the errors
	287	Double_t PionEffiLQ[AliTRDCalPID::kNMom],
	288	PionEffiNN[AliTRDCalPID::kNMom],
	289	PionEffiErrorLQ[AliTRDCalPID::kNMom],
	290	PionEffiErrorNN[AliTRDCalPID::kNMom];
	291
	292
	293	// calculate the pion efficiencies and the errors for 90% electron efficiency (2-dim LQ)
	294	for(Int_t iMom = 0; iMom < AliTRDCalPID::kNMom; iMom++){
	295	PionEffiLQ[iMom] = GetPionEfficiency(iMom+1,iMom+23);
	296	PionEffiErrorLQ[iMom] = GetError(iMom+1,iMom+23);
	297	Printf("Pion Efficiency for 2-dim LQ is : %f +/- %f\n\n", PionEffiLQ[iMom], PionEffiErrorLQ[iMom]);
	298	}
	299
	300	// calculate the pion efficiencies and the errors for 90% electron efficiency (NN)
	301	for(Int_t iMom = 0; iMom < AliTRDCalPID::kNMom; iMom++){
	302	PionEffiNN[iMom] = GetPionEfficiency(iMom+56,iMom+78);
	303	PionEffiErrorNN[iMom] = GetError(iMom+56,iMom+78);
	304	Printf("Pion Efficiency for NN is : %f +/- %f\n\n", PionEffiNN[iMom], PionEffiErrorNN[iMom]);
	305	}
	306
	307
	308	// create TGraph to plot the pion efficiencies
	309	TGraph gEffisLQ=0x0, gEffisNN=0x0;
	310	TGraphErrors gEffisErrLQ=0x0, gEffisErrNN=0x0;
	311
3d86166d	312	gEffisLQ = (TGraph*)fContainer->At(166);
	313	gEffisErrLQ = (TGraphErrors*)fContainer->At(167);
	314	gEffisNN = (TGraph*)fContainer->At(168);
	315	gEffisErrNN = (TGraphErrors*)fContainer->At(169);
773d3f8c	316
	317	for(Int_t iBin = 0; iBin < AliTRDCalPID::kNMom; iBin++){
	318	Float_t momentum = AliTRDCalPID::GetMomentum(iBin);
	319	gEffisLQ->SetPoint(iBin, momentum, PionEffiLQ[iBin]);
	320	gEffisErrLQ->SetPoint(iBin, momentum, PionEffiLQ[iBin]);
	321	gEffisErrLQ->SetPointError(iBin, 0., PionEffiErrorLQ[iBin]);
	322
	323	gEffisNN->SetPoint(iBin, momentum, PionEffiNN[iBin]);
	324	gEffisErrNN->SetPoint(iBin, momentum, PionEffiNN[iBin]);
	325	gEffisErrNN->SetPointError(iBin, 0., PionEffiErrorNN[iBin]);
	326	}
	327
	328	gEffisLQ -> SetNameTitle("gEffisLQ", "Efficiencies of the 2-dim LQ method");
	329	gEffisErrLQ -> SetNameTitle("gEffisLQErr", "Efficiencies and Errors of the 2-dim LQ method");
	330	gEffisNN -> SetNameTitle("gEffisNN", "Efficiencies of the NN method");
	331	gEffisErrNN -> SetNameTitle("gEffisNNErr", "Efficiencies and Errors of the NN method");
	332	}
	333
	334
	335	//________________________________________________________________________
	336	Double_t AliTRDpidChecker::GetPionEfficiency(Int_t Index1, Int_t Index2){
	337
	338	Float_t Multi = 0.9; // electron efficiency
	339	Int_t abin, bbin;
	340	Double_t SumElecs, SumPions; // integrated sum of elecs and pions in the histos
	341	Double_t aBinSum, bBinSum; // content of a single bin in the histos
	342
3d86166d	343	TH1F Histo1 = (TH1F)fContainer->At(Index1); // electron histo
3d86166d	344	TH1F Histo2 = (TH1F)fContainer->At(Index2); // pion histo
773d3f8c	345
	346
	347	SumElecs = 0.;
	348	if(!(Histo1 -> GetEntries() > 0 && Histo2 -> GetEntries() > 0)){
1a6d7c9a	349	Printf("AliTRDpidChecker::GetPionEfficiency : Warning: Histo momentum intervall %d has no Entries!", Index1-1);
773d3f8c	350	return -1.;
	351	}
	352	Histo1 -> Scale(1./Histo1->GetEntries());
	353	Histo2 -> Scale(1./Histo2->GetEntries());
	354
	355
	356	// calculate threshold for pion efficiency
	357	for (abin = (Histo1 -> GetNbinsX()); abin >= 0; abin--){
	358	aBinSum = 0;
	359	aBinSum = Histo1 -> GetBinContent(abin);
	360	if(!(aBinSum == 0)){
	361	SumElecs = SumElecs + aBinSum;
	362	}
	363	if (SumElecs >= Multi){
	364	break;
	365	}
	366	}
	367
	368
	369	// calculate pion efficiency
	370	SumPions = 0.;
	371	for (bbin = (Histo2 -> GetNbinsX()); bbin >= abin; bbin--){
	372	bBinSum = 0;
	373	bBinSum = Histo2 -> GetBinContent(bbin);
	374	if(!(bBinSum == 0)){
	375	SumPions = SumPions + bBinSum;
	376	}
	377	}
	378
	379
	380	// print the electron efficiency and its cuts
	381	Printf("Cut for momentum intervall %d and electron efficiency of %f for: 0.%d", Index1-1, SumElecs, abin-1000);
	382	Printf("(%.0f electrons and %.0f pions)",Histo1 -> GetEntries(), Histo2 -> GetEntries());
	383
	384
	385	// return the pion efficiency
	386	return SumPions;
	387
	388	}
	389
	390
	391	//________________________________________________________________________
	392	Double_t AliTRDpidChecker::GetError(Int_t Index1, Int_t Index2){
	393
	394
	395	const Int_t kBins = 12000; // binning of the histos and eficiency calculation
	396	const Float_t Multi = 0.9; // electron efficiency
	397	Int_t abinE, bbinE, cbinE;
	398	Double_t SumElecsE[kBins], SumPionsE[kBins]; // array of the integrated sum in each bin
	399	Double_t aBinSumE, bBinSumE; // content of a single bin
	400	Double_t EleEffi, PioEffi; // electron and pion efficiency
	401	Bool_t bFirst = 1; // checks if threshold is crossed for the first time
	402	Double_t fError = 0.; // error of the pion efficiency
	403
	404
3d86166d	405	TH1F Histo1 = (TH1F)fContainer->At(Index1); // electron histo
3d86166d	406	TH1F Histo2 = (TH1F)fContainer->At(Index2); // pion histo
773d3f8c	407
1a6d7c9a	408	if(Index1 > 10) // print the correct momentum index for neural networks
1a6d7c9a	409	Index1 = Index1 - 55;
773d3f8c	410	if(!(Histo1 -> GetEntries() > 0 && Histo2 -> GetEntries() > 0)){
	411	Printf("Warning: Histo momentum intervall %d has no Entries!", Index1-1);
	412	return -1.;
	413	}
	414
	415	for(Int_t iBin = 0; iBin < kBins; iBin++){
	416	SumElecsE[iBin] = 0.;
	417	SumPionsE[iBin] = 0.;
	418	}
	419
	420	EleEffi = 0.;
	421	PioEffi = 0.;
	422	cbinE = -1;
	423
	424
	425	// calculate electron efficiency of each bin
	426	for (abinE = (Histo1 -> GetNbinsX())-2; abinE >= 0; abinE--){
	427	aBinSumE = 0;
	428	aBinSumE = Histo1 -> GetBinContent(abinE);
	429
	430	SumElecsE[abinE] = SumElecsE[abinE+1] + aBinSumE;
	431	if((SumElecsE[abinE] >= Multi) && (bFirst == 1)){
	432	bFirst = 0;
	433	cbinE = abinE;
	434	EleEffi = (SumElecsE[cbinE]);
	435	}
	436	}
	437
	438
	439	// calculate pion efficiency of each bin
	440	for (bbinE = (Histo2 -> GetNbinsX())-2; bbinE >= abinE; bbinE--){
	441	bBinSumE = 0;
	442	bBinSumE = Histo2 -> GetBinContent(bbinE);
	443
	444	SumPionsE[bbinE] = SumPionsE[bbinE+1] + bBinSumE;
	445	if(bbinE == cbinE){
	446	PioEffi = (SumPionsE[cbinE]);
	447	}
	448	}
	449
	450
	451	// pion efficiency vs electron efficiency
	452	TGraph *gEffis = new TGraph(kBins, SumElecsE, SumPionsE);
	453
	454	// use fit function to get derivate of the TGraph for error calculation
	455	TF1 f1 = new TF1("f1","[0]xx+[1]x+[2]", Multi-.05, Multi+.05);
	456	gEffis -> Fit("f1","Q","",Multi-.05, Multi+.05);
	457	Printf("Derivative at %4.2f : %f", Multi, f1 -> Derivative(Multi));
	458
	459	// return the error of the pion efficiency
1a6d7c9a	460	if(((1.-PioEffi) < 0) \|\| ((1.-EleEffi) < 0)){
	461	Printf("AliTRDpidChecker::GetError : Warning: EleEffi or PioEffi > 1. Error can not be calculated. Please increase statistics or check your simulation!");
	462	return -1;
	463	}
773d3f8c	464	fError = sqrt(((1/Histo2 -> GetEntries())PioEffi(1-PioEffi))+((f1 -> Derivative(Multi))(f1 -> Derivative(Multi))(1/Histo1 -> GetEntries())EleEffi(1-EleEffi)));
	465	return fError;
	466	}