[u/mrichter/AliRoot.git] / PWG3 / hfe / AliHFEtrdPIDqaV1.cxx

/**************************************************************************
* 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.                  *
**************************************************************************/

/* $Id$ */

//
// Class AliHFEtrdPIDqaV1
// Monitoring TRD PID in the HFE PID montioring framework. The following
// quantities are monitored:
//   TRD electron likelihood
//   TRD dE/dx (Absolute values)
//   TPC dE/dx (Number of sigmas, control histogram)
// (Always as function of momentum, particle species and centrality 
// before and after cut)
// More information about the PID monitoring framework can be found in
// AliHFEpidQAmanager.cxx and AliHFEdetPIDqa.cxx
//
// Author:
//    Markus Fasel <M.Fasel@gsi.de>
//

#include <TAxis.h>
#include <TBrowser.h>
#include <TH2.h>
#include <THnSparse.h>
#include <TString.h>

#include "AliESDtrack.h"
#include "AliLog.h"
#include "AliPID.h"

#include "AliHFEcollection.h"
#include "AliHFEpidBase.h"
#include "AliHFEpidQAmanager.h"
#include "AliHFEpidTPC.h"
#include "AliHFEpidTRD.h"
#include "AliHFEtrdPIDqaV1.h"

ClassImp(AliHFEtrdPIDqaV1)

//____________________________________________________________
AliHFEtrdPIDqaV1::AliHFEtrdPIDqaV1():
    AliHFEdetPIDqa(),
    fHistos(NULL)
{
  //
  // Dummy constructor
  //
}

//____________________________________________________________
AliHFEtrdPIDqaV1::AliHFEtrdPIDqaV1(const Char_t *name):
    AliHFEdetPIDqa(name, "QA for TRD"),
    fHistos(NULL)
{
  //
  // Default constructor
  //
}

//____________________________________________________________
AliHFEtrdPIDqaV1::AliHFEtrdPIDqaV1(const AliHFEtrdPIDqaV1 &o):
    AliHFEdetPIDqa(o),
    fHistos(NULL)
{
  //
  // Copy constructor
  //
}

//____________________________________________________________
AliHFEtrdPIDqaV1 &AliHFEtrdPIDqaV1::operator=(const AliHFEtrdPIDqaV1 &o){
  //
  // Make assignment
  //
  AliHFEdetPIDqa::operator=(o);
  fHistos = o.fHistos;
  
  return *this;
}

//_________________________________________________________
Long64_t AliHFEtrdPIDqaV1::Merge(TCollection *coll){
  //
  // Merge with other objects
  //
  if(!coll) return 0;
  if(coll->IsEmpty()) return 1;

  TIter it(coll);
  AliHFEtrdPIDqaV1 *refQA = NULL;
  TObject *o = NULL;
  Long64_t count = 0;
  TList listHistos;
  while((o = it())){
    refQA = dynamic_cast<AliHFEtrdPIDqaV1 *>(o);
    if(!refQA) continue;

    listHistos.Add(refQA->fHistos);
    count++; 
  }
  fHistos->Merge(&listHistos);
  return count + 1;
}

//_________________________________________________________
void AliHFEtrdPIDqaV1::Browse(TBrowser *b){
  //
  // Browse the PID QA
  //
  if(b){
    if(fHistos){
      b->Add(fHistos, fHistos->GetName());

      // Make Projections of the dE/dx Spectra and add them to a new Folder
      TString specnames[4] = {"All", "Electrons", "Pions", "Protons"};
      Int_t specind[4] = {-1, AliPID::kElectron, AliPID::kPion, AliPID::kProton};
      TList *listTM = new TList;
      listTM->SetOwner();
      TList *listLike = new TList;
      listLike->SetOwner();
      TList *listCharge = new TList;
      listCharge->SetOwner();
      TList *listTPCnsigma = new TList;
      listTPCnsigma->SetOwner();

      TH2 *hptr = NULL; 
      for(Int_t ispec = 0; ispec < 4; ispec++){
        for(Int_t istep = 0; istep < 2; istep++){
          hptr = MakeTRDspectrumTM(static_cast<AliHFEdetPIDqa::EStep_t>(istep), specind[ispec]);
          if(hptr){
            hptr->SetName(Form("hTRDtm%s%s", specnames[ispec].Data(), istep == 0 ? "Before" : "After"));
            listTM->Add(hptr);
          }
          hptr = MakeTRDlikelihoodDistribution(static_cast<AliHFEdetPIDqa::EStep_t>(istep), specind[ispec]);
          hptr->SetName(Form("hTRDlike%s%s", specnames[ispec].Data(), istep == 0 ? "Before" : "After"));
          listLike->Add(hptr);
          hptr = MakeTRDchargeDistribution(static_cast<AliHFEdetPIDqa::EStep_t>(istep), specind[ispec]);
          hptr->SetName(Form("hTRDcharge%s%s", specnames[ispec].Data(), istep == 0 ? "Before" : "After"));
          listCharge->Add(hptr);
          hptr = MakeTPCspectrumNsigma(static_cast<AliHFEdetPIDqa::EStep_t>(istep), specind[ispec]);
          hptr->SetName(Form("hTPCspectrum%s%s", specnames[ispec].Data(), istep == 0 ? "Before" : "After"));
          listTPCnsigma->Add(hptr);
        }
      }
      
      b->Add(listTM, "Projections Truncated Mean");
      b->Add(listLike, "Projections Likelihood distribution");
      b->Add(listCharge, "Projections Tracklet Charge");
      b->Add(listTPCnsigma, "Projections TPC spectra");
    }
  }
}

