[u/mrichter/AliRoot.git] / PWGCF / EBYE / BalanceFunctions / AliAnalysisTaskAODFilterBitQA.cxx

#include "TChain.h"
#include "TH2D.h"
#include "TH3D.h"


#include "AliAODEvent.h"
#include "AliAODTrack.h"
#include "AliLog.h"

#include "AliAnalysisTaskSE.h"
#include "AliMCEvent.h"
#include "AliAODMCParticle.h"


#include "AliAnalysisTaskAODFilterBitQA.h"

// Analysis task for the QA of AOD track filter bits
// Authors: m.weber@cern.ch

ClassImp(AliAnalysisTaskAODFilterBitQA)

//________________________________________________________________________
AliAnalysisTaskAODFilterBitQA::AliAnalysisTaskAODFilterBitQA(const char *name) 
  : AliAnalysisTaskSE(name),
  fArrayMC(0x0),
  fListQA(0x0),
  fillOnlySecondaries(kFALSE),
  fCentralityPercentileMin(0.),
  fCentralityPercentileMax(80.), 
  fPtMin(0),
  fPtMax(1000),
  fEtaMin(-10),
  fEtaMax(10),
  fHistTrackStats(0)
{

  for(Int_t iCharge = 0; iCharge < gNCharge; iCharge++){
    for(Int_t iTrackBit = 0; iTrackBit < gBitMax; iTrackBit++){
      fHistKinematics[iCharge][iTrackBit] = NULL;
      fHistDCAconstrained[iCharge][iTrackBit] = NULL;
      fHistDCAglobal[iCharge][iTrackBit]  = NULL;
      fHistChiClus[iCharge][iTrackBit]    = NULL;
    }
  }
  
  DefineInput(0, TChain::Class());
  // Output slot #0 writes into a TH1 container
  DefineOutput(1, TList::Class());

}

//________________________________________________________________________
AliAnalysisTaskAODFilterBitQA::~AliAnalysisTaskAODFilterBitQA() {
  // Destructor
  // ... not implemented
}

//________________________________________________________________________
void AliAnalysisTaskAODFilterBitQA::UserCreateOutputObjects() {
  // Create histograms
  // Called once

   // global switch disabling the reference 
  // (to avoid "Replacing existing TH1" if several wagons are created in train)
  Bool_t oldStatus = TH1::AddDirectoryStatus();
  TH1::AddDirectory(kFALSE);

  // QA list
  fListQA = new TList();
  fListQA->SetName("listQA");
  fListQA->SetOwner();

  // QA histograms
  fHistTrackStats = new TH2D("fHistTrackStats","Track statistics;Centrality;TrackFilterBit;N_{events}",100,0,100,gBitMax,0,gBitMax);
  fListQA->Add(fHistTrackStats);

  TString sCharge[gNCharge] = {"Plus","Minus"};
  
  for(Int_t iCharge = 0; iCharge < gNCharge; iCharge++){
    for(Int_t iTrackBit = 0; iTrackBit < gBitMax; iTrackBit++){
      fHistKinematics[iCharge][iTrackBit] = new TH3D(Form("Bit%d_%s_Kinematics",iTrackBit,sCharge[iCharge].Data()),Form("Bit%d_%s_Kinematics;#eta;#varphi (rad);p_{T} (GeV/c)",iTrackBit,sCharge[iCharge].Data()),100,-1.0,1.0,100,0,TMath::Pi()*2,100,0,10);
      fHistDCAconstrained[iCharge][iTrackBit] = new TH2D(Form("Bit%d_%s_DCAconstrained",iTrackBit,sCharge[iCharge].Data()),Form("Bit%d_%s_DCAconstrained;DCA XY [Constrained] (cm);DCA Z [Constrained] (cm)",iTrackBit,sCharge[iCharge].Data()),100,-5.0,5.0,100,-5.0,5.0);
      fHistDCAglobal[iCharge][iTrackBit]  = new TH3D(Form("Bit%d_%s_DCAglobal",iTrackBit,sCharge[iCharge].Data()),Form("Bit%d_%s_DCAglobal;DCA X [Global] (cm);DCA Y [Global] (cm);DCA Z [Global] (cm)",iTrackBit,sCharge[iCharge].Data()),100,-5.0,5.0,100,-5.0,5.0,100,-5.0,5.0);
fHistChiClus[iCharge][iTrackBit]    = new TH2D(Form("Bit%d_%s_ChiClus",iTrackBit,sCharge[iCharge].Data()),Form("Bit%d_%s_ChiClus;#chi^{2} [Fit];N_{clus} [TPC]",iTrackBit,sCharge[iCharge].Data()),100,-1.0,5.0,160,0,160.0);
      fListQA->Add(fHistKinematics[iCharge][iTrackBit]);
      fListQA->Add(fHistDCAconstrained[iCharge][iTrackBit]);
      fListQA->Add(fHistDCAglobal[iCharge][iTrackBit]);
      fListQA->Add(fHistChiClus[iCharge][iTrackBit]);
    }
  }

  // Post output data.
  PostData(1, fListQA);

  AliInfo("Finished setting up the Output");
  TH1::AddDirectory(oldStatus);
}

