[u/mrichter/AliRoot.git] / PWG3 / muon / AliAnalysisTaskTrigChEff.cxx

#define AliAnalysisTaskTrigChEff_cxx

// ROOT includes
#include "TChain.h"
#include "TH1.h"
#include "TCanvas.h"
#include "TROOT.h"
#include "TString.h"
#include "TList.h"

// STEER includes
#include "AliLog.h"

#include "AliESDEvent.h"
#include "AliESDMuonTrack.h"
#include "AliESDInputHandler.h"

#include "AliAODEvent.h"
#include "AliAODTrack.h"
#include "AliAODInputHandler.h"

// ANALYSIS includes
#include "AliAnalysisTask.h"
#include "AliAnalysisDataSlot.h"
#include "AliAnalysisManager.h"
#include "AliAnalysisTaskTrigChEff.h"

ClassImp(AliAnalysisTaskTrigChEff)

//________________________________________________________________________
AliAnalysisTaskTrigChEff::AliAnalysisTaskTrigChEff(const char *name) :
  AliAnalysisTask(name,""), 
  fESD(0),
  fAOD(0),
  fAnalysisType("ESD"),
  fList(0)
{
  //
  /// Constructor.
  //
  // Input slot #0 works with an Ntuple
  DefineInput(0, TChain::Class());
  // Output slot #0 writes into a TObjArray container
  DefineOutput(0,  TList::Class());
}

//___________________________________________________________________________
void AliAnalysisTaskTrigChEff::ConnectInputData(Option_t *) {
  //
  /// Connect ESD or AOD here
  /// Called once
  //

  TTree* tree = dynamic_cast<TTree*> (GetInputData(0));
  if (!tree) {
    Printf("ERROR: Could not read chain from input slot 0");
  } else {
    // Disable all branches and enable only the needed ones
    // The next two lines are different when data produced as AliESDEvent is read
    if(fAnalysisType == "ESD") {
      tree->SetBranchStatus("*", kFALSE);
      tree->SetBranchStatus("MuonTracks.*", kTRUE);

      AliESDInputHandler *esdH = dynamic_cast<AliESDInputHandler*> (AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler());
      
      if (!esdH) {
	Printf("ERROR: Could not get ESDInputHandler");
      } else
	fESD = esdH->GetEvent();
    }
    else if(fAnalysisType == "AOD") {
      AliAODInputHandler *aodH = dynamic_cast<AliAODInputHandler*> (AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler());
      
      if (!aodH) {
	Printf("ERROR: Could not get AODInputHandler");
      } else
	fAOD = aodH->GetEvent();
    }
    else 
      Printf("Wrong analysis type: Only ESD and AOD types are allowed!");
  }
}