//____________________________________________________________
void AliHFEtrdPIDqaV1::Initialize(){
  //
  // Initialize QA histos for TRD PID
  //

  AliDebug(1, "Initializing PID QA for TRD");
  // Make common binning
  const Int_t kPIDbins = AliPID::kSPECIES + 1;
  const Int_t kPbins = 100;
  const Int_t kSteps = 2;
  const Double_t kMinPID = -1;
  const Double_t kMinP = 0.;
  const Double_t kMaxPID = (Double_t)AliPID::kSPECIES;
  const Double_t kMaxP = 20.;

  fHistos = new AliHFEcollection("trdqahistos", "Collection of TRD QA histograms");
  
  // Create Control Histogram monitoring the TPC sigma between the selection steps
  const Int_t kTPCSigmaBins = 140;
  Int_t nBinsTPCSigma[4] = {kPIDbins, kPbins, kTPCSigmaBins, kSteps};
  Double_t minTPCSigma[4] = {kMinPID, kMinP, -12., 0};
  Double_t maxTPCSigma[4] = {kMaxPID, kMaxP, 12., 2.};
  fHistos->CreateTHnSparse("hTPCsigma", "TPC sigma; species p [GeV/c]; TPC dEdx - <dE/dx>|_{el} [#sigma]; selection step", 4, nBinsTPCSigma, minTPCSigma, maxTPCSigma);
  fHistos->Sumw2("hTPCsigma");
  // Create Monitoring histogram for the Likelihood distribution
  const Int_t kTRDLikelihoodBins = 100;
  Int_t nBinsTRDlike[4] = {kPIDbins, kPbins, kTRDLikelihoodBins, kSteps};
  Double_t minTRDlike[4] = {kMinPID, kMinP, 0., 0.};
  Double_t maxTRDlike[4] = {kMaxPID, kMaxP, 1., 2.};
  fHistos->CreateTHnSparse("hTRDlikelihood", "TRD Likelihood Distribution; species; p [GeV/c]; TRD electron Likelihood; selection step", 4, nBinsTRDlike, minTRDlike, maxTRDlike);
  fHistos->Sumw2("hTRDlikelihood");
  // Create Monitoring histogram for the TRD total charge
  const Int_t kTRDchargeBins = 1000;
  Int_t nBinsTRDcharge[4] = {kPIDbins, kPbins, kTRDchargeBins, kSteps};
  Double_t minTRDcharge[4] = {kMinPID, kMinP, 0., 0.};
  Double_t maxTRDcharge[4] = {kMaxPID, kMaxP, 100000., 2.};
  fHistos->CreateTHnSparse("hTRDcharge", "Total TRD charge; species; p [GeV/c]; TRD charge [a.u.]; selection step", 4, nBinsTRDcharge, minTRDcharge, maxTRDcharge);
  fHistos->Sumw2("hTRDcharge");
  // Monitoring of the TRD truncated mean according to version 1
  const Int_t kTRDtmBins = 1000;
  Int_t nBinsTRDtm[4] = {kPIDbins, kPbins, kTRDtmBins, kSteps};
  Double_t minTRDtm[4] = {kMinPID, kMinP, 0., 0.};
  Double_t maxTRDtm[4] = {kMaxPID, kMaxP, 20000., 2.};
  fHistos->CreateTHnSparse("hTRDtruncatedMean", "TRD truncated Mean; species; p [GeV/c]; TRD signal [a.u.]; selection step", 4, nBinsTRDtm, minTRDtm, maxTRDtm);
  fHistos->Sumw2("hTRDtruncatedMean");
}

