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

/**************************************************************************
 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
 *                                                                        *
 * Author: The ALICE Off-line Project.                                    *
 * 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.                  *
 **************************************************************************/

/**************************************************************************
 * Analysis Class Implimentation for the MC Truth and responce Matrix 
 *    Auther:    Satyajit Jena, IIT Bombay |  sjena@cern.ch
 * 
 *                     Mon Nov 22 19:54:27 CET 2010
 * 
 * You may need some type of tunning to histograms and filling schemes,
 * Please contact me if you want something to change this code 
 **************************************************************************/

#include "TChain.h"
#include "TList.h"
#include "TFile.h"
#include "TTree.h"
#include "TH1F.h"
#include "TH1D.h"
#include "TH2F.h"
#include "TH3F.h"
#include "TCanvas.h"

#include "AliAnalysisTask.h"
#include "AliAnalysisManager.h"
#include "AliVEvent.h"
#include "AliAODEvent.h"

#include "AliMCEvent.h"
#include "AliStack.h"
#include "AliHeader.h"
#include "AliESD.h"
#include "AliESDEvent.h"

#include "AliMCEventHandler.h"
#include "AliGenEventHeader.h"
#include "AliGenPythiaEventHeader.h"
#include "AliGenDPMjetEventHeader.h"
#include "AliGenCocktailEventHeader.h"
#include "AliPID.h"

#include "AliESDInputHandler.h"
#include "AliPMDAnalysisMCTaskSE.h"

ClassImp(AliPMDAnalysisMCTaskSE)

//________________________________________________________________________
AliPMDAnalysisMCTaskSE::AliPMDAnalysisMCTaskSE(const char *name) 
  : AliAnalysisTaskSE(name),
  fPhysList(0),
//  fhCounter(0),
//  fhVtx(0),
//  fhVty(0),
//  fhVtz(0),
  fhResponseAll(0),
  fhTrueAll(0),
  fhMeasuredAll(0)
    
{
  for (Int_t i = 0; i < 10; i++) {
    fhResponse[i] = NULL;
    fhTrue[i] = NULL;
    fhMeasured[i] = NULL;
  }
  
  DefineOutput(1, TList::Class());
}


//________________________________________________________________________
void AliPMDAnalysisMCTaskSE::UserCreateOutputObjects()
{
  fPhysList   = new TList();
  fPhysList->SetOwner();
    
  //___________ Counter ________________
  // fCounter = new TH1D("hCounter","EVENT COUNTERS", 100, 0, 100);;
  // fPhysList->Add(fCounter);

  Int_t   kMultBin = 100;  
  Float_t gMultMin = -0.5; 
  Float_t gMultMax = 99.5;
  

  fhResponseAll = new TH2F("hResponseAll","Responce Matrix For All",
			   kMultBin,gMultMin,gMultMax,kMultBin,gMultMin,gMultMax);
  
  fhResponseAll->GetXaxis()->SetTitle("Measured Multiplicity #gamma_{like}");
  fhResponseAll->GetXaxis()->SetTitle("Ture Multiplicity #gamma");

  fhTrueAll = new TH2F("hTrueAll","True Multiplicity of #gamma",kMultBin,gMultMin,gMultMax);
  fhTrueAll->GetXaxis()->SetTitle("True Multiplicity #gamma");

  fhMeasuredAll = new TH2F("hMeasuredAll","Measured Multiplicity of #gamma",kMultBin,gMultMin,gMultMax);
  fhMeasuredAll->GetXaxis()->SetTitle("Measured Multiplicity #gamma");
 
  fPhysList->Add(fhResponseAll);
  fPhysList->Add(fhTrueAll);
  fPhysList->Add(fhMeasuredAll);
  
  for(Int_t i = 0; i < 10; i++) {
    Float_t gbin = 2.3 + i*0.2;
    
    char a[40];
    sprintf(a,"hResponseEtaBin%2.1fto%2.1f",gbin, gbin+0.2);
    char b[40];
    sprintf(b,"hTrueEtaBin%2.1fto%2.1f",gbin, gbin+0.2);
    char c[40];
    sprintf(c,"hMeasuredEtaBin%2.1fto%2.1f",gbin, gbin+0.2);

    char aa[40];
    sprintf(aa," Responce Matrix from \eta - %2.1f to %2.1f",gbin, gbin+0.2);
    char bb[40];
    sprintf(bb," True Multiplicity of #gamma from \eta %2.1f to %2.1f",gbin, gbin+0.2);
    char cc[40];
    sprintf(cc," Measured multiplicity of gamma [#gamma_{like}] from \eta %2.1f to %2.1f",gbin, gbin+0.2);
    

    fhResponse[i] = new TH2F(a,aa,kMultBin,gMultMin,gMultMax,kMultBin,gMultMin,gMultMax);
    fhResponse[i]->GetXaxis()->SetTitle("Measured Multiplicity #gamma_{like}");
    fhResponse[i]->GetYaxis()->SetTitle("Ture Multiplicity #gamma");


    fhTrue[i] = new TH1F(b,bb,kMultBin,gMultMin,gMultMax);
    fhResponse[i]->GetXaxis()->SetTitle("True Multiplicity #gamma");

    fhMeasured[i] = new TH1F(c,cc,kMultBin,gMultMin,gMultMax);
    fhMeasured[i]->GetXaxis()->SetTitle("True Multiplicity #gamma");

    fPhysList->Add(fhResponse[i]);
    fPhysList->Add(fhTrue[i]);
    fPhysList->Add(fhMeasured[i]);

  }

  PostData(1,fPhysList);
}