//___________________________________________________________________________
void AliAnalysisTaskTrigChEff::CreateOutputObjects() {
  //
  /// Create histograms
  /// Called once
  //
  Printf("   CreateOutputObjects of task %s\n", GetName());

  TString cathCode[2] = {"bendPlane", "nonBendPlane"};
  TString countTypeName[2] = {"CountInCh", "NonCountInCh"};

  Char_t* yAxisTitle = "counts";

  const Int_t kNboards = 234; //AliMpConstants::NofLocalBoards();
  const Int_t kFirstTrigCh = 11;//AliMpConstants::NofTrackingChambers()+1;

  Int_t chamberBins = kNchambers;
  Float_t chamberLow = kFirstTrigCh-0.5, chamberHigh = kFirstTrigCh+kNchambers-0.5;
  Char_t* chamberName = "chamber";

  Int_t slatBins = kNslats;
  Float_t slatLow = 0-0.5, slatHigh = kNslats-0.5;
  Char_t* slatName = "slat";

  Int_t boardBins = kNboards;
  Float_t boardLow = 1-0.5, boardHigh = kNboards+1.-0.5;
  Char_t* boardName = "board";

  TString baseName, histoName;
  fList = new TList();

  TH1F* histo;

  histo = new TH1F("nTracksInSlat", "Num. of tracks used for efficiency calculation", 
		   slatBins, slatLow, slatHigh);
  histo->GetXaxis()->SetTitle(slatName);
  histo->GetYaxis()->SetTitle("num of used tracks");

  fList->AddAt(histo, kHtracksInSlat);

  histo = new TH1F("nTracksInBoard", "Num. of tracks used for efficiency calculation", 
		   boardBins, boardLow, boardHigh);
  histo->GetXaxis()->SetTitle(boardName);
  histo->GetYaxis()->SetTitle("num of used tracks");

  fList->AddAt(histo, kHtracksInBoard);

  for(Int_t hType=0; hType<kNcounts; hType++){
    Int_t hindex = (hType==0) ? kHchamberAllEff : kHchamberNonEff;
    for(Int_t cath=0; cath<kNcathodes; cath++){
      histoName = Form("%sChamber%s", cathCode[cath].Data(), countTypeName[hType].Data());
      histo = new TH1F(histoName, histoName,
		       chamberBins, chamberLow, chamberHigh);
      histo->GetXaxis()->SetTitle(chamberName);
      histo->GetYaxis()->SetTitle(yAxisTitle);
	
      fList->AddAt(histo, hindex + cath);
    } // loop on cath
  } // loop on counts

  for(Int_t hType=0; hType<kNcounts; hType++){
    Int_t hindex = (hType==0) ? kHslatAllEff : kHslatNonEff;
    for(Int_t cath=0; cath<kNcathodes; cath++){
      for(Int_t ch=0; ch<kNchambers; ch++){
	Int_t chCath = GetPlane(cath, ch);
	histoName = Form("%sSlat%s%i", cathCode[cath].Data(), countTypeName[hType].Data(), kFirstTrigCh+ch);
	histo = new TH1F(histoName, histoName,
			 slatBins, slatLow, slatHigh);
	histo->GetXaxis()->SetTitle(slatName);
	histo->GetYaxis()->SetTitle(yAxisTitle);

	fList->AddAt(histo, hindex + chCath);
      } // loop on chamber
    } // loop on cath
  } // loop on counts

  for(Int_t hType=0; hType<kNcounts; hType++){
    Int_t hindex = (hType==0) ? kHboardAllEff : kHboardNonEff;
    for(Int_t cath=0; cath<kNcathodes; cath++){
      for(Int_t ch=0; ch<kNchambers; ch++){
	Int_t chCath = GetPlane(cath, ch);
	histoName = Form("%sBoard%s%i", cathCode[cath].Data(), countTypeName[hType].Data(), kFirstTrigCh+ch);
	histo = new TH1F(histoName, histoName,
			 boardBins, boardLow, boardHigh);
	histo->GetXaxis()->SetTitle(boardName);
	histo->GetYaxis()->SetTitle(yAxisTitle);

	fList->AddAt(histo, hindex + chCath);
      } // loop on chamber
    } // loop on cath
  } // loop on counts
}