//____________________________________________________________
void AliHFEtrdPIDqaV1::ProcessTrack(const AliHFEpidObject *track, AliHFEdetPIDqa::EStep_t step){
  //
  // Process the track, fill the containers 
  //
  AliDebug(1, Form("QA started for TRD PID for step %d", (Int_t)step));
  Int_t species = track->GetAbInitioPID();
  if(species >= AliPID::kSPECIES) species = -1;
  AliHFEpidObject::AnalysisType_t anatype = track->IsESDanalysis() ? AliHFEpidObject::kESDanalysis : AliHFEpidObject::kAODanalysis;

  AliHFEpidTRD *trdpid = dynamic_cast<AliHFEpidTRD *>(fQAmanager->GetDetectorPID(AliHFEpid::kTRDpid));
  AliHFEpidTPC *tpcpid = dynamic_cast<AliHFEpidTPC *>(fQAmanager->GetDetectorPID(AliHFEpid::kTPCpid));
 
  Double_t container[4];
  container[0] = species;
  container[1] = trdpid ? trdpid->GetP(track->GetRecTrack(), anatype) : 0.;
  container[2] = tpcpid ? tpcpid->NumberOfSigmas(track->GetRecTrack(), AliPID::kElectron, anatype) : 0.;
  container[3] = step;
  fHistos->Fill("hTPCsigma", container);

  container[2] = trdpid->GetElectronLikelihood(track->GetRecTrack(), anatype);
  fHistos->Fill("hTRDlikelihood", container);

  if(track->IsESDanalysis()){
    container[2] = trdpid->GetTRDSignalV1(dynamic_cast<const AliESDtrack *>(track->GetRecTrack()));
    fHistos->Fill("hTRDtruncatedMean", container);
  }
  for(UInt_t ily = 0; ily < 6; ily++){
    container[2] = trdpid->GetChargeLayer(track->GetRecTrack(), ily, anatype);
    fHistos->Fill("hTRDcharge", container);
  }
}

//_________________________________________________________
TH2 *AliHFEtrdPIDqaV1::MakeTPCspectrumNsigma(AliHFEdetPIDqa::EStep_t step, Int_t species){
  //
  // Get the TPC control histogram for the TRD selection step (either before or after PID)
  //
  THnSparseF *histo = dynamic_cast<THnSparseF *>(fHistos->Get("hTPCsigma"));
  if(!histo){
    AliError("QA histogram monitoring TPC nSigma not available");
    return NULL;
  }
  if(species > -1 && species < AliPID::kSPECIES){
    // cut on species (if available)
    histo->GetAxis(0)->SetRange(species + 2, species + 2); // undef + underflow
  }
  histo->GetAxis(3)->SetRange(step + 1, step + 1); 

  TH2 *hSpec = histo->Projection(2, 1);
  // construct title and name
  TString stepname = step == AliHFEdetPIDqa::kBeforePID ? "before" : "after";
  TString speciesname = species > -1 && species < AliPID::kSPECIES ? AliPID::ParticleName(species) : "all Particles";
  TString specID = species > -1 && species < AliPID::kSPECIES ? AliPID::ParticleName(species) : "unid";
  TString histname = Form("hSigmaTPC%s%s", specID.Data(), stepname.Data());
  TString histtitle = Form("TPC Sigma for %s %s PID", speciesname.Data(), stepname.Data());
  hSpec->SetName(histname.Data());
  hSpec->SetTitle(histtitle.Data());

  // Unset range on the original histogram
  histo->GetAxis(0)->SetRange(0, histo->GetAxis(0)->GetNbins());
  histo->GetAxis(2)->SetRange(0, histo->GetAxis(2)->GetNbins());
  return hSpec; 
}

//_________________________________________________________
TH2 *AliHFEtrdPIDqaV1::MakeTRDspectrumTM(AliHFEdetPIDqa::EStep_t step, Int_t species){
  //
  // Get the TPC control histogram for the TRD selection step (either before or after PID)
  //
  THnSparseF *histo = dynamic_cast<THnSparseF *>(fHistos->Get("hTRDtruncatedMean"));
  if(!histo){
    AliError("QA histogram monitoring TPC nSigma not available");
    return NULL;
  }
  if(species > -1 && species < AliPID::kSPECIES){
    // cut on species (if available)
    histo->GetAxis(0)->SetRange(species + 2, species + 2); // undef + underflow
  }
  histo->GetAxis(3)->SetRange(step + 1, step + 1); 

  TH2 *hSpec = histo->Projection(2, 1);
  // construct title and name
  TString stepname = step == AliHFEdetPIDqa::kBeforePID ? "before" : "after";
  TString speciesname = species > -1 && species < AliPID::kSPECIES ? AliPID::ParticleName(species) : "all Particles";
  TString specID = species > -1 && species < AliPID::kSPECIES ? AliPID::ParticleName(species) : "unid";
  TString histname = Form("hTMTRD%s%s", specID.Data(), stepname.Data());
  TString histtitle = Form("TRD Truncated Mean for %s %s PID", speciesname.Data(), stepname.Data());
  hSpec->SetName(histname.Data());
  hSpec->SetTitle(histtitle.Data());

  // Unset range on the original histogram
  histo->GetAxis(0)->SetRange(0, histo->GetAxis(0)->GetNbins());
  histo->GetAxis(2)->SetRange(0, histo->GetAxis(2)->GetNbins());
  return hSpec; 
}

