[u/mrichter/AliRoot.git] / PMD / anal / AliPMDAnalysisTaskPbPb.cxx

/**************************************************************************
 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
 *                                                                        *
 * Author: Satyajit Jena.                                                 *
 * 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.                  *
 **************************************************************************/
/*
  .
  .                  A template class to read tracks (PMD Cluster)
  .		         Runs in Local and Grid Modes
  .			Can be used for PbPb PMD analysis
  .		     Auther: Satyajit Jena <sjena@cern.ch>
  .                               12/04/2012  
  .
  . You are free to ask me anything you want about this code, I will add
  . a readme shortly to this folder about the running and adjusting the
  . code   
**************************************************************************/


#include "TChain.h"
#include "TTree.h"
#include "TH1F.h"
#include "TCanvas.h"

#include "AliAnalysisTask.h"
#include "AliAnalysisTaskSE.h"
#include "AliAnalysisManager.h"

#include "AliESDPmdTrack.h"
#include "AliESDEvent.h"
#include "AliESDInputHandler.h"
#include "AliESDtrackCuts.h"
#include "AliPMDAnalysisTaskPbPb.h"
#include "AliCentrality.h"


ClassImp(AliPMDAnalysisTaskPbPb)

//________________________________________________________________________
AliPMDAnalysisTaskPbPb::AliPMDAnalysisTaskPbPb(const char *name) 
: AliAnalysisTaskSE(name), fOutputList(0), fTrackCuts(0), fESD(0), fHistPt(0),fHistEta(0),
  fhEsdXYP(0),
  fhEsdXYC(0),
  fIsMC(kFALSE){

  DefineOutput(1, TList::Class());
}