//________________________________________________________________________
void AliAnalysisTaskTrigChEff::Exec(Option_t *) {
  //
  /// Main loop
  /// Called for each event
  //
  Int_t nTracks = 0, board = 0;
  UShort_t pattern = 0;
  AliESDMuonTrack *esdTrack = 0x0;
  AliAODTrack* aodTrack = 0x0;

  if(fAnalysisType == "ESD") {
    if (!fESD) {
      Printf("ERROR: fESD not available");
      return;
    }
    nTracks = fESD->GetNumberOfMuonTracks(); 
  }
  else if(fAnalysisType == "AOD") {
    if (!fAOD) {
      Printf("ERROR: fAOD not available");
      return;
    }
    nTracks = fAOD->GetNumberOfTracks();
  }

  // Object declaration
  const Int_t kFirstTrigCh = 11; //AliMpConstants::NofTrackingChambers()+1;

  for (Int_t itrack = 0; itrack < nTracks; itrack++) {
    if(fAnalysisType == "ESD") {
      esdTrack = fESD->GetMuonTrack(itrack);
      pattern =  esdTrack->GetHitsPatternInTrigCh();
      board = esdTrack->LoCircuit();
    }
    else if(fAnalysisType == "AOD") {
      aodTrack = fAOD->GetTrack(itrack);
      if(!aodTrack->IsMuonTrack()) continue;
      pattern =  aodTrack->GetHitsPatternInTrigCh();
      board = 0; // aodTrack->LoCircuit(); Lo Circuit not implemented in AOD
    }

    Int_t effFlag = GetEffFlag(pattern);

    if(effFlag < kChEff) continue; // Track not good for efficiency calculation

    Int_t slat = GetSlat(pattern);

    if(effFlag >= kSlatEff) ((TH1F*)fList->At(kHtracksInSlat))->Fill(slat);
    if(effFlag >= kBoardEff) ((TH1F*)fList->At(kHtracksInBoard))->Fill(board);

    for(Int_t cath=0; cath<kNcathodes; cath++){
      Int_t ineffCh = IsChInefficient(pattern, cath);
      Int_t nChambers = kNchambers;
      for(Int_t ch=0; ch<nChambers; ch++){
	Int_t whichType = kAllChEff;
	Int_t currCh = ch;
	if(ineffCh>=0){
	  whichType = kChNonEff;
	  currCh = ineffCh;
	  nChambers = -1;
	}

	Int_t iChamber = kFirstTrigCh + currCh;
	Int_t hindex = (whichType==kAllChEff) ? kHchamberAllEff : kHchamberNonEff;
	((TH1F*)fList->At(hindex + cath))->Fill(iChamber);

	if(effFlag < kSlatEff) continue; // Track crossed different slats
	Int_t chCath = GetPlane(cath, currCh);
	hindex = (whichType==kAllChEff) ? kHslatAllEff : kHslatNonEff;
	((TH1F*)fList->At(hindex + chCath))->Fill(slat);

	if(effFlag < kBoardEff) continue; // Track crossed different boards
	hindex = (whichType==kAllChEff) ? kHboardAllEff : kHboardNonEff;
	((TH1F*)fList->At(hindex + chCath))->Fill(board);
      } // loop on chambers
    } // loop on cathodes
  }

  // Post final data. It will be written to a file with option "RECREATE"
  PostData(0, fList);
}

//________________________________________________________________________
void AliAnalysisTaskTrigChEff::Terminate(Option_t *) {
  //
  /// Draw result to the screen
  /// Called once at the end of the query.
  //
  if (!gROOT->IsBatch()) {
    TCanvas *can[kNcathodes];
    TH1F *num = 0x0;
    TH1F *den = 0x0;
    for(Int_t cath=0; cath<kNcathodes; cath++){
      TString canName = Form("can%i",cath);
      can[cath] = new TCanvas(canName.Data(),canName.Data(),10*(1+cath),10*(1+cath),310,310);
      can[cath]->SetFillColor(10); can[cath]->SetHighLightColor(10);
      can[cath]->SetLeftMargin(0.15); can[cath]->SetBottomMargin(0.15);  
      can[cath]->Divide(2,2);
      for(Int_t ch=0; ch<kNchambers; ch++){
	Int_t chCath = GetPlane(cath, ch);
	num = (TH1F*)(fList->At(kHboardAllEff + chCath)->Clone());
	den = (TH1F*)(fList->At(kHboardNonEff + chCath)->Clone());
	den->Add(num);
	num->Divide(den);
	can[cath]->cd(ch+1);
	num->DrawCopy("E");
      }
    }
  }
}