//_________________________________________________________
TH2 *AliHFEtrdPIDqaV1::MakeTRDlikelihoodDistribution(AliHFEdetPIDqa::EStep_t step, Int_t species){
  //
  // Make Histogram for TRD Likelihood distribution
  //
  THnSparseF *histo = dynamic_cast<THnSparseF *>(fHistos->Get("hTRDlikelihood"));
  if(!histo){
    AliError("QA histogram monitoring TRD Electron Likelihood not available");
    return NULL;
  }
  if(species > -1 && species < AliPID::kSPECIES){
    // cut on species (if available)
    histo->GetAxis(0)->SetRange(species + 2, species + 2); // undef + underflow
  }
  histo->GetAxis(3)->SetRangeUser(step + 1, step + 1);

  TH2 *hSpec = histo->Projection(2, 1);
  // construct title and name
  TString stepname = step == AliHFEdetPIDqa::kBeforePID ? "before" : "after";
  TString speciesname = species > -1 && species < AliPID::kSPECIES ? AliPID::ParticleName(species) : "all Particles";
  TString specID = species > -1 && species < AliPID::kSPECIES ? AliPID::ParticleName(species) : "unid";
  TString histname = Form("hLikeElTRD%s%s", specID.Data(), stepname.Data());
  TString histtitle = Form("TRD electron Likelihood for %s %s PID", speciesname.Data(), stepname.Data());
  hSpec->SetName(histname.Data());
  hSpec->SetTitle(histtitle.Data());

  // Unset range on the original histogram
  histo->GetAxis(0)->SetRange(0, histo->GetAxis(0)->GetNbins());
  histo->GetAxis(2)->SetRange(0, histo->GetAxis(2)->GetNbins());
  return hSpec; 
}

//_________________________________________________________
TH2 *AliHFEtrdPIDqaV1::MakeTRDchargeDistribution(AliHFEdetPIDqa::EStep_t step, Int_t species){
  //
  // Make Histogram for TRD Likelihood distribution
  //
  THnSparseF *histo = dynamic_cast<THnSparseF *>(fHistos->Get("hTRDcharge"));
  if(!histo){
    AliError("QA histogram monitoring TRD total charge not available");
    return NULL;
  }
  if(species > -1 && species < AliPID::kSPECIES){
    // cut on species (if available)
    histo->GetAxis(0)->SetRange(species + 2, species + 2); // undef + underflow
  }
  histo->GetAxis(3)->SetRange(step + 1, step + 1);

  TH2 *hSpec = histo->Projection(2, 1);
  // construct title and name
  TString stepname = step == AliHFEdetPIDqa::kBeforePID ? "before" : "after";
  TString speciesname = species > -1 && species < AliPID::kSPECIES ? AliPID::ParticleName(species) : "all Particles";
  TString specID = species > -1 && species < AliPID::kSPECIES ? AliPID::ParticleName(species) : "unid";
  TString histname = Form("hChargeTRD%s%s", specID.Data(), stepname.Data());
  TString histtitle = Form("TRD total charge for %s %s PID", speciesname.Data(), stepname.Data());
  hSpec->SetName(histname.Data());
  hSpec->SetTitle(histtitle.Data());

  // Unset range on the original histogram
  histo->GetAxis(0)->SetRange(0, histo->GetAxis(0)->GetNbins());
  histo->GetAxis(2)->SetRange(0, histo->GetAxis(2)->GetNbins());
  return hSpec; 
}
Commit	Line	Data
3a72645a	1	/**************************************************************************
	2	* Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
	3	* *
	4	* Author: The ALICE Off-line Project. *
	5	* Contributors are mentioned in the code where appropriate. *
	6	* *
	7	* Permission to use, copy, modify and distribute this software and its *
	8	* documentation strictly for non-commercial purposes is hereby granted *
	9	* without fee, provided that the above copyright notice appears in all *
	10	* copies and that both the copyright notice and this permission notice *
	11	* appear in the supporting documentation. The authors make no claims *
	12	* about the suitability of this software for any purpose. It is *
	13	* provided "as is" without express or implied warranty. *
	14	**************************************************************************/
27de2dfb	15
	16	/* $Id$ */
	17
3a72645a	18	//
	19	// Class AliHFEtrdPIDqaV1
	20	// Monitoring TRD PID in the HFE PID montioring framework. The following
	21	// quantities are monitored:
	22	// TRD electron likelihood
	23	// TRD dE/dx (Absolute values)
	24	// TPC dE/dx (Number of sigmas, control histogram)
	25	// (Always as function of momentum, particle species and centrality
	26	// before and after cut)
	27	// More information about the PID monitoring framework can be found in
	28	// AliHFEpidQAmanager.cxx and AliHFEdetPIDqa.cxx
	29	//
	30	// Author:
	31	// Markus Fasel <M.Fasel@gsi.de>
	32	//
	33
	34	#include <TAxis.h>
e3ae862b	35	#include <TBrowser.h>
3a72645a	36	#include <TH2.h>
	37	#include <THnSparse.h>
	38	#include <TString.h>
	39