//________________________________________________________________________
void AliAnalysisTaskAODFilterBitQA::UserExec(Option_t *) {
  // Main loop
  // Called for each event

  AliVEvent* event = dynamic_cast<AliVEvent*>(InputEvent());     
  if(!event) {
    AliError("event not available");
    return;
  }

  // MC information (set if available)
  fArrayMC = dynamic_cast<TClonesArray*>(event->FindListObject(AliAODMCParticle::StdBranchName()));   
    
  // check event cuts
  Double_t lMultiplicityVar = -1;
  if((lMultiplicityVar = IsEventAccepted(event)) < 0){ 
    return;
  }

  
  // get the accepted tracks in main event
  GetAcceptedTracks(event,lMultiplicityVar);

}

//________________________________________________________________________
void  AliAnalysisTaskAODFilterBitQA::FinishTaskOutput(){
  // Finish task output
  // not implemented ...

}

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

}


//________________________________________________________________________
Double_t AliAnalysisTaskAODFilterBitQA::IsEventAccepted(AliVEvent *event){
  // Checks the Event cuts

  // still hard coded
  Double_t fVxMax = 0.5;
  Double_t fVyMax = 0.5;
  Double_t fVzMax = 10.0;
  TString fCentralityEstimator = "V0M";

  Float_t gCentrality      = -1.;
  const AliVVertex *vertex = event->GetPrimaryVertex();
  
  if(vertex) {
    Double32_t fCov[6];
    vertex->GetCovarianceMatrix(fCov);
    if(vertex->GetNContributors() > 0) {
      if(fCov[5] != 0) {
	if(TMath::Abs(vertex->GetX()) < fVxMax) {
	  if(TMath::Abs(vertex->GetY()) < fVyMax) {
	    if(TMath::Abs(vertex->GetZ()) < fVzMax) {
	      
	      // get the reference multiplicty or centrality
	      AliAODHeader *header = (AliAODHeader*) event->GetHeader();
	      gCentrality = header->GetCentralityP()->GetCentralityPercentile(fCentralityEstimator.Data());
	      
	      if((gCentrality > fCentralityPercentileMin) && (gCentrality < fCentralityPercentileMax)){
		
		return gCentrality;
		
	      }//centrality range
	    }//Vz cut
	  }//Vy cut
	}//Vx cut
      }//proper vertex resolution
    }//proper number of contributors
  }//vertex object valid
  
  // in all other cases return -1 (event not accepted)
  return -1;
}

