[u/mrichter/AliRoot.git] / PWGLF / FORWARD / photons / AliAnalysisTaskPMD.cxx

#include "TChain.h"
#include "TTree.h"
#include "TH1F.h"
#include "TH2F.h"
#include "TCanvas.h"
#include "AliAnalysisTask.h"
#include "AliAnalysisManager.h"
#include "AliESDEvent.h"
#include "AliESDInputHandler.h"
#include "AliESDPmdTrack.h"
#include "AliESDVertex.h"
#include "AliESDVZERO.h"
#include "AliAnalysisTaskPMD.h"

// AnalysisTask For PMD
// Authors: Sudipan De, Subhash Singha

ClassImp(AliAnalysisTaskPMD)

//________________________________________________________________________
AliAnalysisTaskPMD::AliAnalysisTaskPMD(const char *name) 
: AliAnalysisTaskSE(name), 
  fESD(0), 
  fOutputList(0), 
  fHistTotEvent(0),
  fHistTotEventAfterPhySel(0),
  fHistTotEventAfterVtx(0),
  fHistVtxZ(0),
  fHistXYPre(0),
  fHistEta(0),
  fHistEta1(0),
  fHistMultMeas(0),
  fHistMultMeas1(0)
{
  for(Int_t i=0; i<10; i++){
    fHistMultMeasEtaBinA[i] = 0;
    fHistMultMeasEtaBinA1[i] = 0;
  }
  // Constructor
  
  // Define input and output slots here
  // Input slot #0 works with a TChain
  DefineInput(0, TChain::Class());
  // Output slot #0 id reserved by the base class for AOD
  // Output slot #1 writes into a TH1 container
  DefineOutput(1, TList::Class());
}

//________________________________________________________________________
void AliAnalysisTaskPMD::UserCreateOutputObjects()
{
  // Create histograms
  // Called once
  
  fOutputList = new TList();
  Int_t kNbinsMultA = 50;  Float_t XminMultA = -0.5; Float_t XmaxMultA = 49.5;
  Int_t kNbinsXY = 200; Float_t XminXY = -100.0; Float_t XmaxXY  = 100.0;
  Int_t kNBinsEvent = 10; Float_t XminEvent = 0; Float_t XmaxEvent = 10;
  Int_t kNBinsEta = 10; Float_t XminEta = 2.1; Float_t XmaxEta = 4.1;
 
  fHistTotEvent = new TH1F("TotEvent","TotEvent",
			   kNBinsEvent,XminEvent,XmaxEvent); 
  fOutputList->Add(fHistTotEvent);
  fHistTotEventAfterPhySel = new TH1F("TotEventAfterPhySel","TotEventAfterPhySel",
				      kNBinsEvent,XminEvent,XmaxEvent); 
  fOutputList->Add(fHistTotEventAfterPhySel);
  fHistTotEventAfterVtx = new TH1F("TotEventAfterVtx","TotEventAfterVtx",
				   kNBinsEvent,XminEvent,XmaxEvent); 
  fOutputList->Add(fHistTotEventAfterVtx);
  fHistVtxZ = new TH1F("VtxZ","VtxZ",100,-50,50);
  fOutputList->Add(fHistVtxZ);
  fHistXYPre = new TH2F("XYPre","XYPre",kNbinsXY,XminXY,XmaxXY,
			kNbinsXY,XminXY,XmaxXY);
  fOutputList->Add(fHistXYPre);
  fHistEta = new TH1F ("Eta","Eta",kNBinsEta,XminEta,XmaxEta);
  fOutputList->Add(fHistEta);
  fHistEta1 = new TH1F ("Eta1","Eta1",kNBinsEta,XminEta,XmaxEta);
  fOutputList->Add(fHistEta1);
  fHistMultMeas = new TH1F("MultM","MultM",100,-0.5,99.5);
  fOutputList->Add(fHistMultMeas);
  fHistMultMeas1 = new TH1F("MultM1","MultM1",100,-0.5,99.5);
  fOutputList->Add(fHistMultMeas1);
  
  Char_t nameM[256], nameM1[256];
  for(Int_t i=0; i<10; i++){
    sprintf(nameM,"MultM_EtaBin%d",i+1);
    fHistMultMeasEtaBinA[i] =
      new TH1F(nameM,nameM,kNbinsMultA,
	       XminMultA,XmaxMultA);
    fOutputList->Add(fHistMultMeasEtaBinA[i]);
    sprintf(nameM1,"MultM1_EtaBin%d",i+1);
    fHistMultMeasEtaBinA1[i] =
      new TH1F(nameM1,nameM1,kNbinsMultA,
	       XminMultA,XmaxMultA);
    fOutputList->Add(fHistMultMeasEtaBinA1[i]);
  }//i loop
  
}