e3ae862b	40	#include "AliESDtrack.h"
6555e2ad	41	#include "AliLog.h"
3a72645a	42	#include "AliPID.h"
	43
	44	#include "AliHFEcollection.h"
	45	#include "AliHFEpidBase.h"
6555e2ad	46	#include "AliHFEpidQAmanager.h"
	47	#include "AliHFEpidTPC.h"
	48	#include "AliHFEpidTRD.h"
3a72645a	49	#include "AliHFEtrdPIDqaV1.h"
	50
	51	ClassImp(AliHFEtrdPIDqaV1)
	52
	53	//____________________________________________________________
	54	AliHFEtrdPIDqaV1::AliHFEtrdPIDqaV1():
	55	AliHFEdetPIDqa(),
	56	fHistos(NULL)
	57	{
	58	//
	59	// Dummy constructor
	60	//
	61	}
	62
	63	//____________________________________________________________
	64	AliHFEtrdPIDqaV1::AliHFEtrdPIDqaV1(const Char_t *name):
	65	AliHFEdetPIDqa(name, "QA for TRD"),
	66	fHistos(NULL)
	67	{
	68	//
	69	// Default constructor
	70	//
	71	}
	72
	73	//____________________________________________________________
	74	AliHFEtrdPIDqaV1::AliHFEtrdPIDqaV1(const AliHFEtrdPIDqaV1 &o):
	75	AliHFEdetPIDqa(o),
	76	fHistos(NULL)
	77	{
	78	//
	79	// Copy constructor
	80	//
	81	}
	82
	83	//____________________________________________________________
	84	AliHFEtrdPIDqaV1 &AliHFEtrdPIDqaV1::operator=(const AliHFEtrdPIDqaV1 &o){
	85	//
	86	// Make assignment
	87	//
	88	AliHFEdetPIDqa::operator=(o);
	89	fHistos = o.fHistos;
	90
	91	return *this;
	92	}
	93
	94	//_________________________________________________________
	95	Long64_t AliHFEtrdPIDqaV1::Merge(TCollection *coll){
	96	//
	97	// Merge with other objects
	98	//
	99	if(!coll) return 0;
	100	if(coll->IsEmpty()) return 1;
	101
	102	TIter it(coll);
	103	AliHFEtrdPIDqaV1 *refQA = NULL;
	104	TObject *o = NULL;
	105	Long64_t count = 0;
	106	TList listHistos;
	107	while((o = it())){
	108	refQA = dynamic_cast<AliHFEtrdPIDqaV1 *>(o);
	109	if(!refQA) continue;
	110
	111	listHistos.Add(refQA->fHistos);
	112	count++;
113	}
114	fHistos->Merge(&listHistos);
115	return count + 1;
116	}
117
e3ae862b	118	//_________________________________________________________
	119	void AliHFEtrdPIDqaV1::Browse(TBrowser *b){
	120	//
	121	// Browse the PID QA
	122	//
	123	if(b){
	124	if(fHistos){
	125	b->Add(fHistos, fHistos->GetName());
	126
	127	// Make Projections of the dE/dx Spectra and add them to a new Folder
	128	TString specnames[4] = {"All", "Electrons", "Pions", "Protons"};
	129	Int_t specind[4] = {-1, AliPID::kElectron, AliPID::kPion, AliPID::kProton};
	130	TList *listTM = new TList;
	131	listTM->SetOwner();
	132	TList *listLike = new TList;
	133	listLike->SetOwner();
	134	TList *listCharge = new TList;
	135	listCharge->SetOwner();
	136	TList *listTPCnsigma = new TList;
	137	listTPCnsigma->SetOwner();
	138
	139	TH2 *hptr = NULL;
	140	for(Int_t ispec = 0; ispec < 4; ispec++){
	141	for(Int_t istep = 0; istep < 2; istep++){
	142	hptr = MakeTRDspectrumTM(static_cast<AliHFEdetPIDqa::EStep_t>(istep), specind[ispec]);
ccc37cdc	143	if(hptr){
	144	hptr->SetName(Form("hTRDtm%s%s", specnames[ispec].Data(), istep == 0 ? "Before" : "After"));
	145	listTM->Add(hptr);
	146	}
e3ae862b	147	hptr = MakeTRDlikelihoodDistribution(static_cast<AliHFEdetPIDqa::EStep_t>(istep), specind[ispec]);
	148	hptr->SetName(Form("hTRDlike%s%s", specnames[ispec].Data(), istep == 0 ? "Before" : "After"));
	149	listLike->Add(hptr);
	150	hptr = MakeTRDchargeDistribution(static_cast<AliHFEdetPIDqa::EStep_t>(istep), specind[ispec]);
	151	hptr->SetName(Form("hTRDcharge%s%s", specnames[ispec].Data(), istep == 0 ? "Before" : "After"));
	152	listCharge->Add(hptr);
	153	hptr = MakeTPCspectrumNsigma(static_cast<AliHFEdetPIDqa::EStep_t>(istep), specind[ispec]);
	154	hptr->SetName(Form("hTPCspectrum%s%s", specnames[ispec].Data(), istep == 0 ? "Before" : "After"));
	155	listTPCnsigma->Add(hptr);
	156	}
	157	}
	158
	159	b->Add(listTM, "Projections Truncated Mean");
	160	b->Add(listLike, "Projections Likelihood distribution");
	161	b->Add(listCharge, "Projections Tracklet Charge");
	162	b->Add(listTPCnsigma, "Projections TPC spectra");
	163	}
	164	}
	165	}
	166