//________________________________________________________________________
void AliPMDAnalysisMCTaskSE::UserExec(Option_t *) 
{
  
  fCntr++;
  
  AliESDEvent* esd = dynamic_cast<AliESDEvent*>(InputEvent());
  AliStack *stack = 0;
  AliMCEventHandler *mcH = 0;
  AliMCEvent* mcEvent = 0;
  
  if(MCEvent()){
    stack = MCEvent()->Stack();
    mcH = (AliMCEventHandler*) ((AliAnalysisManager::GetAnalysisManager())->GetMCtruthEventHandler());
    mcEvent = mcH->MCEvent();
    AliLog::SetClassDebugLevel("AliMCEvent", AliLog::kWarning - AliLog::kDebug + 1);
  }
  if(!stack) return;
  if(!esd) return;
  if(!mcEvent)return;
  
  EventByEvent(esd,mcEvent);
  
}

//----------------------------
void AliPMDAnalysisMCTaskSE::EventByEvent(AliESDEvent* esd, AliMCEvent* mcEvent)
{ 
  Int_t trnn = mcEvent->GetNumberOfTracks();

  Float_t McEta = 0.;
 
  Int_t kcntr[0] ={0};
  Int_t kcntrr = 0;

  AliStack *stack = mcEvent->Stack();
  Int_t sln = 0;
  for (Int_t iTracks = 0; iTracks < trnn; iTracks++) {
    AliMCParticle* track = (AliMCParticle*)mcEvent->GetTrack(iTracks);
    // TParticle* part = track->Particle();

    TParticle* part = stack->Particle(iTracks);
    if (!part) continue;
    Int_t pid = part->GetPdgCode();
    TDatabasePDG* pdgDB = TDatabasePDG::Instance();
    TParticlePDG* pdgPart = pdgDB->GetParticle(pid);
    Float_t charge = (pdgPart ? pdgPart->Charge() : 0);
    
    McEta = part->Eta();

    if(McEta < 2.0 || McEta > 4.2) continue;

    kcntrr++;
    if (McEta >= 2.1 && McEta < 2.3 ) kcntr[0]++;
    else if (McEta >= 2.3 && McEta < 2.5 ) kcntr[1]++;
    else if (McEta >= 2.5 && McEta < 2.7 ) kcntr[2]++;
    else if (McEta >= 2.7 && McEta < 2.9 ) kcntr[3]++;
    else if (McEta >= 2.9 && McEta < 3.1 ) kcntr[4]++;
    else if (McEta >= 3.1 && McEta < 3.3 ) kcntr[5]++;
    else if (McEta >= 3.3 && McEta < 3.5 ) kcntr[6]++;
    else if (McEta >= 3.5 && McEta < 3.7 ) kcntr[7]++;
    else if (McEta >= 3.7 && McEta < 3.9 ) kcntr[8]++;  
    else if (McEta >= 3.9 && McEta < 4.1 ) kcntr[9]++;
    
  }
  
  // ESD
  
  Int_t npmdcl   = esd->GetNumberOfPmdTracks();
  Float_t rpxpy  = 0.;
  Float_t theta  = 0.; 
  Float_t pEta   = 0.;
 
  Int_t cntr[10] = {0};
  Int_t cntrr = 0;

  while (npmdcl--)
    {
      
      AliESDPmdTrack *pmdtr = esd->GetPmdTrack(npmdcl);
      Int_t det     = pmdtr->GetDetector();
      if( det == 1) continue;
      Int_t   ncell = pmdtr->GetClusterCells();
      if(ncell < 2) continue;
      Float_t adc   = pmdtr->GetClusterADC();
      if( adc < 216 ) continue;

      Float_t clsX  = pmdtr->GetClusterX();
      Float_t clsY  = pmdtr->GetClusterY();
      Float_t clsZ  = pmdtr->GetClusterZ();
      
      rpxpy  =  TMath::Sqrt(clsX*clsX + clsY*clsY);
      theta  =  TMath::ATan2(rpxpy,clsZ);
      pEta   = -TMath::Log(TMath::Tan(0.5*theta));
      cntrr++;
      if (pEta >= 2.1 && pEta < 2.3 ) cntr[0]++;
      else if (pEta >= 2.3 && pEta < 2.5 ) cntr[1]++;
      else if (pEta >= 2.5 && pEta < 2.7 ) cntr[2]++;
      else if (pEta >= 2.7 && pEta < 2.9 ) cntr[3]++;
      else if (pEta >= 2.9 && pEta < 3.1 ) cntr[4]++;
      else if (pEta >= 3.1 && pEta < 3.3 ) cntr[5]++;
      else if (pEta >= 3.3 && pEta < 3.5 ) cntr[6]++;
      else if (pEta >= 3.5 && pEta < 3.7 ) cntr[7]++;
      else if (pEta >= 3.7 && pEta < 3.9 ) cntr[8]++;  
      else if (pEta >= 3.9 && pEta < 4.1 ) cntr[9]++;
    }
  
    fhResponseAll->Fill(cntrr,kcntrr);
    fhTrueAll->Fill(kcntrr);
    fhMeasuredAll->Fill(cntrr);

  for(int j = 0; j < 10; j++ ) {
    fhResponse[i]->Fill(cntr[i],kcntr[i]);
    fhTrue[i]->Fill(kcntr[i]);
    fhMeasured[i]->Fill(cntr[i]);
  }
}