//________________________________________________________________________
void AliPMDAnalysisTaskPbPb::CreateOutputObjects()
{
  fOutputList = new TList();
  fOutputList->SetOwner();

  //  fTrackCuts = AliESDtrackCuts::GetStandardITSTPCTrackCuts2010(kTRUE);

  // Track Cut Setup
  // Tune it accoringly
  fTrackCuts = new  AliESDtrackCuts("MyTrack Cuts","tracks");
  
  fTrackCuts->SetMinNClustersTPC(70);
  fTrackCuts->SetMaxChi2PerClusterTPC(4);
  fTrackCuts->SetAcceptKinkDaughters(kFALSE);
  fTrackCuts->SetRequireTPCRefit(kTRUE);
  fTrackCuts->SetRequireITSRefit(kTRUE);
  fTrackCuts->SetClusterRequirementITS(AliESDtrackCuts::kSPD,AliESDtrackCuts::kAny);
  fTrackCuts->SetMaxDCAToVertexZ(2);

  fTrackCuts->SetDCAToVertex2D(kFALSE);
  fTrackCuts->SetRequireSigmaToVertex(kFALSE);
    
  // Create the histograms as you like. 
  Int_t etabins = 40;
  Float_t etalow = -2.0, etaup = 2.0;
  fHistEta = new TH1F("fHistEta","#eta distribution for reconstructed",etabins, etalow, etaup);
  fHistEta->GetXaxis()->SetTitle("#eta");
  fHistEta->GetYaxis()->SetTitle("counts");
  fHistEta->SetMarkerStyle(kFullCircle);
  fHistEta->SetMarkerColor(kRed);
  fOutputList->Add(fHistEta);

  fHistPt = new TH1F("fHistPt", "P_{T} distribution", 100, 0.1, 10.1);
  fHistPt->GetXaxis()->SetTitle("P_{T} (GeV/c)");
  fHistPt->GetYaxis()->SetTitle("dN/dP_{T} (c/GeV)");
  fHistPt->SetMarkerStyle(kFullCircle);
  fHistPt->SetMarkerColor(kRed);
  fOutputList->Add(fHistPt);
  
  fhEsdXYC = new TH2F("hEsdXYC"," Scattered Plot for CPV (ESD) ",2000,-100.,100.,2000,-100.,100.);
  fhEsdXYC->GetXaxis()->SetTitle("X-axis");
  fhEsdXYC->GetYaxis()->SetTitle("Y-axis");
  fOutputList->Add(fhEsdXYC);
  
  fhEsdXYP = new TH2F("hEsdXYP"," Scattered Plot for PRE (ESD) ",2000,-100.,100.,2000,-100.,100.);
  fhEsdXYP->GetXaxis()->SetTitle("X-axis");
  fhEsdXYP->GetYaxis()->SetTitle("Y-axis");
  fOutputList->Add(fhEsdXYP);


  PostData(1, fOutputList);
}

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

  fESD = dynamic_cast<AliESDEvent*>(InputEvent());
  if (!fESD) {
    Printf("ERROR: fESD not available");
    return;
  }
  
  Bool_t isSelected = (((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->IsEventSelected() & AliVEvent::kMB);
  if(!isSelected) return;

   
  if(!(fESD->GetPrimaryVertex()->GetStatus())) return;
  // if vertex is from spd vertexZ, require more stringent cut
  if (fESD->GetPrimaryVertex()->IsFromVertexerZ()) {
    if (fESD->GetPrimaryVertex()->GetDispersion()>0.02 ||  fESD->GetPrimaryVertex()->GetZRes()>0.25 ) return; // bad vertex from VertexerZ
  }
  

 AliCentrality *centrality = fESD->GetCentrality();

 /* If you want to make centrality bining use following */
 // Double_t nCentrality = -1;
 // nCentrality = centrality->GetCentralityPercentile("V0M");
 // printf(" %f \n", nCentrality);
 // if (nCentrality < 0) return;

 /* otherwise 
    if you want to cuttoff other centrality and keep a window 
    for example this is for 0 - 10% centrality with V0M estimator
*/
 if(!centrality->IsEventInCentralityClass(0,10,"V0M")) return;  
  
  
  // Track loop for reconstructed event - General purpose
  // any central analysis can be done.

  Int_t   nGoodTracks = 0;
  Int_t ntracks = fESD->GetNumberOfTracks();
  for(Int_t i = 0; i < ntracks; i++) {
    AliESDtrack* esdtrack = fESD->GetTrack(i); // pointer to reconstructed to track          
    if(!esdtrack) {
      AliError(Form("ERROR: Could not retrieve esdtrack %d",i));
      continue;
    }
    
    if(!fTrackCuts->AcceptTrack(esdtrack)) continue;
    nGoodTracks++;
    fHistPt->Fill(esdtrack->Pt());
    fHistEta->Fill(esdtrack->Eta());
  } 
  
   
  //_______________________________________________________
  // Only for PMD part
  
  Int_t npmdcl = fESD->GetNumberOfPmdTracks();
  printf ("Number of PMD Clusters: %8d \n",npmdcl) ;  

  for (Int_t i = 0; i < npmdcl; i++) {
    AliESDPmdTrack *pmdtr = fESD->GetPmdTrack(i);
    Int_t   det     = pmdtr->GetDetector();
    Float_t clsX    = pmdtr->GetClusterX();
    Float_t clsY    = pmdtr->GetClusterY();
    
    if ( det == 0 ) {
      fhEsdXYP->Fill(clsX,clsY);
    }
    else if ( det == 1 ) {
      fhEsdXYC->Fill(clsX,clsY);
    }
  }


  PostData(1, fOutputList);
}      

//________________________________________________________________________
void AliPMDAnalysisTaskPbPb::Terminate(Option_t *) {
  // Called once at the end of the query

  fOutputList = dynamic_cast<TList*> (GetOutputData(1));
  if (!fOutputList) {
    Printf("ERROR: fHistPt not available");
    return;
  }

  fHistPt = dynamic_cast<TH1F*>(fOutputList->FindObject("fHistPt"));
 TCanvas *c1 = new TCanvas("AliAnalysisTaskPt","Pt",10,10,745,412);
 c1->SetFillColor(10);
 c1->SetBorderMode(0);
 c1->SetBorderSize(0);
 c1->SetFrameBorderMode(0);
 c1->cd(1)->SetLogy();
 fHistPt->DrawCopy("E");


  fHistEta = dynamic_cast<TH1F*>(fOutputList->FindObject("fHistEta"));
  TCanvas *c2 = new TCanvas("eta","Pt",10,10,745,412);
  c2->SetFillColor(10);
  c2->SetBorderMode(0);
  c2->SetBorderSize(0);
  c2->SetFrameBorderMode(0);
  c2->cd();
  fHistEta->DrawCopy("E");


  fhEsdXYP = dynamic_cast<TH2F*>(fOutputList->FindObject("hEsdXYP"));
  TCanvas *c3 = new TCanvas("XYP","Pt",10,10,745,745);
  c3->SetFillColor(10);
  c3->SetBorderMode(0);
  c3->SetBorderSize(0);
  c3->SetFrameBorderMode(0);
  c3->cd();
  fhEsdXYP->DrawCopy();


fhEsdXYC = dynamic_cast<TH2F*>(fOutputList->FindObject("hEsdXYC"));
  TCanvas *c4 = new TCanvas("XYC","Pt",10,10,745,745);
  c4->SetFillColor(10);
  c4->SetBorderMode(0);
  c4->SetBorderSize(0);
  c4->SetFrameBorderMode(0);
  c4->cd();
  fhEsdXYC->DrawCopy();

 
}
Commit	Line	Data
1bb16ba5	1	/**************************************************************************
	2	* Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
	3	* *
	4	* Author: Satyajit Jena. *
	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	**************************************************************************/
	15	/*
	16	.
	17	. A template class to read tracks (PMD Cluster)
	18	. Runs in Local and Grid Modes
	19	. Can be used for PbPb PMD analysis
	20	. Auther: Satyajit Jena <sjena@cern.ch>
	21	. 12/04/2012
	22	.
	23	. You are free to ask me anything you want about this code, I will add
	24	. a readme shortly to this folder about the running and adjusting the
	25	. code
	26	**************************************************************************/
	27
	28
	29	#include "TChain.h"
	30	#include "TTree.h"
	31	#include "TH1F.h"
	32	#include "TCanvas.h"
	33
	34	#include "AliAnalysisTask.h"
	35	#include "AliAnalysisTaskSE.h"
	36	#include "AliAnalysisManager.h"
	37
	38	#include "AliESDPmdTrack.h"
	39	#include "AliESDEvent.h"
	40	#include "AliESDInputHandler.h"
	41	#include "AliESDtrackCuts.h"
	42	#include "AliPMDAnalysisTaskPbPb.h"
	43	#include "AliCentrality.h"
	44
	45
	46	ClassImp(AliPMDAnalysisTaskPbPb)
	47
	48	//________________________________________________________________________
	49	AliPMDAnalysisTaskPbPb::AliPMDAnalysisTaskPbPb(const char *name)
	50	: AliAnalysisTaskSE(name), fOutputList(0), fTrackCuts(0), fESD(0), fHistPt(0),fHistEta(0),
	51	fhEsdXYP(0),
	52	fhEsdXYC(0),
	53	fIsMC(kFALSE){
	54
	55	DefineOutput(1, TList::Class());
	56	}
	57
	58	//________________________________________________________________________
	59	void AliPMDAnalysisTaskPbPb::CreateOutputObjects()
	60	{
	61	fOutputList = new TList();
	62	fOutputList->SetOwner();
	63
	64	// fTrackCuts = AliESDtrackCuts::GetStandardITSTPCTrackCuts2010(kTRUE);
65
66	// Track Cut Setup
67	// Tune it accoringly
68	fTrackCuts = new AliESDtrackCuts("MyTrack Cuts","tracks");
69
70	fTrackCuts->SetMinNClustersTPC(70);
71	fTrackCuts->SetMaxChi2PerClusterTPC(4);
72	fTrackCuts->SetAcceptKinkDaughters(kFALSE);
73	fTrackCuts->SetRequireTPCRefit(kTRUE);
74	fTrackCuts->SetRequireITSRefit(kTRUE);
75	fTrackCuts->SetClusterRequirementITS(AliESDtrackCuts::kSPD,AliESDtrackCuts::kAny);
76	fTrackCuts->SetMaxDCAToVertexZ(2);
77
78	fTrackCuts->SetDCAToVertex2D(kFALSE);
79	fTrackCuts->SetRequireSigmaToVertex(kFALSE);
80
81	// Create the histograms as you like.
82	Int_t etabins = 40;
83	Float_t etalow = -2.0, etaup = 2.0;
84	fHistEta = new TH1F("fHistEta","#eta distribution for reconstructed",etabins, etalow, etaup);
85	fHistEta->GetXaxis()->SetTitle("#eta");
86	fHistEta->GetYaxis()->SetTitle("counts");
87	fHistEta->SetMarkerStyle(kFullCircle);
88	fHistEta->SetMarkerColor(kRed);
89	fOutputList->Add(fHistEta);
90
91	fHistPt = new TH1F("fHistPt", "P_{T} distribution", 100, 0.1, 10.1);
92	fHistPt->GetXaxis()->SetTitle("P_{T} (GeV/c)");
93	fHistPt->GetYaxis()->SetTitle("dN/dP_{T} (c/GeV)");
94	fHistPt->SetMarkerStyle(kFullCircle);
95	fHistPt->SetMarkerColor(kRed);
96	fOutputList->Add(fHistPt);
97
98	fhEsdXYC = new TH2F("hEsdXYC"," Scattered Plot for CPV (ESD) ",2000,-100.,100.,2000,-100.,100.);
99	fhEsdXYC->GetXaxis()->SetTitle("X-axis");
100	fhEsdXYC->GetYaxis()->SetTitle("Y-axis");
101	fOutputList->Add(fhEsdXYC);
102
103	fhEsdXYP = new TH2F("hEsdXYP"," Scattered Plot for PRE (ESD) ",2000,-100.,100.,2000,-100.,100.);
104	fhEsdXYP->GetXaxis()->SetTitle("X-axis");
105	fhEsdXYP->GetYaxis()->SetTitle("Y-axis");
106	fOutputList->Add(fhEsdXYP);
107
108
109	PostData(1, fOutputList);
110	}
111
112	//________________________________________________________________________
113	void AliPMDAnalysisTaskPbPb::Exec(Option_t *)
114	{
115	// Main loop
116	// Called for each event
117
118	fESD = dynamic_cast<AliESDEvent*>(InputEvent());
119	if (!fESD) {
120	Printf("ERROR: fESD not available");
121	return;
122	}
123
124	Bool_t isSelected = (((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->IsEventSelected() & AliVEvent::kMB);
125	if(!isSelected) return;
126
127
128	if(!(fESD->GetPrimaryVertex()->GetStatus())) return;
129	// if vertex is from spd vertexZ, require more stringent cut
130	if (fESD->GetPrimaryVertex()->IsFromVertexerZ()) {
131	if (fESD->GetPrimaryVertex()->GetDispersion()>0.02 \|\| fESD->GetPrimaryVertex()->GetZRes()>0.25 ) return; // bad vertex from VertexerZ
132	}
133
134
135	AliCentrality *centrality = fESD->GetCentrality();
136
137	/* If you want to make centrality bining use following */
138	// Double_t nCentrality = -1;
139	// nCentrality = centrality->GetCentralityPercentile("V0M");
140	// printf(" %f \n", nCentrality);
141	// if (nCentrality < 0) return;
142
143	/* otherwise
144	if you want to cuttoff other centrality and keep a window
145	for example this is for 0 - 10% centrality with V0M estimator
146	*/
147	if(!centrality->IsEventInCentralityClass(0,10,"V0M")) return;
148
149
150	// Track loop for reconstructed event - General purpose
151	// any central analysis can be done.
152
153	Int_t nGoodTracks = 0;
154	Int_t ntracks = fESD->GetNumberOfTracks();
155	for(Int_t i = 0; i < ntracks; i++) {
156	AliESDtrack* esdtrack = fESD->GetTrack(i); // pointer to reconstructed to track
157	if(!esdtrack) {
158	AliError(Form("ERROR: Could not retrieve esdtrack %d",i));
159	continue;
160	}
161
162	if(!fTrackCuts->AcceptTrack(esdtrack)) continue;
163	nGoodTracks++;
164	fHistPt->Fill(esdtrack->Pt());
165	fHistEta->Fill(esdtrack->Eta());
166	}
167
168
169	//_______________________________________________________
170	// Only for PMD part
171
172	Int_t npmdcl = fESD->GetNumberOfPmdTracks();
173	printf ("Number of PMD Clusters: %8d \n",npmdcl) ;
174
175	for (Int_t i = 0; i < npmdcl; i++) {
176	AliESDPmdTrack *pmdtr = fESD->GetPmdTrack(i);
177	Int_t det = pmdtr->GetDetector();
178	Float_t clsX = pmdtr->GetClusterX();
179	Float_t clsY = pmdtr->GetClusterY();
180
181	if ( det == 0 ) {
182	fhEsdXYP->Fill(clsX,clsY);
183	}
184	else if ( det == 1 ) {
185	fhEsdXYC->Fill(clsX,clsY);
186	}
187	}
188
189
190	PostData(1, fOutputList);
191	}
192
193	//________________________________________________________________________
194	void AliPMDAnalysisTaskPbPb::Terminate(Option_t *) {
195	// Called once at the end of the query
196
197	fOutputList = dynamic_cast<TList*> (GetOutputData(1));
198	if (!fOutputList) {
199	Printf("ERROR: fHistPt not available");
200	return;
201	}
202
203	fHistPt = dynamic_cast<TH1F*>(fOutputList->FindObject("fHistPt"));
204	TCanvas *c1 = new TCanvas("AliAnalysisTaskPt","Pt",10,10,745,412);
205	c1->SetFillColor(10);
206	c1->SetBorderMode(0);
207	c1->SetBorderSize(0);
208	c1->SetFrameBorderMode(0);
209	c1->cd(1)->SetLogy();
210	fHistPt->DrawCopy("E");
211
212
213	fHistEta = dynamic_cast<TH1F*>(fOutputList->FindObject("fHistEta"));
214	TCanvas *c2 = new TCanvas("eta","Pt",10,10,745,412);
215	c2->SetFillColor(10);
216	c2->SetBorderMode(0);
217	c2->SetBorderSize(0);
218	c2->SetFrameBorderMode(0);
219	c2->cd();
220	fHistEta->DrawCopy("E");
221
222
223	fhEsdXYP = dynamic_cast<TH2F*>(fOutputList->FindObject("hEsdXYP"));
224	TCanvas *c3 = new TCanvas("XYP","Pt",10,10,745,745);
225	c3->SetFillColor(10);
226	c3->SetBorderMode(0);
227	c3->SetBorderSize(0);
228	c3->SetFrameBorderMode(0);
229	c3->cd();
230	fhEsdXYP->DrawCopy();
231
232
233	fhEsdXYC = dynamic_cast<TH2F*>(fOutputList->FindObject("hEsdXYC"));
234	TCanvas *c4 = new TCanvas("XYC","Pt",10,10,745,745);
235	c4->SetFillColor(10);
236	c4->SetBorderMode(0);
237	c4->SetBorderSize(0);
238	c4->SetFrameBorderMode(0);
239	c4->cd();
240	fhEsdXYC->DrawCopy();
241
242
243
244
245	}
246