3a72645a	167	//____________________________________________________________
	168	void AliHFEtrdPIDqaV1::Initialize(){
	169	//
	170	// Initialize QA histos for TRD PID
	171	//
	172
	173	AliDebug(1, "Initializing PID QA for TRD");
	174	// Make common binning
	175	const Int_t kPIDbins = AliPID::kSPECIES + 1;
	176	const Int_t kPbins = 100;
	177	const Int_t kSteps = 2;
	178	const Double_t kMinPID = -1;
	179	const Double_t kMinP = 0.;
	180	const Double_t kMaxPID = (Double_t)AliPID::kSPECIES;
	181	const Double_t kMaxP = 20.;
	182
	183	fHistos = new AliHFEcollection("trdqahistos", "Collection of TRD QA histograms");
	184
	185	// Create Control Histogram monitoring the TPC sigma between the selection steps
	186	const Int_t kTPCSigmaBins = 140;
	187	Int_t nBinsTPCSigma[4] = {kPIDbins, kPbins, kTPCSigmaBins, kSteps};
	188	Double_t minTPCSigma[4] = {kMinPID, kMinP, -12., 0};
	189	Double_t maxTPCSigma[4] = {kMaxPID, kMaxP, 12., 2.};
	190	fHistos->CreateTHnSparse("hTPCsigma", "TPC sigma; species p [GeV/c]; TPC dEdx - <dE/dx>\|_{el} [#sigma]; selection step", 4, nBinsTPCSigma, minTPCSigma, maxTPCSigma);
bf892a6a	191	fHistos->Sumw2("hTPCsigma");
3a72645a	192	// Create Monitoring histogram for the Likelihood distribution
	193	const Int_t kTRDLikelihoodBins = 100;
	194	Int_t nBinsTRDlike[4] = {kPIDbins, kPbins, kTRDLikelihoodBins, kSteps};
	195	Double_t minTRDlike[4] = {kMinPID, kMinP, 0., 0.};
	196	Double_t maxTRDlike[4] = {kMaxPID, kMaxP, 1., 2.};
	197	fHistos->CreateTHnSparse("hTRDlikelihood", "TRD Likelihood Distribution; species; p [GeV/c]; TRD electron Likelihood; selection step", 4, nBinsTRDlike, minTRDlike, maxTRDlike);
bf892a6a	198	fHistos->Sumw2("hTRDlikelihood");
3a72645a	199	// Create Monitoring histogram for the TRD total charge
	200	const Int_t kTRDchargeBins = 1000;
	201	Int_t nBinsTRDcharge[4] = {kPIDbins, kPbins, kTRDchargeBins, kSteps};
	202	Double_t minTRDcharge[4] = {kMinPID, kMinP, 0., 0.};
	203	Double_t maxTRDcharge[4] = {kMaxPID, kMaxP, 100000., 2.};
	204	fHistos->CreateTHnSparse("hTRDcharge", "Total TRD charge; species; p [GeV/c]; TRD charge [a.u.]; selection step", 4, nBinsTRDcharge, minTRDcharge, maxTRDcharge);
bf892a6a	205	fHistos->Sumw2("hTRDcharge");
e3ae862b	206	// Monitoring of the TRD truncated mean according to version 1
	207	const Int_t kTRDtmBins = 1000;
	208	Int_t nBinsTRDtm[4] = {kPIDbins, kPbins, kTRDtmBins, kSteps};
	209	Double_t minTRDtm[4] = {kMinPID, kMinP, 0., 0.};
	210	Double_t maxTRDtm[4] = {kMaxPID, kMaxP, 20000., 2.};
	211	fHistos->CreateTHnSparse("hTRDtruncatedMean", "TRD truncated Mean; species; p [GeV/c]; TRD signal [a.u.]; selection step", 4, nBinsTRDtm, minTRDtm, maxTRDtm);
ccc37cdc	212	fHistos->Sumw2("hTRDtruncatedMean");
3a72645a	213	}
	214
	215	//____________________________________________________________