//________________________________________________________________________
void AliAnalysisTaskPMD::UserExec(Option_t *) 
{
  // Main loop
  // Called for each event
  Float_t MipCut1 = 432;//MPV=72=> 6*MPV=432
  Float_t MipCut2 = 648;//MPV=72=> 9*MPV=648
  Float_t etacls=0., theta=0., rdist=0.;
  Int_t PhotonCls = 0;//# of photon measured within 2.3 to 3.9 #eta
  Int_t PhotonCls1 = 0;//# of photon measured within 2.3 to 3.9 #eta
  Int_t PhotonClsAEtaBin[10] = {0};//# of photon measured with diff Eta bin
  Int_t PhotonClsAEtaBin1[10] = {0};//# of photon measured with diff Eta bin
  // Post output data.
  fESD = dynamic_cast<AliESDEvent*>(InputEvent());
  if (!fESD) {
    printf("ERROR: fESD not available\n");
    return;
  }
  fHistTotEvent->Fill(5);
  //event selection
  Bool_t isSelected = (((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->IsEventSelected() & AliVEvent::kMB);
  if (! isSelected) return;
  fHistTotEventAfterPhySel->Fill(5);
  //printf("There are %d tracks in this event\n", fESD->GetNumberOfPmdTracks());
  //vertex selection
  const AliESDVertex *vertex = fESD->GetPrimaryVertex();
  Float_t Vz = vertex->GetZ();    
  Bool_t zVerStatus = vertex->GetStatus();
  if(zVerStatus){
    fHistVtxZ->Fill(Vz);
    if(TMath::Abs(Vz)<10){
      fHistTotEventAfterVtx->Fill(5);
      Int_t ptracks = fESD->GetNumberOfPmdTracks();
      //PMDtrack Loop
      for(Int_t kk=0;kk<ptracks;kk++){
	AliESDPmdTrack *pmdtr = fESD->GetPmdTrack(kk);
	Int_t   det   = pmdtr->GetDetector();
	Float_t clsX  = pmdtr->GetClusterX();
	Float_t clsY  = pmdtr->GetClusterY();
	Float_t clsZ  = pmdtr->GetClusterZ();
	clsZ -= Vz;
	Float_t ncell = pmdtr->GetClusterCells();
	Float_t adc   = pmdtr->GetClusterADC();
	//Int_t smn = pmdtr->GetSmn();
	//Calculation of #eta
	rdist = TMath::Sqrt(clsX*clsX + clsY*clsY);
	if(clsZ!=0) theta = TMath::ATan2(rdist,clsZ);
	etacls  = -TMath::Log(TMath::Tan(0.5*theta));
	
	if(det==0 && adc>72)fHistXYPre->Fill(clsX,clsY);
	if(det==0 && adc>MipCut1 && ncell>2){
	  if(etacls > 2.1 && etacls<= 2.3)PhotonClsAEtaBin[0]++;
	  if(etacls > 2.3 && etacls<= 2.5)PhotonClsAEtaBin[1]++;
	  if(etacls > 2.5 && etacls<= 2.7)PhotonClsAEtaBin[2]++;
	  if(etacls > 2.7 && etacls<= 2.9)PhotonClsAEtaBin[3]++;
	  if(etacls > 2.9 && etacls<= 3.1)PhotonClsAEtaBin[4]++;
	  if(etacls > 3.1 && etacls<= 3.3)PhotonClsAEtaBin[5]++;
	  if(etacls > 3.3 && etacls<= 3.5)PhotonClsAEtaBin[6]++;
	  if(etacls > 3.5 && etacls<= 3.7)PhotonClsAEtaBin[7]++;
	  if(etacls > 3.7 && etacls<= 3.9)PhotonClsAEtaBin[8]++;
	  if(etacls > 3.9 && etacls<= 4.1)PhotonClsAEtaBin[9]++;
	  if(etacls > 2.3 && etacls<= 3.9)PhotonCls++;
	  fHistEta->Fill(etacls);
	}//if MipCut1
	if(det==0 && adc>MipCut2 && ncell>2){
	  if(etacls > 2.1 && etacls<= 2.3)PhotonClsAEtaBin1[0]++;
	  if(etacls > 2.3 && etacls<= 2.5)PhotonClsAEtaBin1[1]++;
	  if(etacls > 2.5 && etacls<= 2.7)PhotonClsAEtaBin1[2]++;
	  if(etacls > 2.7 && etacls<= 2.9)PhotonClsAEtaBin1[3]++;
	  if(etacls > 2.9 && etacls<= 3.1)PhotonClsAEtaBin1[4]++;
	  if(etacls > 3.1 && etacls<= 3.3)PhotonClsAEtaBin1[5]++;
	  if(etacls > 3.3 && etacls<= 3.5)PhotonClsAEtaBin1[6]++;
	  if(etacls > 3.5 && etacls<= 3.7)PhotonClsAEtaBin1[7]++;
	  if(etacls > 3.7 && etacls<= 3.9)PhotonClsAEtaBin1[8]++;
	  if(etacls > 3.9 && etacls<= 4.1)PhotonClsAEtaBin1[9]++;
	  if(etacls > 2.3 && etacls<= 3.9)PhotonCls1++;
	  fHistEta1->Fill(etacls);
	}//if Mipcut2
      } //PMDtrack loop 
      for(Int_t i=0; i<10; i++){
	fHistMultMeasEtaBinA[i]->Fill(PhotonClsAEtaBin[i]);
	fHistMultMeasEtaBinA1[i]->Fill(PhotonClsAEtaBin1[i]);
      }//i loop
      fHistMultMeas->Fill(PhotonCls);
      fHistMultMeas1->Fill(PhotonCls1);
    }// Vz loop 
  }// Vzcut loop
  PostData(1, fOutputList);
}      

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

  fOutputList = dynamic_cast<TList*> (GetOutputData(1));
  if (!fOutputList) {
    printf("ERROR: Output list not available\n");
    return;
  }
}
Commit	Line	Data
d4806ef4	1	#include "TChain.h"
	2	#include "TTree.h"
	3	#include "TH1F.h"
	4	#include "TH2F.h"
	5	#include "TCanvas.h"
	6	#include "AliAnalysisTask.h"
	7	#include "AliAnalysisManager.h"
	8	#include "AliESDEvent.h"
	9	#include "AliESDInputHandler.h"
	10	#include "AliESDPmdTrack.h"
	11	#include "AliESDVertex.h"
	12	#include "AliESDVZERO.h"
	13	#include "AliAnalysisTaskPMD.h"
	14
	15	// AnalysisTask For PMD
	16	// Authors: Sudipan De, Subhash Singha
	17
	18	ClassImp(AliAnalysisTaskPMD)
	19
	20	//________________________________________________________________________
	21	AliAnalysisTaskPMD::AliAnalysisTaskPMD(const char *name)
	22	: AliAnalysisTaskSE(name),
	23	fESD(0),
	24	fOutputList(0),
	25	fHistTotEvent(0),
	26	fHistTotEventAfterPhySel(0),
	27	fHistTotEventAfterVtx(0),
	28	fHistVtxZ(0),
	29	fHistXYPre(0),
	30	fHistEta(0),
	31	fHistEta1(0),
	32	fHistMultMeas(0),
	33	fHistMultMeas1(0)
	34	{
	35	for(Int_t i=0; i<10; i++){
	36	fHistMultMeasEtaBinA[i] = 0;
	37	fHistMultMeasEtaBinA1[i] = 0;
	38	}
	39	// Constructor
	40
	41	// Define input and output slots here
	42	// Input slot #0 works with a TChain
	43	DefineInput(0, TChain::Class());
	44	// Output slot #0 id reserved by the base class for AOD
	45	// Output slot #1 writes into a TH1 container
	46	DefineOutput(1, TList::Class());
	47	}
	48
	49	//________________________________________________________________________
	50	void AliAnalysisTaskPMD::UserCreateOutputObjects()
	51	{
	52	// Create histograms
	53	// Called once
	54
	55	fOutputList = new TList();
	56	Int_t kNbinsMultA = 50; Float_t XminMultA = -0.5; Float_t XmaxMultA = 49.5;
	57	Int_t kNbinsXY = 200; Float_t XminXY = -100.0; Float_t XmaxXY = 100.0;
	58	Int_t kNBinsEvent = 10; Float_t XminEvent = 0; Float_t XmaxEvent = 10;
	59	Int_t kNBinsEta = 10; Float_t XminEta = 2.1; Float_t XmaxEta = 4.1;
	60
	61	fHistTotEvent = new TH1F("TotEvent","TotEvent",
	62	kNBinsEvent,XminEvent,XmaxEvent);
	63	fOutputList->Add(fHistTotEvent);
	64	fHistTotEventAfterPhySel = new TH1F("TotEventAfterPhySel","TotEventAfterPhySel",
65	kNBinsEvent,XminEvent,XmaxEvent);
66	fOutputList->Add(fHistTotEventAfterPhySel);
67	fHistTotEventAfterVtx = new TH1F("TotEventAfterVtx","TotEventAfterVtx",
68	kNBinsEvent,XminEvent,XmaxEvent);
69	fOutputList->Add(fHistTotEventAfterVtx);
70	fHistVtxZ = new TH1F("VtxZ","VtxZ",100,-50,50);
71	fOutputList->Add(fHistVtxZ);
72	fHistXYPre = new TH2F("XYPre","XYPre",kNbinsXY,XminXY,XmaxXY,
73	kNbinsXY,XminXY,XmaxXY);
74	fOutputList->Add(fHistXYPre);
75	fHistEta = new TH1F ("Eta","Eta",kNBinsEta,XminEta,XmaxEta);
76	fOutputList->Add(fHistEta);
77	fHistEta1 = new TH1F ("Eta1","Eta1",kNBinsEta,XminEta,XmaxEta);
78	fOutputList->Add(fHistEta1);
79	fHistMultMeas = new TH1F("MultM","MultM",100,-0.5,99.5);
80	fOutputList->Add(fHistMultMeas);
81	fHistMultMeas1 = new TH1F("MultM1","MultM1",100,-0.5,99.5);
82	fOutputList->Add(fHistMultMeas1);
83
84	Char_t nameM[256], nameM1[256];
85	for(Int_t i=0; i<10; i++){
86	sprintf(nameM,"MultM_EtaBin%d",i+1);
87	fHistMultMeasEtaBinA[i] =
88	new TH1F(nameM,nameM,kNbinsMultA,
89	XminMultA,XmaxMultA);
90	fOutputList->Add(fHistMultMeasEtaBinA[i]);
91	sprintf(nameM1,"MultM1_EtaBin%d",i+1);
92	fHistMultMeasEtaBinA1[i] =
93	new TH1F(nameM1,nameM1,kNbinsMultA,
94	XminMultA,XmaxMultA);
95	fOutputList->Add(fHistMultMeasEtaBinA1[i]);
96	}//i loop
97
98	}
99
100	//________________________________________________________________________
101	void AliAnalysisTaskPMD::UserExec(Option_t *)
102	{
103	// Main loop
104	// Called for each event
105	Float_t MipCut1 = 432;//MPV=72=> 6*MPV=432
106	Float_t MipCut2 = 648;//MPV=72=> 9*MPV=648
759ae157	107	Float_t etacls=0., theta=0., rdist=0.;
d4806ef4	108	Int_t PhotonCls = 0;//# of photon measured within 2.3 to 3.9 #eta
	109	Int_t PhotonCls1 = 0;//# of photon measured within 2.3 to 3.9 #eta
	110	Int_t PhotonClsAEtaBin[10] = {0};//# of photon measured with diff Eta bin
	111	Int_t PhotonClsAEtaBin1[10] = {0};//# of photon measured with diff Eta bin
	112	// Post output data.
	113	fESD = dynamic_cast<AliESDEvent*>(InputEvent());
	114	if (!fESD) {
	115	printf("ERROR: fESD not available\n");
	116	return;
	117	}
	118	fHistTotEvent->Fill(5);
	119	//event selection
	120	Bool_t isSelected = (((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->IsEventSelected() & AliVEvent::kMB);
	121	if (! isSelected) return;
	122	fHistTotEventAfterPhySel->Fill(5);
	123	//printf("There are %d tracks in this event\n", fESD->GetNumberOfPmdTracks());
	124	//vertex selection
	125	const AliESDVertex *vertex = fESD->GetPrimaryVertex();
e690d4d0	126	Float_t Vz = vertex->GetZ();
d4806ef4	127	Bool_t zVerStatus = vertex->GetStatus();
	128	if(zVerStatus){
	129	fHistVtxZ->Fill(Vz);
	130	if(TMath::Abs(Vz)<10){
	131	fHistTotEventAfterVtx->Fill(5);
	132	Int_t ptracks = fESD->GetNumberOfPmdTracks();
	133	//PMDtrack Loop
	134	for(Int_t kk=0;kk<ptracks;kk++){
	135	AliESDPmdTrack *pmdtr = fESD->GetPmdTrack(kk);
	136	Int_t det = pmdtr->GetDetector();
	137	Float_t clsX = pmdtr->GetClusterX();
	138	Float_t clsY = pmdtr->GetClusterY();
	139	Float_t clsZ = pmdtr->GetClusterZ();
	140	clsZ -= Vz;
	141	Float_t ncell = pmdtr->GetClusterCells();
	142	Float_t adc = pmdtr->GetClusterADC();
	143	//Int_t smn = pmdtr->GetSmn();
	144	//Calculation of #eta
	145	rdist = TMath::Sqrt(clsXclsX + clsYclsY);
	146	if(clsZ!=0) theta = TMath::ATan2(rdist,clsZ);
	147	etacls = -TMath::Log(TMath::Tan(0.5*theta));
	148
	149	if(det==0 && adc>72)fHistXYPre->Fill(clsX,clsY);
	150	if(det==0 && adc>MipCut1 && ncell>2){
	151	if(etacls > 2.1 && etacls<= 2.3)PhotonClsAEtaBin[0]++;
	152	if(etacls > 2.3 && etacls<= 2.5)PhotonClsAEtaBin[1]++;
	153	if(etacls > 2.5 && etacls<= 2.7)PhotonClsAEtaBin[2]++;
	154	if(etacls > 2.7 && etacls<= 2.9)PhotonClsAEtaBin[3]++;
	155	if(etacls > 2.9 && etacls<= 3.1)PhotonClsAEtaBin[4]++;
	156	if(etacls > 3.1 && etacls<= 3.3)PhotonClsAEtaBin[5]++;
	157	if(etacls > 3.3 && etacls<= 3.5)PhotonClsAEtaBin[6]++;
	158	if(etacls > 3.5 && etacls<= 3.7)PhotonClsAEtaBin[7]++;
	159	if(etacls > 3.7 && etacls<= 3.9)PhotonClsAEtaBin[8]++;
	160	if(etacls > 3.9 && etacls<= 4.1)PhotonClsAEtaBin[9]++;
	161	if(etacls > 2.3 && etacls<= 3.9)PhotonCls++;
	162	fHistEta->Fill(etacls);
	163	}//if MipCut1
	164	if(det==0 && adc>MipCut2 && ncell>2){
	165	if(etacls > 2.1 && etacls<= 2.3)PhotonClsAEtaBin1[0]++;
	166	if(etacls > 2.3 && etacls<= 2.5)PhotonClsAEtaBin1[1]++;
	167	if(etacls > 2.5 && etacls<= 2.7)PhotonClsAEtaBin1[2]++;
	168	if(etacls > 2.7 && etacls<= 2.9)PhotonClsAEtaBin1[3]++;
	169	if(etacls > 2.9 && etacls<= 3.1)PhotonClsAEtaBin1[4]++;
	170	if(etacls > 3.1 && etacls<= 3.3)PhotonClsAEtaBin1[5]++;
	171	if(etacls > 3.3 && etacls<= 3.5)PhotonClsAEtaBin1[6]++;
	172	if(etacls > 3.5 && etacls<= 3.7)PhotonClsAEtaBin1[7]++;
	173	if(etacls > 3.7 && etacls<= 3.9)PhotonClsAEtaBin1[8]++;
	174	if(etacls > 3.9 && etacls<= 4.1)PhotonClsAEtaBin1[9]++;
	175	if(etacls > 2.3 && etacls<= 3.9)PhotonCls1++;
	176	fHistEta1->Fill(etacls);
	177	}//if Mipcut2
	178	} //PMDtrack loop
	179	for(Int_t i=0; i<10; i++){
	180	fHistMultMeasEtaBinA[i]->Fill(PhotonClsAEtaBin[i]);
	181	fHistMultMeasEtaBinA1[i]->Fill(PhotonClsAEtaBin1[i]);
	182	}//i loop
	183	fHistMultMeas->Fill(PhotonCls);
	184	fHistMultMeas1->Fill(PhotonCls1);
	185	}// Vz loop
	186	}// Vzcut loop
	187	PostData(1, fOutputList);
	188	}
	189
	190	//________________________________________________________________________
191	void AliAnalysisTaskPMD::Terminate(Option_t *)
192	{
193	// Draw result to the screen
194	// Called once at the end of the query
195
196	fOutputList = dynamic_cast<TList*> (GetOutputData(1));
197	if (!fOutputList) {
198	printf("ERROR: Output list not available\n");
199	return;
200	}
201	}