//________________________________________________________________________
void AliAnalysisTaskAODFilterBitQA::GetAcceptedTracks(AliVEvent *event, Double_t gCentrality){
  // Checks track cuts (filter bits)
  // Fills QA histograms


  Double_t vEta;
  Double_t vPhi;
  Double_t vPt;
  Double_t vDCAconstrainedxy;
  Double_t vDCAconstrainedz; 
  Double_t vDCAglobalx;
  Double_t vDCAglobaly;
  Double_t vDCAglobalz;
  Double_t vChi2;
  Double_t vClus;

  Double_t pos[3];
  Double_t v[3];

  const AliVVertex *vertex = event->GetPrimaryVertex();
  vertex->GetXYZ(v);

  // Loop over tracks in event
  for (Int_t iTracks = 0; iTracks < event->GetNumberOfTracks(); iTracks++) {
    AliAODTrack* aodTrack = dynamic_cast<AliAODTrack *>(event->GetTrack(iTracks));
    if (!aodTrack) {
      AliError(Form("Could not receive track %d", iTracks));
      continue;
    }

    // get MC information (if available)
    if(fArrayMC && fillOnlySecondaries){
      
      Int_t label = aodTrack->GetLabel();
      AliAODMCParticle *mcTrack = (AliAODMCParticle *)fArrayMC->At(TMath::Abs(label));

      if(mcTrack->IsPhysicalPrimary())
	continue;      
    }

    // track parameters
    vEta    = aodTrack->Eta();
    vPhi    = aodTrack->Phi();// * TMath::RadToDeg();
    vPt     = aodTrack->Pt();
    vDCAconstrainedxy  = aodTrack->DCA();     
    vDCAconstrainedz   = aodTrack->ZAtDCA();   
    vChi2   = aodTrack->Chi2perNDF(); 
    vClus   = aodTrack->GetTPCNcls(); 

    // kinematic cuts
    if( vPt > fPtMax || vPt < fPtMin )
      continue;
    if( vEta > fEtaMax || vEta < fEtaMin )
      continue;

    // if not constrained track the position is stored (primary vertex to be subtracted)
    aodTrack->GetXYZ(pos);
    vDCAglobalx  = pos[0] - v[0];
    vDCAglobaly  = pos[1] - v[1];
    vDCAglobalz  = pos[2] - v[2];
    
    // fill for separately for positive and negative charges
    Int_t iCharge = -1;
    // positive 
    if(aodTrack->Charge() > 0)
      iCharge = 0;
    else if(aodTrack->Charge() < 0)
      iCharge = 1;
    else{
      AliError("Charge==0?");
      iCharge = -1;
    }


    // AOD track cuts
    if(iCharge > -1){
      for(Int_t iTrackBit = 0; iTrackBit < gBitMax; iTrackBit++){
	fHistTrackStats->Fill(gCentrality,iTrackBit,aodTrack->TestFilterBit(1<<iTrackBit));
	
	if(aodTrack->TestFilterBit(1<<iTrackBit)){
	  fHistKinematics[iCharge][iTrackBit]->Fill(vEta,vPhi,vPt);
	  fHistDCAconstrained[iCharge][iTrackBit]->Fill(vDCAconstrainedxy,vDCAconstrainedz);
	  fHistDCAglobal[iCharge][iTrackBit]->Fill(vDCAglobalx,vDCAglobaly,vDCAglobalz);
	  fHistChiClus[iCharge][iTrackBit]->Fill(vChi2,vClus);
	} 
      }//bit loop
    }//charge positive or negative
  }//track loop
  
  return;  
}
Commit	Line	Data
bf6a8159	1	#include "TChain.h"
	2	#include "TH2D.h"
	3	#include "TH3D.h"
	4
	5
	6	#include "AliAODEvent.h"
	7	#include "AliAODTrack.h"
	8	#include "AliLog.h"
	9
a3d2ee41	10	#include "AliAnalysisTaskSE.h"
	11	#include "AliMCEvent.h"
	12	#include "AliAODMCParticle.h"
	13
	14
bf6a8159	15
	16	#include "AliAnalysisTaskAODFilterBitQA.h"
	17
	18	// Analysis task for the QA of AOD track filter bits
	19	// Authors: m.weber@cern.ch
	20
	21	ClassImp(AliAnalysisTaskAODFilterBitQA)
	22
	23	//________________________________________________________________________
	24	AliAnalysisTaskAODFilterBitQA::AliAnalysisTaskAODFilterBitQA(const char *name)
	25	: AliAnalysisTaskSE(name),