6555e2ad	216	void AliHFEtrdPIDqaV1::ProcessTrack(const AliHFEpidObject *track, AliHFEdetPIDqa::EStep_t step){
3a72645a	217	//
	218	// Process the track, fill the containers
	219	//
	220	AliDebug(1, Form("QA started for TRD PID for step %d", (Int_t)step));
	221	Int_t species = track->GetAbInitioPID();
e3ae862b	222	if(species >= AliPID::kSPECIES) species = -1;
6555e2ad	223	AliHFEpidObject::AnalysisType_t anatype = track->IsESDanalysis() ? AliHFEpidObject::kESDanalysis : AliHFEpidObject::kAODanalysis;
3a72645a	224
6555e2ad	225	AliHFEpidTRD trdpid = dynamic_cast<AliHFEpidTRD >(fQAmanager->GetDetectorPID(AliHFEpid::kTRDpid));
	226	AliHFEpidTPC tpcpid = dynamic_cast<AliHFEpidTPC >(fQAmanager->GetDetectorPID(AliHFEpid::kTPCpid));
	227
3a72645a	228	Double_t container[4];
3a72645a	229	container[0] = species;
bf892a6a	230	container[1] = trdpid ? trdpid->GetP(track->GetRecTrack(), anatype) : 0.;
bf892a6a	231	container[2] = tpcpid ? tpcpid->NumberOfSigmas(track->GetRecTrack(), AliPID::kElectron, anatype) : 0.;
3a72645a	232	container[3] = step;
	233	fHistos->Fill("hTPCsigma", container);
	234
6555e2ad	235	container[2] = trdpid->GetElectronLikelihood(track->GetRecTrack(), anatype);
3a72645a	236	fHistos->Fill("hTRDlikelihood", container);
3a72645a	237
e3ae862b	238	if(track->IsESDanalysis()){
	239	container[2] = trdpid->GetTRDSignalV1(dynamic_cast<const AliESDtrack *>(track->GetRecTrack()));
	240	fHistos->Fill("hTRDtruncatedMean", container);
	241	}
6555e2ad	242	for(UInt_t ily = 0; ily < 6; ily++){
6555e2ad	243	container[2] = trdpid->GetChargeLayer(track->GetRecTrack(), ily, anatype);
3a72645a	244	fHistos->Fill("hTRDcharge", container);
3a72645a	245	}
3a72645a	246	}
	247
	248	//_________________________________________________________
	249	TH2 *AliHFEtrdPIDqaV1::MakeTPCspectrumNsigma(AliHFEdetPIDqa::EStep_t step, Int_t species){
	250	//
	251	// Get the TPC control histogram for the TRD selection step (either before or after PID)
	252	//
3a72645a	253	THnSparseF histo = dynamic_cast<THnSparseF >(fHistos->Get("hTPCsigma"));
	254	if(!histo){
	255	AliError("QA histogram monitoring TPC nSigma not available");
	256	return NULL;
	257	}
	258	if(species > -1 && species < AliPID::kSPECIES){
	259	// cut on species (if available)
	260	histo->GetAxis(0)->SetRange(species + 2, species + 2); // undef + underflow
	261	}
6555e2ad	262	histo->GetAxis(3)->SetRange(step + 1, step + 1);
3a72645a	263
	264	TH2 *hSpec = histo->Projection(2, 1);
	265	// construct title and name
6555e2ad	266	TString stepname = step == AliHFEdetPIDqa::kBeforePID ? "before" : "after";
3a72645a	267	TString speciesname = species > -1 && species < AliPID::kSPECIES ? AliPID::ParticleName(species) : "all Particles";
	268	TString specID = species > -1 && species < AliPID::kSPECIES ? AliPID::ParticleName(species) : "unid";
	269	TString histname = Form("hSigmaTPC%s%s", specID.Data(), stepname.Data());
	270	TString histtitle = Form("TPC Sigma for %s %s PID", speciesname.Data(), stepname.Data());
	271	hSpec->SetName(histname.Data());
	272	hSpec->SetTitle(histtitle.Data());
	273
	274	// Unset range on the original histogram
e3ae862b	275	histo->GetAxis(0)->SetRange(0, histo->GetAxis(0)->GetNbins());
	276	histo->GetAxis(2)->SetRange(0, histo->GetAxis(2)->GetNbins());
	277	return hSpec;
	278	}
	279
	280	//_________________________________________________________
	281	TH2 *AliHFEtrdPIDqaV1::MakeTRDspectrumTM(AliHFEdetPIDqa::EStep_t step, Int_t species){
	282	//
	283	// Get the TPC control histogram for the TRD selection step (either before or after PID)
	284	//
	285	THnSparseF histo = dynamic_cast<THnSparseF >(fHistos->Get("hTRDtruncatedMean"));
	286	if(!histo){
	287	AliError("QA histogram monitoring TPC nSigma not available");
	288	return NULL;
	289	}
	290	if(species > -1 && species < AliPID::kSPECIES){
	291	// cut on species (if available)
	292	histo->GetAxis(0)->SetRange(species + 2, species + 2); // undef + underflow
	293	}
	294	histo->GetAxis(3)->SetRange(step + 1, step + 1);
	295
	296	TH2 *hSpec = histo->Projection(2, 1);
	297	// construct title and name
	298	TString stepname = step == AliHFEdetPIDqa::kBeforePID ? "before" : "after";
	299	TString speciesname = species > -1 && species < AliPID::kSPECIES ? AliPID::ParticleName(species) : "all Particles";
	300	TString specID = species > -1 && species < AliPID::kSPECIES ? AliPID::ParticleName(species) : "unid";
	301	TString histname = Form("hTMTRD%s%s", specID.Data(), stepname.Data());
	302	TString histtitle = Form("TRD Truncated Mean for %s %s PID", speciesname.Data(), stepname.Data());
	303	hSpec->SetName(histname.Data());
	304	hSpec->SetTitle(histtitle.Data());
	305
	306	// Unset range on the original histogram
3a72645a	307	histo->GetAxis(0)->SetRange(0, histo->GetAxis(0)->GetNbins());
	308	histo->GetAxis(2)->SetRange(0, histo->GetAxis(2)->GetNbins());
	309	return hSpec;
	310	}
	311
	312	//_________________________________________________________
	313	TH2 *AliHFEtrdPIDqaV1::MakeTRDlikelihoodDistribution(AliHFEdetPIDqa::EStep_t step, Int_t species){
	314	//
	315	// Make Histogram for TRD Likelihood distribution
	316	//
	317	THnSparseF histo = dynamic_cast<THnSparseF >(fHistos->Get("hTRDlikelihood"));
	318	if(!histo){
	319	AliError("QA histogram monitoring TRD Electron Likelihood not available");
	320	return NULL;
	321	}
	322	if(species > -1 && species < AliPID::kSPECIES){
	323	// cut on species (if available)
	324	histo->GetAxis(0)->SetRange(species + 2, species + 2); // undef + underflow
	325	}
6555e2ad	326	histo->GetAxis(3)->SetRangeUser(step + 1, step + 1);
3a72645a	327
	328	TH2 *hSpec = histo->Projection(2, 1);
	329	// construct title and name
6555e2ad	330	TString stepname = step == AliHFEdetPIDqa::kBeforePID ? "before" : "after";
3a72645a	331	TString speciesname = species > -1 && species < AliPID::kSPECIES ? AliPID::ParticleName(species) : "all Particles";
	332	TString specID = species > -1 && species < AliPID::kSPECIES ? AliPID::ParticleName(species) : "unid";
	333	TString histname = Form("hLikeElTRD%s%s", specID.Data(), stepname.Data());
	334	TString histtitle = Form("TRD electron Likelihood for %s %s PID", speciesname.Data(), stepname.Data());
	335	hSpec->SetName(histname.Data());
	336	hSpec->SetTitle(histtitle.Data());
	337
	338	// Unset range on the original histogram
	339	histo->GetAxis(0)->SetRange(0, histo->GetAxis(0)->GetNbins());
	340	histo->GetAxis(2)->SetRange(0, histo->GetAxis(2)->GetNbins());
	341	return hSpec;
	342	}
	343
	344	//_________________________________________________________
	345	TH2 *AliHFEtrdPIDqaV1::MakeTRDchargeDistribution(AliHFEdetPIDqa::EStep_t step, Int_t species){
	346	//
	347	// Make Histogram for TRD Likelihood distribution
	348	//
	349	THnSparseF histo = dynamic_cast<THnSparseF >(fHistos->Get("hTRDcharge"));
	350	if(!histo){
	351	AliError("QA histogram monitoring TRD total charge not available");
	352	return NULL;
	353	}
	354	if(species > -1 && species < AliPID::kSPECIES){
	355	// cut on species (if available)
	356	histo->GetAxis(0)->SetRange(species + 2, species + 2); // undef + underflow
	357	}
6555e2ad	358	histo->GetAxis(3)->SetRange(step + 1, step + 1);
3a72645a	359
	360	TH2 *hSpec = histo->Projection(2, 1);
	361	// construct title and name
6555e2ad	362	TString stepname = step == AliHFEdetPIDqa::kBeforePID ? "before" : "after";
3a72645a	363	TString speciesname = species > -1 && species < AliPID::kSPECIES ? AliPID::ParticleName(species) : "all Particles";
	364	TString specID = species > -1 && species < AliPID::kSPECIES ? AliPID::ParticleName(species) : "unid";
	365	TString histname = Form("hChargeTRD%s%s", specID.Data(), stepname.Data());
	366	TString histtitle = Form("TRD total charge for %s %s PID", speciesname.Data(), stepname.Data());
	367	hSpec->SetName(histname.Data());
	368	hSpec->SetTitle(histtitle.Data());
	369
	370	// Unset range on the original histogram
	371	histo->GetAxis(0)->SetRange(0, histo->GetAxis(0)->GetNbins());
	372	histo->GetAxis(2)->SetRange(0, histo->GetAxis(2)->GetNbins());
	373	return hSpec;
	374	}
	375