//________________________________________________________________________
void AliPMDAnalysisMCTaskSE::Terminate(Option_t *) 
{

  printf("Event number %d\n", fCntr);

     
}
Commit	Line	Data
642a954f	1	/**************************************************************************
	2	* Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
	3	* *
	4	* Author: The ALICE Off-line Project. *
	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	* Analysis Class Implimentation for the MC Truth and responce Matrix
	18	* Auther: Satyajit Jena, IIT Bombay \| sjena@cern.ch
	19	*
	20	* Mon Nov 22 19:54:27 CET 2010
	21	*
	22	* You may need some type of tunning to histograms and filling schemes,
	23	* Please contact me if you want something to change this code
	24	**************************************************************************/
	25
	26	#include "TChain.h"
	27	#include "TList.h"
	28	#include "TFile.h"
	29	#include "TTree.h"
	30	#include "TH1F.h"
	31	#include "TH1D.h"
	32	#include "TH2F.h"
	33	#include "TH3F.h"
	34	#include "TCanvas.h"
	35
	36	#include "AliAnalysisTask.h"
	37	#include "AliAnalysisManager.h"
	38	#include "AliVEvent.h"
	39	#include "AliAODEvent.h"
	40
	41	#include "AliMCEvent.h"
	42	#include "AliStack.h"
	43	#include "AliHeader.h"
	44	#include "AliESD.h"
	45	#include "AliESDEvent.h"
	46
	47	#include "AliMCEventHandler.h"
	48	#include "AliGenEventHeader.h"
	49	#include "AliGenPythiaEventHeader.h"
	50	#include "AliGenDPMjetEventHeader.h"
	51	#include "AliGenCocktailEventHeader.h"
	52	#include "AliPID.h"
	53
	54	#include "AliESDInputHandler.h"
	55	#include "AliPMDAnalysisMCTaskSE.h"
	56
	57	ClassImp(AliPMDAnalysisMCTaskSE)
	58
	59	//________________________________________________________________________
	60	AliPMDAnalysisMCTaskSE::AliPMDAnalysisMCTaskSE(const char *name)
	61	: AliAnalysisTaskSE(name),
	62	fPhysList(0),
	63	// fhCounter(0),
	64	// fhVtx(0),