//________________________________________________________________________
Int_t AliAnalysisTaskTrigChEff::IsChInefficient(UShort_t pattern, 
						Int_t cathode)
{
  //
  /// Check which chamber was inefficient.
  //
  Int_t ineffCh = -999;
  for(Int_t ch=0; ch<kNchambers; ch++){
    Int_t chCath = GetPlane(cathode, ch);
    Int_t invert = kNplanes - chCath - 1;
    Int_t response = (pattern >> invert) & 0x01;
    if(!response) ineffCh = ch;
  }
  return ineffCh;
}
Commit	Line	Data
59156747	1	#define AliAnalysisTaskTrigChEff_cxx
	2
	3	// ROOT includes
	4	#include "TChain.h"
	5	#include "TH1.h"
	6	#include "TCanvas.h"
	7	#include "TROOT.h"
	8	#include "TString.h"
	9	#include "TList.h"
	10
	11	// STEER includes
	12	#include "AliLog.h"
	13
	14	#include "AliESDEvent.h"
	15	#include "AliESDMuonTrack.h"
	16	#include "AliESDInputHandler.h"
	17
	18	#include "AliAODEvent.h"
	19	#include "AliAODTrack.h"
	20	#include "AliAODInputHandler.h"
	21
	22	// ANALYSIS includes
	23	#include "AliAnalysisTask.h"
	24	#include "AliAnalysisDataSlot.h"
	25	#include "AliAnalysisManager.h"
	26	#include "AliAnalysisTaskTrigChEff.h"
	27
	28	ClassImp(AliAnalysisTaskTrigChEff)
	29
	30	//________________________________________________________________________
	31	AliAnalysisTaskTrigChEff::AliAnalysisTaskTrigChEff(const char *name) :
	32	AliAnalysisTask(name,""),
	33	fESD(0),
	34	fAOD(0),
	35	fAnalysisType("ESD"),
	36	fList(0)
	37	{
	38	//
	39	/// Constructor.
	40	//
	41	// Input slot #0 works with an Ntuple
	42	DefineInput(0, TChain::Class());
	43	// Output slot #0 writes into a TObjArray container
	44	DefineOutput(0, TList::Class());
	45	}
	46
	47	//___________________________________________________________________________
	48	void AliAnalysisTaskTrigChEff::ConnectInputData(Option_t *) {
	49	//
	50	/// Connect ESD or AOD here
	51	/// Called once
	52	//
	53
	54	TTree* tree = dynamic_cast<TTree*> (GetInputData(0));
	55	if (!tree) {
	56	Printf("ERROR: Could not read chain from input slot 0");
	57	} else {
	58	// Disable all branches and enable only the needed ones
	59	// The next two lines are different when data produced as AliESDEvent is read
	60	if(fAnalysisType == "ESD") {
	61	tree->SetBranchStatus("*", kFALSE);
	62	tree->SetBranchStatus("MuonTracks.*", kTRUE);
	63
	64	AliESDInputHandler esdH = dynamic_cast<AliESDInputHandler> (AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler());
65
66	if (!esdH) {
67	Printf("ERROR: Could not get ESDInputHandler");
68	} else
69	fESD = esdH->GetEvent();
70	}
71	else if(fAnalysisType == "AOD") {
72	AliAODInputHandler aodH = dynamic_cast<AliAODInputHandler> (AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler());
73
74	if (!aodH) {
75	Printf("ERROR: Could not get AODInputHandler");
76	} else
77	fAOD = aodH->GetEvent();
78	}
79	else
80	Printf("Wrong analysis type: Only ESD and AOD types are allowed!");
81	}
82	}
83
84	//___________________________________________________________________________
85	void AliAnalysisTaskTrigChEff::CreateOutputObjects() {
86	//
87	/// Create histograms
88	/// Called once
89	//
90	Printf(" CreateOutputObjects of task %s\n", GetName());
91
92	TString cathCode[2] = {"bendPlane", "nonBendPlane"};
93	TString countTypeName[2] = {"CountInCh", "NonCountInCh"};
94
95	Char_t* yAxisTitle = "counts";
96
97	const Int_t kNboards = 234; //AliMpConstants::NofLocalBoards();
98	const Int_t kFirstTrigCh = 11;//AliMpConstants::NofTrackingChambers()+1;
99
100	Int_t chamberBins = kNchambers;
101	Float_t chamberLow = kFirstTrigCh-0.5, chamberHigh = kFirstTrigCh+kNchambers-0.5;
102	Char_t* chamberName = "chamber";
103
104	Int_t slatBins = kNslats;
105	Float_t slatLow = 0-0.5, slatHigh = kNslats-0.5;
106	Char_t* slatName = "slat";
107
108	Int_t boardBins = kNboards;
109	Float_t boardLow = 1-0.5, boardHigh = kNboards+1.-0.5;
110	Char_t* boardName = "board";
111
112	TString baseName, histoName;
113	fList = new TList();
114
115	TH1F* histo;
116
117	histo = new TH1F("nTracksInSlat", "Num. of tracks used for efficiency calculation",
118	slatBins, slatLow, slatHigh);
119	histo->GetXaxis()->SetTitle(slatName);
120	histo->GetYaxis()->SetTitle("num of used tracks");
121
122	fList->AddAt(histo, kHtracksInSlat);
123
124	histo = new TH1F("nTracksInBoard", "Num. of tracks used for efficiency calculation",
125	boardBins, boardLow, boardHigh);
126	histo->GetXaxis()->SetTitle(boardName);
127	histo->GetYaxis()->SetTitle("num of used tracks");
128
129	fList->AddAt(histo, kHtracksInBoard);
130
131	for(Int_t hType=0; hType<kNcounts; hType++){
132	Int_t hindex = (hType==0) ? kHchamberAllEff : kHchamberNonEff;
133	for(Int_t cath=0; cath<kNcathodes; cath++){
134	histoName = Form("%sChamber%s", cathCode[cath].Data(), countTypeName[hType].Data());
135	histo = new TH1F(histoName, histoName,
136	chamberBins, chamberLow, chamberHigh);
137	histo->GetXaxis()->SetTitle(chamberName);
138	histo->GetYaxis()->SetTitle(yAxisTitle);
139
140	fList->AddAt(histo, hindex + cath);
141	} // loop on cath
142	} // loop on counts
143
144	for(Int_t hType=0; hType<kNcounts; hType++){
145	Int_t hindex = (hType==0) ? kHslatAllEff : kHslatNonEff;
146	for(Int_t cath=0; cath<kNcathodes; cath++){
147	for(Int_t ch=0; ch<kNchambers; ch++){
148	Int_t chCath = GetPlane(cath, ch);
149	histoName = Form("%sSlat%s%i", cathCode[cath].Data(), countTypeName[hType].Data(), kFirstTrigCh+ch);
150	histo = new TH1F(histoName, histoName,
151	slatBins, slatLow, slatHigh);
152	histo->GetXaxis()->SetTitle(slatName);
153	histo->GetYaxis()->SetTitle(yAxisTitle);
154
155	fList->AddAt(histo, hindex + chCath);
156	} // loop on chamber
157	} // loop on cath
158	} // loop on counts
159
160	for(Int_t hType=0; hType<kNcounts; hType++){
161	Int_t hindex = (hType==0) ? kHboardAllEff : kHboardNonEff;
162	for(Int_t cath=0; cath<kNcathodes; cath++){
163	for(Int_t ch=0; ch<kNchambers; ch++){
164	Int_t chCath = GetPlane(cath, ch);
165	histoName = Form("%sBoard%s%i", cathCode[cath].Data(), countTypeName[hType].Data(), kFirstTrigCh+ch);
166	histo = new TH1F(histoName, histoName,
167	boardBins, boardLow, boardHigh);
168	histo->GetXaxis()->SetTitle(boardName);
169	histo->GetYaxis()->SetTitle(yAxisTitle);
170
171	fList->AddAt(histo, hindex + chCath);
172	} // loop on chamber
173	} // loop on cath
174	} // loop on counts
175	}
176
177	//________________________________________________________________________
178	void AliAnalysisTaskTrigChEff::Exec(Option_t *) {
179	//
180	/// Main loop
181	/// Called for each event
182	//
183	Int_t nTracks = 0, board = 0;
184	UShort_t pattern = 0;
185	AliESDMuonTrack *esdTrack = 0x0;
186	AliAODTrack* aodTrack = 0x0;
187
188	if(fAnalysisType == "ESD") {
189	if (!fESD) {
190	Printf("ERROR: fESD not available");
191	return;
192	}
193	nTracks = fESD->GetNumberOfMuonTracks();
194	}
195	else if(fAnalysisType == "AOD") {
196	if (!fAOD) {
197	Printf("ERROR: fAOD not available");
198	return;
199	}
200	nTracks = fAOD->GetNumberOfTracks();
201	}
202
203	// Object declaration
204	const Int_t kFirstTrigCh = 11; //AliMpConstants::NofTrackingChambers()+1;
205
206	for (Int_t itrack = 0; itrack < nTracks; itrack++) {
207	if(fAnalysisType == "ESD") {
208	esdTrack = fESD->GetMuonTrack(itrack);
209	pattern = esdTrack->GetHitsPatternInTrigCh();
210	board = esdTrack->LoCircuit();
211	}
212	else if(fAnalysisType == "AOD") {
213	aodTrack = fAOD->GetTrack(itrack);
214	if(!aodTrack->IsMuonTrack()) continue;
215	pattern = aodTrack->GetHitsPatternInTrigCh();
216	board = 0; // aodTrack->LoCircuit(); Lo Circuit not implemented in AOD
217	}
218
219	Int_t effFlag = GetEffFlag(pattern);
220
221	if(effFlag < kChEff) continue; // Track not good for efficiency calculation
222
223	Int_t slat = GetSlat(pattern);
224
225	if(effFlag >= kSlatEff) ((TH1F*)fList->At(kHtracksInSlat))->Fill(slat);
226	if(effFlag >= kBoardEff) ((TH1F*)fList->At(kHtracksInBoard))->Fill(board);
227
228	for(Int_t cath=0; cath<kNcathodes; cath++){
229	Int_t ineffCh = IsChInefficient(pattern, cath);
230	Int_t nChambers = kNchambers;
231	for(Int_t ch=0; ch<nChambers; ch++){
232	Int_t whichType = kAllChEff;
233	Int_t currCh = ch;
234	if(ineffCh>=0){
235	whichType = kChNonEff;
236	currCh = ineffCh;
237	nChambers = -1;
238	}
239
240	Int_t iChamber = kFirstTrigCh + currCh;
241	Int_t hindex = (whichType==kAllChEff) ? kHchamberAllEff : kHchamberNonEff;
242	((TH1F*)fList->At(hindex + cath))->Fill(iChamber);
243
244	if(effFlag < kSlatEff) continue; // Track crossed different slats
245	Int_t chCath = GetPlane(cath, currCh);
246	hindex = (whichType==kAllChEff) ? kHslatAllEff : kHslatNonEff;
247	((TH1F*)fList->At(hindex + chCath))->Fill(slat);
248
249	if(effFlag < kBoardEff) continue; // Track crossed different boards
250	hindex = (whichType==kAllChEff) ? kHboardAllEff : kHboardNonEff;
251	((TH1F*)fList->At(hindex + chCath))->Fill(board);
252	} // loop on chambers
253	} // loop on cathodes
254	}
255
256	// Post final data. It will be written to a file with option "RECREATE"
257	PostData(0, fList);
258	}
259
260	//________________________________________________________________________
261	void AliAnalysisTaskTrigChEff::Terminate(Option_t *) {
262	//
263	/// Draw result to the screen
264	/// Called once at the end of the query.
265	//
266	if (!gROOT->IsBatch()) {
267	TCanvas *can[kNcathodes];
268	TH1F *num = 0x0;
269	TH1F *den = 0x0;
270	for(Int_t cath=0; cath<kNcathodes; cath++){
271	TString canName = Form("can%i",cath);
272	can[cath] = new TCanvas(canName.Data(),canName.Data(),10(1+cath),10(1+cath),310,310);
273	can[cath]->SetFillColor(10); can[cath]->SetHighLightColor(10);
274	can[cath]->SetLeftMargin(0.15); can[cath]->SetBottomMargin(0.15);
275	can[cath]->Divide(2,2);
276	for(Int_t ch=0; ch<kNchambers; ch++){
277	Int_t chCath = GetPlane(cath, ch);
278	num = (TH1F*)(fList->At(kHboardAllEff + chCath)->Clone());
279	den = (TH1F*)(fList->At(kHboardNonEff + chCath)->Clone());
280	den->Add(num);
281	num->Divide(den);
282	can[cath]->cd(ch+1);
283	num->DrawCopy("E");
284	}
285	}
286	}
287	}
288
289	//________________________________________________________________________
290	Int_t AliAnalysisTaskTrigChEff::IsChInefficient(UShort_t pattern,
291	Int_t cathode)
292	{
293	//
294	/// Check which chamber was inefficient.
295	//
296	Int_t ineffCh = -999;
297	for(Int_t ch=0; ch<kNchambers; ch++){
298	Int_t chCath = GetPlane(cathode, ch);
299	Int_t invert = kNplanes - chCath - 1;
300	Int_t response = (pattern >> invert) & 0x01;
301	if(!response) ineffCh = ch;
302	}
303	return ineffCh;
304	}