[u/mrichter/AliRoot.git] / PWGCF / EBYE / Fluctuations / AliEbyEFluctuationAnalysisTaskTrain.cxx

#include "TChain.h"\r
#include "TTree.h"\r
#include "TH1F.h"\r
#include "TH2F.h"\r
#include "TCanvas.h"\r
#include "TParticle.h"\r
#include "TObjArray.h"\r
\r
#include "AliAnalysisTask.h"\r
#include "AliAnalysisManager.h"\r
\r
#include "AliStack.h"\r
#include "AliMCEvent.h"\r
#include "AliESDEvent.h"\r
#include "AliESDInputHandler.h"\r
//#include "AliESDtrackCuts.h"\r
#include "AliCentrality.h"\r
\r
#include "AliEbyEFluctuationAnalysisTaskTrain.h"\r
\r
// Event by event charge fluctuation analysis\r
// Authors: Satyajit Jena and Panos Cristakoglou\r
// \r
\r
ClassImp(AliEbyEFluctuationAnalysisTaskTrain)\r
\r
//________________________________________________________________________\r
AliEbyEFluctuationAnalysisTaskTrain::AliEbyEFluctuationAnalysisTaskTrain(const char *name) \r
  : AliAnalysisTaskSE(name), fESD(0), fOutputList(0), \r
    fHistEventStats(0), fHistCentrality(0),\r
    fHistNMult(0), fHistNPlusNMinus(0),\r
    fHistNMultMC(0), fHistNPlusNMinusMC(0),\r
    fAnalysisType("ESD"), fAnalysisMode("TPC"),\r
    fCentralityEstimator("V0M"), \r
    fCentralityPercentileMin(0.), fCentralityPercentileMax(5.),\r
    fVxMax(3.), fVyMax(3.), fVzMax(10.),\r
    fMinNumberOfTPCClusters(80), fMaxChi2PerTPCCluster(4.0), \r
    fMaxDCAxy(3.0), fMaxDCAz(3.0) {\r
  // Constructor\r
  \r
  // Define input and output slots here\r
  // Input slot #0 works with a TChain\r
  DefineInput(0, TChain::Class());\r
  // Output slot #0 id reserved by the base class for AOD\r
  // Output slot #1 writes into a TH1 container\r
  DefineOutput(1, TList::Class());\r
}\r
\r
//________________________________________________________________________\r
void AliEbyEFluctuationAnalysisTaskTrain::UserCreateOutputObjects() {\r
  // Create histograms\r
  // Called once\r
\r
  fOutputList = new TList();\r
  fOutputList->SetOwner();\r
\r
  //Event stats.\r
  TString gCutName[4] = {"Total","Offline trigger",\r
                         "Vertex","Analyzed"};\r
  fHistEventStats = new TH1F("fHistEventStats",\r
                             "Event statistics;;N_{events}",\r
                             4,0.5,4.5);\r
  for(Int_t i = 1; i <= 4; i++)\r
    fHistEventStats->GetXaxis()->SetBinLabel(i,gCutName[i-1].Data());\r
  fOutputList->Add(fHistEventStats);\r
\r
  //ESD analysis\r
  if(fAnalysisType.CompareTo("ESD") == 0 ) {\r
    fHistCentrality = new TH1F("fHistCentrality",";Centrality bin;Events",\r
			       20,0.5,20.5);\r
    fOutputList->Add(fHistCentrality);\r
    \r
    TString histName;\r
    histName = "fHistNMult";\r
    fHistNMult = new TH1F(histName.Data(), \r
			  ";N_{mult.}",\r
			  800,0,4000);\r
    fOutputList->Add(fHistNMult);\r
\r
    histName = "fHistNPlusNMinus";\r
    fHistNPlusNMinus = new TH2F(histName.Data(), \r
				";N_{+};N_{-}",\r
				2000,0.5,2000.5,2000,0.5,2000.5);  \r
    fOutputList->Add(fHistNPlusNMinus);\r
  }//ESD analysis\r
  else if(fAnalysisType.CompareTo("MC") == 0) {\r
    TString histName = "fHistNMultMC";\r
    fHistNMultMC = new TH1F(histName.Data(), \r
			    ";N_{mult.}",\r
			    800,0,8000);\r
    fOutputList->Add(fHistNMultMC);\r
    \r
    histName = "fHistNPlusNMinusMC";\r
    fHistNPlusNMinusMC = new TH2F(histName.Data(), \r
				  ";N_{+};N_{-}",\r
				  3000,0.5,3000.5,3000,0.5,3000.5);  \r
    fOutputList->Add(fHistNPlusNMinusMC);\r
  }//MC analysis\r
\r
  PostData(1, fOutputList);\r
}\r
\r
//________________________________________________________________________\r
void AliEbyEFluctuationAnalysisTaskTrain::UserExec(Option_t *) {\r
  // Main loop\r
  // Called for each event\r
  \r
  Int_t nPlus = 0, nMinus = 0;\r
\r
  // Post output data.\r
  //ESD analysis\r
  if(fAnalysisType.CompareTo("ESD") == 0) {\r
    fESD = dynamic_cast<AliESDEvent*>(InputEvent());\r
    if (!fESD) {\r
      printf("ERROR: fESD not available\n");\r
      return;\r
    }\r
  \r
    fHistEventStats->Fill(1); //all events\r
    \r
    //Bool_t isSelected = (((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->IsEventSelected() & AliVEvent::kMB);\r
    //if(isSelected) {\r
      fHistEventStats->Fill(2); //triggered + centrality\r
      \r
\r
      //Centrality stuff\r
      AliCentrality *centrality = fESD->GetCentrality();\r
      Int_t nCentrality = 0;\r
      \r
      if(centrality->IsEventInCentralityClass(fCentralityPercentileMin,\r
					      fCentralityPercentileMax,\r
					      fCentralityEstimator.Data())) {\r
	if(fAnalysisMode.CompareTo("TPC") == 0 ) {\r
	  const AliESDVertex *vertex = fESD->GetPrimaryVertexTPC();\r
	  if(vertex) {\r
	    if(vertex->GetNContributors() > 0) {\r
	      if(vertex->GetZRes() != 0) {\r
		fHistEventStats->Fill(3); //events with a proper vertex\r
		if(TMath::Abs(vertex->GetXv()) < fVxMax) {\r
		  if(TMath::Abs(vertex->GetYv()) < fVyMax) {\r
		    if(TMath::Abs(vertex->GetZv()) < fVzMax) {\r
		      fHistEventStats->Fill(4); //analyzed events\r
		      \r
		      Printf("Centrality percentile: %lf - Centrality: %d - Total tracks: %d",\r
                             centrality->GetCentralityPercentile(fCentralityEstimator.Data()),\r
                             nCentrality,fESD->GetNumberOfTracks());\r
		      \r
\r
		      Int_t nAcceptedTracks = 0;\r
		      \r
		      nAcceptedTracks = fESD->GetNumberOfTracks();\r
		      AliESDtrack* track = 0;\r
		      \r
		      Float_t dcaXY = 0.0, dcaZ = 0.0;\r
		      \r
		      for (Int_t iTracks = 0; iTracks < nAcceptedTracks; iTracks++) \r
			{\r
			  track = dynamic_cast<AliESDtrack *>(fESD->GetTrack(iTracks));\r
			  if (!track) {\r
			    printf("ERROR: Could not receive track %d\n", iTracks);\r
			    continue;\r
			  }\r
			  \r
			  AliESDtrack *tpcOnlyTrack = new AliESDtrack();\r
			  \r
			  if (!track->FillTPCOnlyTrack(*tpcOnlyTrack)) {\r
			    delete tpcOnlyTrack;\r
			    continue;\r
			  }\r
			  \r
			  tpcOnlyTrack->GetImpactParameters(dcaXY,dcaZ);\r
			  			  \r
			  Int_t nClustersITS = track->GetITSclusters(0x0);\r
			  Int_t nClustersTPC = track->GetTPCclusters(0x0);\r
			  \r
			  Float_t chi2PerClusterITS = -1;\r
			  if (nClustersITS!=0)\r
			    chi2PerClusterITS = track->GetITSchi2()/Float_t(nClustersITS);\r
			  Float_t chi2PerClusterTPC = -1;\r
			  if (nClustersTPC!=0)\r
			    chi2PerClusterTPC = track->GetTPCchi2()/Float_t(nClustersTPC);\r
			  \r
			  //	\r
			  \r
			  if(nClustersTPC < fMinNumberOfTPCClusters) continue;\r
			  if(chi2PerClusterTPC > fMaxChi2PerTPCCluster) continue;\r
			  if(TMath::Abs(dcaXY) > fMaxDCAxy) continue;\r
			  if(TMath::Abs(dcaZ) > fMaxDCAz) continue;\r
			  \r
			  Short_t gCharge = tpcOnlyTrack->Charge();\r
			  if(gCharge > 0) nPlus += 1;\r
			  if(gCharge < 0) nMinus += 1;\r
			  delete tpcOnlyTrack;\r
			} // track loop\r
		      \r
		      if( nPlus == 0 && nMinus == 0) return;\r
\r
		      //  cout<<"=========NPlus: "<<nPlus<<"=========NMinus: "<<nMinus<<endl;\r
		      fHistCentrality->Fill(nCentrality);\r
		      fHistNPlusNMinus->Fill(nPlus,nMinus);\r
		      fHistNMult->Fill(nPlus+nMinus);\r
		      \r
		      \r
		    }//Vz cut\r
		  }//Vy cut\r
		}//Vx cut\r
	      }//Vz resolution\r
	    }//number of contributors\r
	  }//valid vertex\r
	}//TPC analysis mode\r
      }//centrality\r
      //}//physics selection\r
  }//ESD analysis level\r
  \r
      //MC analysis\r
  if(fAnalysisType.CompareTo("MC") == 0 ) {\r
    AliMCEvent* mcEvent = MCEvent();\r
    if (!mcEvent) {\r
      Printf("ERROR: Could not retrieve MC event");\r
      return;\r
    }\r
    AliStack* stack = mcEvent->Stack();\r
    if (!stack) {\r
      Printf("ERROR: Could not retrieve MC stack");\r
      return;\r
    }\r
    \r
    fHistEventStats->Fill(1); \r
    fHistEventStats->Fill(2); \r
    fHistEventStats->Fill(3); \r
    fHistEventStats->Fill(4); //analyzed events\r
    for (Int_t iTracks = 0; iTracks < mcEvent->GetNumberOfTracks(); iTracks++) {\r
      AliVParticle* track = mcEvent->GetTrack(iTracks);\r
      if (!track) {\r
	Printf("ERROR: Could not receive track %d (mc loop)", iTracks);\r
	continue;\r
      }\r
      \r
      if(!stack->IsPhysicalPrimary(iTracks)) continue;\r
      if((track->Pt() < 0.3) && (track->Pt() > 1.5)) continue;\r
      if(TMath::Abs(track->Eta()) > 0.8) continue;\r
      \r
      Short_t gCharge = track->Charge();\r
      if(gCharge > 0) nPlus += 1;\r
      if(gCharge < 0) nMinus += 1;\r
    }//particle loop\r
    fHistNPlusNMinusMC->Fill(nPlus,nMinus);\r
    fHistNMultMC->Fill(nPlus+nMinus);\r
    \r
  }//MC analysis level\r
\r
\r
}      \r
\r
//________________________________________________________________________\r
void AliEbyEFluctuationAnalysisTaskTrain::Terminate(Option_t *) {\r
  // Draw result to the screen\r
  // Called once at the end of the query\r
\r
  fOutputList = dynamic_cast<TList*> (GetOutputData(1));\r
  if (!fOutputList) {\r
    printf("ERROR: Output list not available\n");\r
    return;\r
  }\r
}\r
Commit	Line	Data
d43b78c9	1	#include "TChain.h"\r
	2	#include "TTree.h"\r
	3	#include "TH1F.h"\r
	4	#include "TH2F.h"\r
	5	#include "TCanvas.h"\r
	6	#include "TParticle.h"\r
	7	#include "TObjArray.h"\r
	8	\r
	9	#include "AliAnalysisTask.h"\r
	10	#include "AliAnalysisManager.h"\r
	11	\r
	12	#include "AliStack.h"\r
	13	#include "AliMCEvent.h"\r
	14	#include "AliESDEvent.h"\r
	15	#include "AliESDInputHandler.h"\r
	16	//#include "AliESDtrackCuts.h"\r
	17	#include "AliCentrality.h"\r
	18	\r
	19	#include "AliEbyEFluctuationAnalysisTaskTrain.h"\r
	20	\r
	21	// Event by event charge fluctuation analysis\r
	22	// Authors: Satyajit Jena and Panos Cristakoglou\r
	23	// \r
	24	\r
	25	ClassImp(AliEbyEFluctuationAnalysisTaskTrain)\r
	26	\r
	27	//________________________________________________________________________\r
	28	AliEbyEFluctuationAnalysisTaskTrain::AliEbyEFluctuationAnalysisTaskTrain(const char *name) \r
	29	: AliAnalysisTaskSE(name), fESD(0), fOutputList(0), \r
	30	fHistEventStats(0), fHistCentrality(0),\r
	31	fHistNMult(0), fHistNPlusNMinus(0),\r
	32	fHistNMultMC(0), fHistNPlusNMinusMC(0),\r
	33	fAnalysisType("ESD"), fAnalysisMode("TPC"),\r
	34	fCentralityEstimator("V0M"), \r
	35	fCentralityPercentileMin(0.), fCentralityPercentileMax(5.),\r
	36	fVxMax(3.), fVyMax(3.), fVzMax(10.),\r
	37	fMinNumberOfTPCClusters(80), fMaxChi2PerTPCCluster(4.0), \r
	38	fMaxDCAxy(3.0), fMaxDCAz(3.0) {\r
	39	// Constructor\r
	40	\r
	41	// Define input and output slots here\r
	42	// Input slot #0 works with a TChain\r
	43	DefineInput(0, TChain::Class());\r
	44	// Output slot #0 id reserved by the base class for AOD\r
	45	// Output slot #1 writes into a TH1 container\r
	46	DefineOutput(1, TList::Class());\r
	47	}\r
	48	\r
	49	//________________________________________________________________________\r
	50	void AliEbyEFluctuationAnalysisTaskTrain::UserCreateOutputObjects() {\r
	51	// Create histograms\r
	52	// Called once\r
	53	\r
	54	fOutputList = new TList();\r
	55	fOutputList->SetOwner();\r
	56	\r
	57	//Event stats.\r
	58	TString gCutName[4] = {"Total","Offline trigger",\r
	59	"Vertex","Analyzed"};\r
	60	fHistEventStats = new TH1F("fHistEventStats",\r
	61	"Event statistics;;N_{events}",\r
	62	4,0.5,4.5);\r
	63	for(Int_t i = 1; i <= 4; i++)\r
	64	fHistEventStats->GetXaxis()->SetBinLabel(i,gCutName[i-1].Data());\r
65	fOutputList->Add(fHistEventStats);\r
66	\r
67	//ESD analysis\r
68	if(fAnalysisType.CompareTo("ESD") == 0 ) {\r
69	fHistCentrality = new TH1F("fHistCentrality",";Centrality bin;Events",\r
70	20,0.5,20.5);\r
71	fOutputList->Add(fHistCentrality);\r
72	\r
73	TString histName;\r
74	histName = "fHistNMult";\r
75	fHistNMult = new TH1F(histName.Data(), \r
76	";N_{mult.}",\r
77	800,0,4000);\r
78	fOutputList->Add(fHistNMult);\r
79	\r
80	histName = "fHistNPlusNMinus";\r
81	fHistNPlusNMinus = new TH2F(histName.Data(), \r
82	";N_{+};N_{-}",\r
83	2000,0.5,2000.5,2000,0.5,2000.5); \r
84	fOutputList->Add(fHistNPlusNMinus);\r
85	}//ESD analysis\r
86	else if(fAnalysisType.CompareTo("MC") == 0) {\r
87	TString histName = "fHistNMultMC";\r
88	fHistNMultMC = new TH1F(histName.Data(), \r
89	";N_{mult.}",\r
90	800,0,8000);\r
91	fOutputList->Add(fHistNMultMC);\r
92	\r
93	histName = "fHistNPlusNMinusMC";\r
94	fHistNPlusNMinusMC = new TH2F(histName.Data(), \r
95	";N_{+};N_{-}",\r
96	3000,0.5,3000.5,3000,0.5,3000.5); \r
97	fOutputList->Add(fHistNPlusNMinusMC);\r
98	}//MC analysis\r
99	\r
100	PostData(1, fOutputList);\r
101	}\r
102	\r
103	//________________________________________________________________________\r
104	void AliEbyEFluctuationAnalysisTaskTrain::UserExec(Option_t *) {\r
105	// Main loop\r
106	// Called for each event\r
107	\r
108	Int_t nPlus = 0, nMinus = 0;\r
109	\r
110	// Post output data.\r
111	//ESD analysis\r
112	if(fAnalysisType.CompareTo("ESD") == 0) {\r
113	fESD = dynamic_cast<AliESDEvent*>(InputEvent());\r
114	if (!fESD) {\r
115	printf("ERROR: fESD not available\n");\r
116	return;\r
117	}\r
118	\r
119	fHistEventStats->Fill(1); //all events\r
120	\r
121	//Bool_t isSelected = (((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->IsEventSelected() & AliVEvent::kMB);\r
122	//if(isSelected) {\r
123	fHistEventStats->Fill(2); //triggered + centrality\r
124	\r
125	\r
126	//Centrality stuff\r
127	AliCentrality *centrality = fESD->GetCentrality();\r
128	Int_t nCentrality = 0;\r
129	\r
130	if(centrality->IsEventInCentralityClass(fCentralityPercentileMin,\r
131	fCentralityPercentileMax,\r
132	fCentralityEstimator.Data())) {\r
133	if(fAnalysisMode.CompareTo("TPC") == 0 ) {\r
134	const AliESDVertex *vertex = fESD->GetPrimaryVertexTPC();\r
135	if(vertex) {\r
136	if(vertex->GetNContributors() > 0) {\r
137	if(vertex->GetZRes() != 0) {\r
138	fHistEventStats->Fill(3); //events with a proper vertex\r
139	if(TMath::Abs(vertex->GetXv()) < fVxMax) {\r
140	if(TMath::Abs(vertex->GetYv()) < fVyMax) {\r
141	if(TMath::Abs(vertex->GetZv()) < fVzMax) {\r
142	fHistEventStats->Fill(4); //analyzed events\r
143	\r
144	Printf("Centrality percentile: %lf - Centrality: %d - Total tracks: %d",\r
145	centrality->GetCentralityPercentile(fCentralityEstimator.Data()),\r
146	nCentrality,fESD->GetNumberOfTracks());\r
147	\r
148	\r
149	Int_t nAcceptedTracks = 0;\r
150	\r
151	nAcceptedTracks = fESD->GetNumberOfTracks();\r
152	AliESDtrack* track = 0;\r
153	\r
154	Float_t dcaXY = 0.0, dcaZ = 0.0;\r
155	\r
156	for (Int_t iTracks = 0; iTracks < nAcceptedTracks; iTracks++) \r
157	{\r
158	track = dynamic_cast<AliESDtrack *>(fESD->GetTrack(iTracks));\r
159	if (!track) {\r
160	printf("ERROR: Could not receive track %d\n", iTracks);\r
161	continue;\r
162	}\r
163	\r
164	AliESDtrack *tpcOnlyTrack = new AliESDtrack();\r
165	\r
166	if (!track->FillTPCOnlyTrack(*tpcOnlyTrack)) {\r
167	delete tpcOnlyTrack;\r
168	continue;\r
169	}\r
170	\r
171	tpcOnlyTrack->GetImpactParameters(dcaXY,dcaZ);\r
172	\r
173	Int_t nClustersITS = track->GetITSclusters(0x0);\r
174	Int_t nClustersTPC = track->GetTPCclusters(0x0);\r
175	\r
176	Float_t chi2PerClusterITS = -1;\r
177	if (nClustersITS!=0)\r
178	chi2PerClusterITS = track->GetITSchi2()/Float_t(nClustersITS);\r
179	Float_t chi2PerClusterTPC = -1;\r
180	if (nClustersTPC!=0)\r
181	chi2PerClusterTPC = track->GetTPCchi2()/Float_t(nClustersTPC);\r
182	\r
183	// \r
184	\r
185	if(nClustersTPC < fMinNumberOfTPCClusters) continue;\r
186	if(chi2PerClusterTPC > fMaxChi2PerTPCCluster) continue;\r
187	if(TMath::Abs(dcaXY) > fMaxDCAxy) continue;\r
188	if(TMath::Abs(dcaZ) > fMaxDCAz) continue;\r
189	\r
190	Short_t gCharge = tpcOnlyTrack->Charge();\r
191	if(gCharge > 0) nPlus += 1;\r
192	if(gCharge < 0) nMinus += 1;\r
193	delete tpcOnlyTrack;\r
194	} // track loop\r
195	\r
196	if( nPlus == 0 && nMinus == 0) return;\r
197	\r
198	// cout<<"=========NPlus: "<<nPlus<<"=========NMinus: "<<nMinus<<endl;\r
199	fHistCentrality->Fill(nCentrality);\r
200	fHistNPlusNMinus->Fill(nPlus,nMinus);\r
201	fHistNMult->Fill(nPlus+nMinus);\r
202	\r
203	\r
204	}//Vz cut\r
205	}//Vy cut\r
206	}//Vx cut\r
207	}//Vz resolution\r
208	}//number of contributors\r
209	}//valid vertex\r
210	}//TPC analysis mode\r
211	}//centrality\r
212	//}//physics selection\r
213	}//ESD analysis level\r
214	\r
215	//MC analysis\r
216	if(fAnalysisType.CompareTo("MC") == 0 ) {\r
217	AliMCEvent* mcEvent = MCEvent();\r
218	if (!mcEvent) {\r
219	Printf("ERROR: Could not retrieve MC event");\r
220	return;\r
221	}\r
222	AliStack* stack = mcEvent->Stack();\r
223	if (!stack) {\r
224	Printf("ERROR: Could not retrieve MC stack");\r
225	return;\r
226	}\r
227	\r
228	fHistEventStats->Fill(1); \r
229	fHistEventStats->Fill(2); \r
230	fHistEventStats->Fill(3); \r
231	fHistEventStats->Fill(4); //analyzed events\r
232	for (Int_t iTracks = 0; iTracks < mcEvent->GetNumberOfTracks(); iTracks++) {\r
233	AliVParticle* track = mcEvent->GetTrack(iTracks);\r
234	if (!track) {\r
235	Printf("ERROR: Could not receive track %d (mc loop)", iTracks);\r
236	continue;\r
237	}\r
238	\r
239	if(!stack->IsPhysicalPrimary(iTracks)) continue;\r
240	if((track->Pt() < 0.3) && (track->Pt() > 1.5)) continue;\r
241	if(TMath::Abs(track->Eta()) > 0.8) continue;\r
242	\r
243	Short_t gCharge = track->Charge();\r
244	if(gCharge > 0) nPlus += 1;\r
245	if(gCharge < 0) nMinus += 1;\r
246	}//particle loop\r
247	fHistNPlusNMinusMC->Fill(nPlus,nMinus);\r
248	fHistNMultMC->Fill(nPlus+nMinus);\r
249	\r
250	}//MC analysis level\r
251	\r
252	\r
253	} \r
254	\r
255	//________________________________________________________________________\r
256	void AliEbyEFluctuationAnalysisTaskTrain::Terminate(Option_t *) {\r
257	// Draw result to the screen\r
258	// Called once at the end of the query\r
259	\r
260	fOutputList = dynamic_cast<TList*> (GetOutputData(1));\r
261	if (!fOutputList) {\r
262	printf("ERROR: Output list not available\n");\r
263	return;\r
264	}\r
265	}\r