a3d2ee41	26	fArrayMC(0x0),
	27	fListQA(0x0),
	28	fillOnlySecondaries(kFALSE),
	29	fCentralityPercentileMin(0.),
	30	fCentralityPercentileMax(80.),
	31	fPtMin(0),
	32	fPtMax(1000),
	33	fEtaMin(-10),
	34	fEtaMax(10),
bf6a8159	35	fHistTrackStats(0)
	36	{
	37
a3d2ee41	38	for(Int_t iCharge = 0; iCharge < gNCharge; iCharge++){
	39	for(Int_t iTrackBit = 0; iTrackBit < gBitMax; iTrackBit++){
	40	fHistKinematics[iCharge][iTrackBit] = NULL;
	41	fHistDCAconstrained[iCharge][iTrackBit] = NULL;
	42	fHistDCAglobal[iCharge][iTrackBit] = NULL;
	43	fHistChiClus[iCharge][iTrackBit] = NULL;
	44	}
bf6a8159	45	}
	46
	47	DefineInput(0, TChain::Class());
	48	// Output slot #0 writes into a TH1 container
	49	DefineOutput(1, TList::Class());
	50
	51	}
	52
	53	//________________________________________________________________________
	54	AliAnalysisTaskAODFilterBitQA::~AliAnalysisTaskAODFilterBitQA() {
	55	// Destructor
	56	// ... not implemented
	57	}
	58
	59	//________________________________________________________________________
	60	void AliAnalysisTaskAODFilterBitQA::UserCreateOutputObjects() {
	61	// Create histograms
	62	// Called once
	63
	64	// global switch disabling the reference
	65	// (to avoid "Replacing existing TH1" if several wagons are created in train)
	66	Bool_t oldStatus = TH1::AddDirectoryStatus();
	67	TH1::AddDirectory(kFALSE);
	68
	69	// QA list
	70	fListQA = new TList();
	71	fListQA->SetName("listQA");
	72	fListQA->SetOwner();
	73
	74	// QA histograms
	75	fHistTrackStats = new TH2D("fHistTrackStats","Track statistics;Centrality;TrackFilterBit;N_{events}",100,0,100,gBitMax,0,gBitMax);
	76	fListQA->Add(fHistTrackStats);
	77
a3d2ee41	78	TString sCharge[gNCharge] = {"Plus","Minus"};
	79
	80	for(Int_t iCharge = 0; iCharge < gNCharge; iCharge++){
	81	for(Int_t iTrackBit = 0; iTrackBit < gBitMax; iTrackBit++){
	82	fHistKinematics[iCharge][iTrackBit] = new TH3D(Form("Bit%d_%s_Kinematics",iTrackBit,sCharge[iCharge].Data()),Form("Bit%d_%s_Kinematics;#eta;#varphi (rad);p_{T} (GeV/c)",iTrackBit,sCharge[iCharge].Data()),100,-1.0,1.0,100,0,TMath::Pi()*2,100,0,10);
	83	fHistDCAconstrained[iCharge][iTrackBit] = new TH2D(Form("Bit%d_%s_DCAconstrained",iTrackBit,sCharge[iCharge].Data()),Form("Bit%d_%s_DCAconstrained;DCA XY [Constrained] (cm);DCA Z [Constrained] (cm)",iTrackBit,sCharge[iCharge].Data()),100,-5.0,5.0,100,-5.0,5.0);
	84	fHistDCAglobal[iCharge][iTrackBit] = new TH3D(Form("Bit%d_%s_DCAglobal",iTrackBit,sCharge[iCharge].Data()),Form("Bit%d_%s_DCAglobal;DCA X [Global] (cm);DCA Y [Global] (cm);DCA Z [Global] (cm)",iTrackBit,sCharge[iCharge].Data()),100,-5.0,5.0,100,-5.0,5.0,100,-5.0,5.0);
	85	fHistChiClus[iCharge][iTrackBit] = new TH2D(Form("Bit%d_%s_ChiClus",iTrackBit,sCharge[iCharge].Data()),Form("Bit%d_%s_ChiClus;#chi^{2} [Fit];N_{clus} [TPC]",iTrackBit,sCharge[iCharge].Data()),100,-1.0,5.0,160,0,160.0);
	86	fListQA->Add(fHistKinematics[iCharge][iTrackBit]);
	87	fListQA->Add(fHistDCAconstrained[iCharge][iTrackBit]);
	88	fListQA->Add(fHistDCAglobal[iCharge][iTrackBit]);
	89	fListQA->Add(fHistChiClus[iCharge][iTrackBit]);
	90	}
bf6a8159	91	}
	92
	93	// Post output data.
	94	PostData(1, fListQA);
	95
	96	AliInfo("Finished setting up the Output");
	97	TH1::AddDirectory(oldStatus);
	98	}
	99
	100	//________________________________________________________________________
	101	void AliAnalysisTaskAODFilterBitQA::UserExec(Option_t *) {
	102	// Main loop
	103	// Called for each event
	104
	105	AliVEvent* event = dynamic_cast<AliVEvent*>(InputEvent());
	106	if(!event) {
	107	AliError("event not available");
	108	return;
	109	}
	110
a3d2ee41	111	// MC information (set if available)
	112	fArrayMC = dynamic_cast<TClonesArray*>(event->FindListObject(AliAODMCParticle::StdBranchName()));
	113
bf6a8159	114	// check event cuts
	115	Double_t lMultiplicityVar = -1;
	116	if((lMultiplicityVar = IsEventAccepted(event)) < 0){
	117	return;
	118	}
	119
	120
	121	// get the accepted tracks in main event
	122	GetAcceptedTracks(event,lMultiplicityVar);
	123
	124	}
	125
	126	//________________________________________________________________________
	127	void AliAnalysisTaskAODFilterBitQA::FinishTaskOutput(){
	128	// Finish task output
	129	// not implemented ...
	130
	131	}
	132
	133	//________________________________________________________________________
	134	void AliAnalysisTaskAODFilterBitQA::Terminate(Option_t *) {
	135	// Draw result to the screen
	136	// Called once at the end of the query
	137	// not implemented ...
	138
	139	}
	140
	141
	142	//________________________________________________________________________
	143	Double_t AliAnalysisTaskAODFilterBitQA::IsEventAccepted(AliVEvent *event){
	144	// Checks the Event cuts
	145
	146	// still hard coded
bf6a8159	147	Double_t fVxMax = 0.5;
	148	Double_t fVyMax = 0.5;
	149	Double_t fVzMax = 10.0;
	150	TString fCentralityEstimator = "V0M";
	151
	152	Float_t gCentrality = -1.;
	153	const AliVVertex *vertex = event->GetPrimaryVertex();
	154
	155	if(vertex) {
	156	Double32_t fCov[6];
	157	vertex->GetCovarianceMatrix(fCov);
	158	if(vertex->GetNContributors() > 0) {
	159	if(fCov[5] != 0) {
	160	if(TMath::Abs(vertex->GetX()) < fVxMax) {
	161	if(TMath::Abs(vertex->GetY()) < fVyMax) {
	162	if(TMath::Abs(vertex->GetZ()) < fVzMax) {
	163
	164	// get the reference multiplicty or centrality
	165	AliAODHeader header = (AliAODHeader) event->GetHeader();
	166	gCentrality = header->GetCentralityP()->GetCentralityPercentile(fCentralityEstimator.Data());
	167
	168	if((gCentrality > fCentralityPercentileMin) && (gCentrality < fCentralityPercentileMax)){
	169
	170	return gCentrality;
	171
	172	}//centrality range
	173	}//Vz cut
	174	}//Vy cut
	175	}//Vx cut
	176	}//proper vertex resolution
	177	}//proper number of contributors
	178	}//vertex object valid
	179
	180	// in all other cases return -1 (event not accepted)
	181	return -1;
	182	}
	183
	184	//________________________________________________________________________
	185	void AliAnalysisTaskAODFilterBitQA::GetAcceptedTracks(AliVEvent *event, Double_t gCentrality){
	186	// Checks track cuts (filter bits)
	187	// Fills QA histograms
	188
	189
bf6a8159	190	Double_t vEta;
bf6a8159	191	Double_t vPhi;
bf6a8159	192	Double_t vPt;
ddf67b83	193	Double_t vDCAconstrainedxy;
	194	Double_t vDCAconstrainedz;
	195	Double_t vDCAglobalx;
	196	Double_t vDCAglobaly;
	197	Double_t vDCAglobalz;
bf6a8159	198	Double_t vChi2;
	199	Double_t vClus;
	200
ddf67b83	201	Double_t pos[3];
	202	Double_t v[3];
	203
bf6a8159	204	const AliVVertex *vertex = event->GetPrimaryVertex();
ddf67b83	205	vertex->GetXYZ(v);
bf6a8159	206
	207	// Loop over tracks in event
	208	for (Int_t iTracks = 0; iTracks < event->GetNumberOfTracks(); iTracks++) {
	209	AliAODTrack* aodTrack = dynamic_cast<AliAODTrack *>(event->GetTrack(iTracks));
	210	if (!aodTrack) {
	211	AliError(Form("Could not receive track %d", iTracks));
	212	continue;
	213	}
	214
a3d2ee41	215	// get MC information (if available)
	216	if(fArrayMC && fillOnlySecondaries){
	217
	218	Int_t label = aodTrack->GetLabel();
	219	AliAODMCParticle mcTrack = (AliAODMCParticle )fArrayMC->At(TMath::Abs(label));
	220
	221	if(mcTrack->IsPhysicalPrimary())
	222	continue;
	223	}
	224
bf6a8159	225	// track parameters
bf6a8159	226	vEta = aodTrack->Eta();
bf6a8159	227	vPhi = aodTrack->Phi();// * TMath::RadToDeg();
bf6a8159	228	vPt = aodTrack->Pt();
ddf67b83	229	vDCAconstrainedxy = aodTrack->DCA();
ddf67b83	230	vDCAconstrainedz = aodTrack->ZAtDCA();
bf6a8159	231	vChi2 = aodTrack->Chi2perNDF();
	232	vClus = aodTrack->GetTPCNcls();
	233
a3d2ee41	234	// kinematic cuts
	235	if( vPt > fPtMax \|\| vPt < fPtMin )
	236	continue;
	237	if( vEta > fEtaMax \|\| vEta < fEtaMin )
	238	continue;
	239
ddf67b83	240	// if not constrained track the position is stored (primary vertex to be subtracted)
	241	aodTrack->GetXYZ(pos);
	242	vDCAglobalx = pos[0] - v[0];
	243	vDCAglobaly = pos[1] - v[1];
	244	vDCAglobalz = pos[2] - v[2];
bf6a8159	245
a3d2ee41	246	// fill for separately for positive and negative charges
	247	Int_t iCharge = -1;
	248	// positive
	249	if(aodTrack->Charge() > 0)
	250	iCharge = 0;
	251	else if(aodTrack->Charge() < 0)
	252	iCharge = 1;
	253	else{
	254	AliError("Charge==0?");
	255	iCharge = -1;
	256	}
	257
ddf67b83	258
bf6a8159	259	// AOD track cuts
d0b37020	260	if(iCharge > -1){
	261	for(Int_t iTrackBit = 0; iTrackBit < gBitMax; iTrackBit++){
	262	fHistTrackStats->Fill(gCentrality,iTrackBit,aodTrack->TestFilterBit(1<<iTrackBit));
	263
	264	if(aodTrack->TestFilterBit(1<<iTrackBit)){
	265	fHistKinematics[iCharge][iTrackBit]->Fill(vEta,vPhi,vPt);
	266	fHistDCAconstrained[iCharge][iTrackBit]->Fill(vDCAconstrainedxy,vDCAconstrainedz);
	267	fHistDCAglobal[iCharge][iTrackBit]->Fill(vDCAglobalx,vDCAglobaly,vDCAglobalz);
	268	fHistChiClus[iCharge][iTrackBit]->Fill(vChi2,vClus);
	269	}
	270	}//bit loop
	271	}//charge positive or negative
bf6a8159	272	}//track loop
d0b37020	273
bf6a8159	274	return;
bf6a8159	275	}