65	// fhVty(0),
66	// fhVtz(0),
67	fhResponseAll(0),
68	fhTrueAll(0),
69	fhMeasuredAll(0)
70
71	{
72	for (Int_t i = 0; i < 10; i++) {
73	fhResponse[i] = NULL;
74	fhTrue[i] = NULL;
75	fhMeasured[i] = NULL;
76	}
77
78	DefineOutput(1, TList::Class());
79	}
80
81
82	//________________________________________________________________________
83	void AliPMDAnalysisMCTaskSE::UserCreateOutputObjects()
84	{
85	fPhysList = new TList();
86	fPhysList->SetOwner();
87
88	//___________ Counter ________________
89	// fCounter = new TH1D("hCounter","EVENT COUNTERS", 100, 0, 100);;
90	// fPhysList->Add(fCounter);
91
92	Int_t kMultBin = 100;
93	Float_t gMultMin = -0.5;
94	Float_t gMultMax = 99.5;
95
96
97	fhResponseAll = new TH2F("hResponseAll","Responce Matrix For All",
98	kMultBin,gMultMin,gMultMax,kMultBin,gMultMin,gMultMax);
99
100	fhResponseAll->GetXaxis()->SetTitle("Measured Multiplicity #gamma_{like}");
101	fhResponseAll->GetXaxis()->SetTitle("Ture Multiplicity #gamma");
102
103	fhTrueAll = new TH2F("hTrueAll","True Multiplicity of #gamma",kMultBin,gMultMin,gMultMax);
104	fhTrueAll->GetXaxis()->SetTitle("True Multiplicity #gamma");
105
106	fhMeasuredAll = new TH2F("hMeasuredAll","Measured Multiplicity of #gamma",kMultBin,gMultMin,gMultMax);
107	fhMeasuredAll->GetXaxis()->SetTitle("Measured Multiplicity #gamma");
108
109	fPhysList->Add(fhResponseAll);
110	fPhysList->Add(fhTrueAll);
111	fPhysList->Add(fhMeasuredAll);
112
113	for(Int_t i = 0; i < 10; i++) {
114	Float_t gbin = 2.3 + i*0.2;
115
116	char a[40];
117	sprintf(a,"hResponseEtaBin%2.1fto%2.1f",gbin, gbin+0.2);
118	char b[40];
119	sprintf(b,"hTrueEtaBin%2.1fto%2.1f",gbin, gbin+0.2);
120	char c[40];
121	sprintf(c,"hMeasuredEtaBin%2.1fto%2.1f",gbin, gbin+0.2);
122
123	char aa[40];
124	sprintf(aa," Responce Matrix from \eta - %2.1f to %2.1f",gbin, gbin+0.2);
125	char bb[40];
126	sprintf(bb," True Multiplicity of #gamma from \eta %2.1f to %2.1f",gbin, gbin+0.2);
127	char cc[40];
128	sprintf(cc," Measured multiplicity of gamma [#gamma_{like}] from \eta %2.1f to %2.1f",gbin, gbin+0.2);
129
130
131	fhResponse[i] = new TH2F(a,aa,kMultBin,gMultMin,gMultMax,kMultBin,gMultMin,gMultMax);
132	fhResponse[i]->GetXaxis()->SetTitle("Measured Multiplicity #gamma_{like}");
133	fhResponse[i]->GetYaxis()->SetTitle("Ture Multiplicity #gamma");
134
135
136	fhTrue[i] = new TH1F(b,bb,kMultBin,gMultMin,gMultMax);
137	fhResponse[i]->GetXaxis()->SetTitle("True Multiplicity #gamma");
138
139	fhMeasured[i] = new TH1F(c,cc,kMultBin,gMultMin,gMultMax);
140	fhMeasured[i]->GetXaxis()->SetTitle("True Multiplicity #gamma");
141
142	fPhysList->Add(fhResponse[i]);
143	fPhysList->Add(fhTrue[i]);
144	fPhysList->Add(fhMeasured[i]);
145
146	}
147
148	PostData(1,fPhysList);
149	}
150
151	//________________________________________________________________________
152	void AliPMDAnalysisMCTaskSE::UserExec(Option_t *)
153	{
154
155	fCntr++;
156
157	AliESDEvent* esd = dynamic_cast<AliESDEvent*>(InputEvent());
158	AliStack *stack = 0;
159	AliMCEventHandler *mcH = 0;
160	AliMCEvent* mcEvent = 0;
161
162	if(MCEvent()){
163	stack = MCEvent()->Stack();
164	mcH = (AliMCEventHandler*) ((AliAnalysisManager::GetAnalysisManager())->GetMCtruthEventHandler());
165	mcEvent = mcH->MCEvent();
166	AliLog::SetClassDebugLevel("AliMCEvent", AliLog::kWarning - AliLog::kDebug + 1);
167	}
168	if(!stack) return;
169	if(!esd) return;
170	if(!mcEvent)return;
171
172	EventByEvent(esd,mcEvent);
173
174	}
175
176	//----------------------------
177	void AliPMDAnalysisMCTaskSE::EventByEvent(AliESDEvent* esd, AliMCEvent* mcEvent)
178	{
179	Int_t trnn = mcEvent->GetNumberOfTracks();
180
181	Float_t McEta = 0.;
182
183	Int_t kcntr[0] ={0};
184	Int_t kcntrr = 0;
185
186	AliStack *stack = mcEvent->Stack();
187	Int_t sln = 0;
188	for (Int_t iTracks = 0; iTracks < trnn; iTracks++) {
189	AliMCParticle* track = (AliMCParticle*)mcEvent->GetTrack(iTracks);
190	// TParticle* part = track->Particle();
191
192	TParticle* part = stack->Particle(iTracks);
193	if (!part) continue;
194	Int_t pid = part->GetPdgCode();
195	TDatabasePDG* pdgDB = TDatabasePDG::Instance();
196	TParticlePDG* pdgPart = pdgDB->GetParticle(pid);
197	Float_t charge = (pdgPart ? pdgPart->Charge() : 0);
198
199	McEta = part->Eta();
200
201	if(McEta < 2.0 \|\| McEta > 4.2) continue;
202
203	kcntrr++;
204	if (McEta >= 2.1 && McEta < 2.3 ) kcntr[0]++;
205	else if (McEta >= 2.3 && McEta < 2.5 ) kcntr[1]++;
206	else if (McEta >= 2.5 && McEta < 2.7 ) kcntr[2]++;
207	else if (McEta >= 2.7 && McEta < 2.9 ) kcntr[3]++;
208	else if (McEta >= 2.9 && McEta < 3.1 ) kcntr[4]++;
209	else if (McEta >= 3.1 && McEta < 3.3 ) kcntr[5]++;
210	else if (McEta >= 3.3 && McEta < 3.5 ) kcntr[6]++;
211	else if (McEta >= 3.5 && McEta < 3.7 ) kcntr[7]++;
212	else if (McEta >= 3.7 && McEta < 3.9 ) kcntr[8]++;
213	else if (McEta >= 3.9 && McEta < 4.1 ) kcntr[9]++;
214
215	}
216
217	// ESD
218
219	Int_t npmdcl = esd->GetNumberOfPmdTracks();
220	Float_t rpxpy = 0.;
221	Float_t theta = 0.;
222	Float_t pEta = 0.;
223
224	Int_t cntr[10] = {0};
225	Int_t cntrr = 0;
226
227	while (npmdcl--)
228	{
229
230	AliESDPmdTrack *pmdtr = esd->GetPmdTrack(npmdcl);
231	Int_t det = pmdtr->GetDetector();
232	if( det == 1) continue;
233	Int_t ncell = pmdtr->GetClusterCells();
234	if(ncell < 2) continue;
235	Float_t adc = pmdtr->GetClusterADC();
236	if( adc < 216 ) continue;
237
238	Float_t clsX = pmdtr->GetClusterX();
239	Float_t clsY = pmdtr->GetClusterY();
240	Float_t clsZ = pmdtr->GetClusterZ();
241
242	rpxpy = TMath::Sqrt(clsXclsX + clsYclsY);
243	theta = TMath::ATan2(rpxpy,clsZ);
244	pEta = -TMath::Log(TMath::Tan(0.5*theta));
245	cntrr++;
246	if (pEta >= 2.1 && pEta < 2.3 ) cntr[0]++;
247	else if (pEta >= 2.3 && pEta < 2.5 ) cntr[1]++;
248	else if (pEta >= 2.5 && pEta < 2.7 ) cntr[2]++;
249	else if (pEta >= 2.7 && pEta < 2.9 ) cntr[3]++;
250	else if (pEta >= 2.9 && pEta < 3.1 ) cntr[4]++;
251	else if (pEta >= 3.1 && pEta < 3.3 ) cntr[5]++;
252	else if (pEta >= 3.3 && pEta < 3.5 ) cntr[6]++;
253	else if (pEta >= 3.5 && pEta < 3.7 ) cntr[7]++;
254	else if (pEta >= 3.7 && pEta < 3.9 ) cntr[8]++;
255	else if (pEta >= 3.9 && pEta < 4.1 ) cntr[9]++;
256	}
257
258	fhResponseAll->Fill(cntrr,kcntrr);
259	fhTrueAll->Fill(kcntrr);
260	fhMeasuredAll->Fill(cntrr);
261
262	for(int j = 0; j < 10; j++ ) {
263	fhResponse[i]->Fill(cntr[i],kcntr[i]);
264	fhTrue[i]->Fill(kcntr[i]);
265	fhMeasured[i]->Fill(cntr[i]);
266	}
267	}
268
269
270	//________________________________________________________________________
271	void AliPMDAnalysisMCTaskSE::Terminate(Option_t *)
272	{
273
274	printf("Event number %d\n", fCntr);
275
276
277	}