[u/mrichter/AliRoot.git] / T0 / AliT0MultiplicityTask.cxx

#include "TChain.h"
#include "TTree.h"
#include "TH1F.h"
#include "TH2F.h"
#include "TList.h"
#include "TClonesArray.h"
#include "TCanvas.h"
#include "TStyle.h"
#include "TLine.h"
#include "TText.h"
#include "TFile.h"
#include "TBenchmark.h"

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

#include "AliLog.h"
#include "AliVEvent.h"
#include "AliESDEvent.h"
#include "AliMCEvent.h"
#include "AliESDtrack.h"
#include "AliStack.h"
#include "AliTrackReference.h"

#include "AliT0MultiplicityTask.h"

ClassImp(AliT0MultiplicityTask)


AliT0MultiplicityTask::AliT0MultiplicityTask()
  : AliAnalysisTaskSE("AliT0MultiplicityTask"),
    fListOfHistos(0),
    fOrA(0),
    fOrC(0),
    fMean(0),
    fVertex(0),
    fTime(0),
    fAmp(0),
    fTotalMult(0),
    fMultRecTotal(0),
    fMultRecRealA(0),
    fMultRecRealC(0),
    fPrim(0),
    fRef(0),
    fRec(0)
{
  // Default constructor
  AliInfo("Default constructor AliT0MultiplicityTask");
  // Define input and output slots here
  // Input slot #0 works with a TChain
  DefineInput(0, TChain::Class());
  // Output slot #1 TList
  DefineOutput(1, TList::Class());
}


AliT0MultiplicityTask::AliT0MultiplicityTask(const char* name)
  : AliAnalysisTaskSE(name),
    fListOfHistos(0),
    fOrA(0),
    fOrC(0),
    fMean(0),
    fVertex(0),
    fTime(0),
    fAmp(0),   
    fTotalMult(0),
    fMultRecTotal(0),
    fMultRecRealA(0),
    fMultRecRealC(0),
    fPrim(0),
    fRef(0),
    fRec(0)
{
  // Constructor
  // Define input and output slots here
  // Input slot #0 works with a TChain
  DefineInput(0, TChain::Class());
  // Output slot #1 TList
  DefineOutput(1, TList::Class());
}


void AliT0MultiplicityTask::UserCreateOutputObjects() 
{
  // Create histograms
  // Called once
  // Create output container
  fListOfHistos = new TList();
  
 
  fOrA = new TH1F("hOrA"," T0 A", 100,2450,2700);
  fOrC = new TH1F("hOrC"," T0 C", 100,2450,2700);
  fMean= new TH1F("hMean"," T0 ",100, 12000, 13000);
  fVertex = new TH1F("hVertex","Z position of vertex",   100,-30, 30);
   
  fTime = new TH1F("Time", " Amp vs Time",100, 12000, 13000);
  fAmp = new TH1F("fAmp", " Amplitude", 100, 0, 200);

  //  fHighMult  = new TH1F("fHighMult", " events with high amp ",100, 0, 10);
  
  fTotalMult = new TH1F("fTotalMult","total multiplicity",500,0,5000);
  
  fMultRecRealA = new TH2F("fMultRecRealA"," ",100,0.,200,100,0,200); 
  fMultRecRealC = new TH2F("fMultRecRealC"," ",100,0,200,100,0,200); 
  fMultRecTotal = new TH2F("fMultRecTotal"," ",100,0,200,100,0,200); 
  
  fPrim = new TH1F("fPrim", " primary",100, 0, 200);
  fRef  = new TH1F("fRef", " from TR ",100, 0, 200);
  fRec  = new TH1F("fRec", " in ESD ",100, 0, 200);
   
  fListOfHistos->Add(fOrA);
  fListOfHistos->Add(fOrC);
  fListOfHistos->Add(fMean);
  fListOfHistos->Add(fVertex);
  fListOfHistos->Add(fAmp);
  fListOfHistos->Add(fTime);
  fListOfHistos->Add(fTotalMult);
  fListOfHistos->Add(fMultRecRealA);
  fListOfHistos->Add(fMultRecRealC);
  fListOfHistos->Add(fMultRecTotal);
  fListOfHistos->Add(fPrim);
  fListOfHistos->Add(fRef);
  fListOfHistos->Add(fRec);
}  


void AliT0MultiplicityTask::UserExec(Option_t *) 
{
  // Main loop
  // Called for each event

    
  // MC information
  AliMCEvent* mcEvent = MCEvent();
  if (!mcEvent) {
    Printf("ERROR: Could not retrieve MC event");
    return;
  }
  
  
  Int_t primaryA=0;
  Int_t primaryC=0;
  Int_t numPrim = 0;
  Int_t refT0=0;

  AliStack* stack = mcEvent->Stack();
  // printf("AliT0MultiplicityAnalysis: Number of tracks on stack %5d\n", stack->GetNtrack());
  Int_t nTracks  = mcEvent->GetNumberOfTracks();
  for (Int_t ipart=0; ipart<nTracks; ipart++)
    {
      //    TParticle* particle = stack->Particle(ipart);
      AliMCParticle* track = mcEvent->GetTrack(ipart);
      if (!track) {
	Printf("ERROR: Could not receive track %d (mc loop)", ipart);
	continue;
      }
      Int_t label = track->GetLabel();
      if(stack->IsPhysicalPrimary(label) == kFALSE)
	continue;
      
      numPrim++;
      Double_t eta=track->Eta();
      
      if(eta<-2.97 && eta>-3.28)  primaryC++;
      if (eta >4.61 && eta<4.92)  primaryA++;
    
       // Loop over Track References
       AliTrackReference* trackRef = 0;

      for (Int_t iTrackRef = 0; iTrackRef  < track->GetNumberOfTrackReferences(); iTrackRef++) {
	trackRef = track->GetTrackReference(iTrackRef);
	if(trackRef) {
	  Int_t detectorId = trackRef->DetectorId(); 
	  if (detectorId == 7)  refT0++;
	}      
      }
    }
  fTotalMult->Fill(numPrim);
  fPrim->Fill(primaryC+primaryA);
  fRef->Fill(refT0);
  //     printf (" tracks %i primaries %d \n",nTracks, numPrim);
 

  // ESD information  
  AliVEvent* event = InputEvent();
  if (!event) {
    Printf("ERROR: Could not retrieve event");
    return;
  }
    
  AliESDEvent* esd = dynamic_cast<AliESDEvent*>(event);

  Bool_t eventTriggered = false;
  ULong64_t triggerMask = esd->GetTriggerMask();
  // definitions from p-p.cfg
  ULong64_t spdFO = (1 << 14);
  ULong64_t v0left = (1 << 11);
  ULong64_t v0right = (1 << 12);

  if (triggerMask & spdFO || ((triggerMask & v0left) || (triggerMask & v0right)))
    eventTriggered == true;
 
    //= AliPWG0Helper::IsEventTriggered(esd, AliPWG0Helper::kMB1);
  //  printf("!!!!! eventTriggered %i",eventTriggered);
  // if(!eventTriggered) return;
   Double_t besttimeC=99999.;
  Double_t besttimeA=99999.;
  
//  Float_t coefA = 0.891466;
//  Float_t coefC = 0.922772;
  Float_t coefA = 1;
  Float_t coefC = 1;
  Float_t sumampA=0;
  Float_t sumampC=0;
  Int_t highA=0, highC=0;
  const Double_t *amp = esd->GetT0amplitude();
  const Double_t *time = esd->GetT0time();
  Float_t vertex = esd->GetT0zVertex();
  if(vertex<999) fVertex->Fill(vertex);
  Float_t   mean = esd->GetT0();
  fMean->Fill(mean);
  
  for (Int_t i=0; i<12; i++) {
    sumampC += amp[i];  
    fTime->Fill(amp[i],time[i]); 
    fAmp->Fill(amp[i]); 
    if(time[i]<besttimeC && time[i]>0) besttimeC=time[i]; //timeC
  }
  fOrC->Fill(besttimeC);
  
  for (Int_t i=12; i<24; i++){
    sumampA += amp[i];  
    fAmp->Fill(amp[i]); 
    fTime->Fill(time[i]);
    if(time[i]<besttimeA && time[i]>0) besttimeA=time[i]; //timeC
  }

  fOrA->Fill(besttimeA);
  
  fMultRecRealA ->Fill(primaryA, sumampA*coefA);
  fMultRecRealC ->Fill(primaryC, sumampC*coefC);
  fMultRecTotal->Fill(numPrim, sumampA*coefA +  sumampC*coefC);
  fRec->Fill(Int_t(sumampA+sumampC));
 
  // Post output data.
  PostData(1, fListOfHistos);
}      


void AliT0MultiplicityTask::Terminate(Option_t *) 
{
  // Draw result to the screen
  // Called once at the end of the query 
  //  printf(" AliT0MultiplicityTask::Terminate ");
  
  fListOfHistos = dynamic_cast<TList*>(GetOutputData(1));
  if (!fListOfHistos) {
    Printf("ERROR: fListOfHistos not available");
    return;
  }
  //  printf(" before cast ");
  fOrA = dynamic_cast<TH1F*>(fListOfHistos->At(0));
  fOrC = dynamic_cast<TH1F*>(fListOfHistos->At(1));
  fMean  = dynamic_cast<TH1F*>(fListOfHistos->At(2));  
  fVertex  = dynamic_cast<TH1F*>(fListOfHistos->At(3));
  fAmp = dynamic_cast<TH1F*>(fListOfHistos->At(4));
  fTime = dynamic_cast<TH1F*>(fListOfHistos->At(5));
  fTotalMult = dynamic_cast<TH1F*>(fListOfHistos->At(6));
  fMultRecTotal = dynamic_cast<TH2F*>(fListOfHistos->At(7));
  fMultRecRealA = dynamic_cast<TH2F*>(fListOfHistos->At(8));
  fMultRecRealC= dynamic_cast<TH2F*>(fListOfHistos->At(9));
  fPrim= dynamic_cast<TH1F*>(fListOfHistos->At(10));
  fRef= dynamic_cast<TH1F*>(fListOfHistos->At(11));
  fRec= dynamic_cast<TH1F*>(fListOfHistos->At(12));

  TFile fc("MultHist.root\n", "RECREATE");
  //  printf(" File MultHist.root recreated\n");
  fOrC->Write();
  //  printf("1\n");
  fOrA->Write();
  //  printf("2\n");
  fMean->Write();
  //  printf("3\n");
  fVertex->Write();
  //  printf("4\n");
  fAmp->Write();
  //  printf("5\n");
  fTime->Write();
 fTotalMult->Write();
 //  printf("6\n");
 fMultRecTotal->Write();
 //  printf("7\n");
 fMultRecRealA->Write();
 // printf("8\n");
  fMultRecRealC->Write();
  // printf("9\n");
  fPrim->Write();
  fRef->Write();
  fRec->Write();
  
  fc.Close();
 printf(" fc.Close()\n");

}
Commit	Line	Data
c137086f	1	#include "TChain.h"
	2	#include "TTree.h"
	3	#include "TH1F.h"
	4	#include "TH2F.h"
	5	#include "TList.h"
	6	#include "TClonesArray.h"
	7	#include "TCanvas.h"
	8	#include "TStyle.h"
	9	#include "TLine.h"
	10	#include "TText.h"
	11	#include "TFile.h"
	12	#include "TBenchmark.h"
	13
	14	#include "AliAnalysisTask.h"
	15	#include "AliAnalysisManager.h"
	16
	17	#include "AliLog.h"
	18	#include "AliVEvent.h"
	19	#include "AliESDEvent.h"
	20	#include "AliMCEvent.h"
	21	#include "AliESDtrack.h"
	22	#include "AliStack.h"
	23	#include "AliTrackReference.h"
	24
	25	#include "AliT0MultiplicityTask.h"
	26
	27	ClassImp(AliT0MultiplicityTask)
	28
	29
	30	AliT0MultiplicityTask::AliT0MultiplicityTask()
	31	: AliAnalysisTaskSE("AliT0MultiplicityTask"),
	32	fListOfHistos(0),
	33	fOrA(0),
	34	fOrC(0),
	35	fMean(0),
	36	fVertex(0),
	37	fTime(0),
	38	fAmp(0),
	39	fTotalMult(0),
	40	fMultRecTotal(0),
	41	fMultRecRealA(0),
	42	fMultRecRealC(0),
	43	fPrim(0),
	44	fRef(0),
	45	fRec(0)
	46	{
	47	// Default constructor
	48	AliInfo("Default constructor AliT0MultiplicityTask");
	49	// Define input and output slots here
	50	// Input slot #0 works with a TChain
	51	DefineInput(0, TChain::Class());
	52	// Output slot #1 TList
	53	DefineOutput(1, TList::Class());
	54	}
	55
	56
	57	AliT0MultiplicityTask::AliT0MultiplicityTask(const char* name)
	58	: AliAnalysisTaskSE(name),
	59	fListOfHistos(0),
	60	fOrA(0),
	61	fOrC(0),
	62	fMean(0),
	63	fVertex(0),
	64	fTime(0),
65	fAmp(0),
66	fTotalMult(0),
67	fMultRecTotal(0),
68	fMultRecRealA(0),
69	fMultRecRealC(0),
70	fPrim(0),
71	fRef(0),
72	fRec(0)
73	{
74	// Constructor
75	// Define input and output slots here
76	// Input slot #0 works with a TChain
77	DefineInput(0, TChain::Class());
78	// Output slot #1 TList
79	DefineOutput(1, TList::Class());
80	}
81
82
83	void AliT0MultiplicityTask::UserCreateOutputObjects()
84	{
85	// Create histograms
86	// Called once
87	// Create output container
88	fListOfHistos = new TList();
89
90
91	fOrA = new TH1F("hOrA"," T0 A", 100,2450,2700);
92	fOrC = new TH1F("hOrC"," T0 C", 100,2450,2700);
93	fMean= new TH1F("hMean"," T0 ",100, 12000, 13000);
94	fVertex = new TH1F("hVertex","Z position of vertex", 100,-30, 30);
95
96	fTime = new TH1F("Time", " Amp vs Time",100, 12000, 13000);
97	fAmp = new TH1F("fAmp", " Amplitude", 100, 0, 200);
98
99	// fHighMult = new TH1F("fHighMult", " events with high amp ",100, 0, 10);
100
101	fTotalMult = new TH1F("fTotalMult","total multiplicity",500,0,5000);
102
103	fMultRecRealA = new TH2F("fMultRecRealA"," ",100,0.,200,100,0,200);
104	fMultRecRealC = new TH2F("fMultRecRealC"," ",100,0,200,100,0,200);
105	fMultRecTotal = new TH2F("fMultRecTotal"," ",100,0,200,100,0,200);
106
107	fPrim = new TH1F("fPrim", " primary",100, 0, 200);
108	fRef = new TH1F("fRef", " from TR ",100, 0, 200);
109	fRec = new TH1F("fRec", " in ESD ",100, 0, 200);
110
111	fListOfHistos->Add(fOrA);
112	fListOfHistos->Add(fOrC);
113	fListOfHistos->Add(fMean);
114	fListOfHistos->Add(fVertex);
115	fListOfHistos->Add(fAmp);
116	fListOfHistos->Add(fTime);
117	fListOfHistos->Add(fTotalMult);
118	fListOfHistos->Add(fMultRecRealA);
119	fListOfHistos->Add(fMultRecRealC);
120	fListOfHistos->Add(fMultRecTotal);
121	fListOfHistos->Add(fPrim);
122	fListOfHistos->Add(fRef);
123	fListOfHistos->Add(fRec);
124	}
125
126
127	void AliT0MultiplicityTask::UserExec(Option_t *)
128	{
129	// Main loop
130	// Called for each event
131
132
133	// MC information
134	AliMCEvent* mcEvent = MCEvent();
135	if (!mcEvent) {
136	Printf("ERROR: Could not retrieve MC event");
137	return;
138	}
139
140
141	Int_t primaryA=0;
142	Int_t primaryC=0;
143	Int_t numPrim = 0;
144	Int_t refT0=0;
145
146	AliStack* stack = mcEvent->Stack();
147	// printf("AliT0MultiplicityAnalysis: Number of tracks on stack %5d\n", stack->GetNtrack());
148	Int_t nTracks = mcEvent->GetNumberOfTracks();
149	for (Int_t ipart=0; ipart<nTracks; ipart++)
150	{
151	// TParticle* particle = stack->Particle(ipart);
152	AliMCParticle* track = mcEvent->GetTrack(ipart);
153	if (!track) {
154	Printf("ERROR: Could not receive track %d (mc loop)", ipart);
155	continue;
156	}
157	Int_t label = track->GetLabel();
158	if(stack->IsPhysicalPrimary(label) == kFALSE)
159	continue;
160
161	numPrim++;
162	Double_t eta=track->Eta();
163
164	if(eta<-2.97 && eta>-3.28) primaryC++;
165	if (eta >4.61 && eta<4.92) primaryA++;
166
167	// Loop over Track References
168	AliTrackReference* trackRef = 0;
169
170	for (Int_t iTrackRef = 0; iTrackRef < track->GetNumberOfTrackReferences(); iTrackRef++) {
171	trackRef = track->GetTrackReference(iTrackRef);
172	if(trackRef) {
173	Int_t detectorId = trackRef->DetectorId();
174	if (detectorId == 7) refT0++;
175	}
176	}
177	}
178	fTotalMult->Fill(numPrim);
179	fPrim->Fill(primaryC+primaryA);
180	fRef->Fill(refT0);
181	// printf (" tracks %i primaries %d \n",nTracks, numPrim);
182
183
184	// ESD information
185	AliVEvent* event = InputEvent();
186	if (!event) {
187	Printf("ERROR: Could not retrieve event");
188	return;
189	}
190
191	AliESDEvent* esd = dynamic_cast<AliESDEvent*>(event);
192
193	Bool_t eventTriggered = false;
194	ULong64_t triggerMask = esd->GetTriggerMask();
195	// definitions from p-p.cfg
196	ULong64_t spdFO = (1 << 14);
197	ULong64_t v0left = (1 << 11);
198	ULong64_t v0right = (1 << 12);
199
200	if (triggerMask & spdFO \|\| ((triggerMask & v0left) \|\| (triggerMask & v0right)))
201	eventTriggered == true;
202
203	//= AliPWG0Helper::IsEventTriggered(esd, AliPWG0Helper::kMB1);
204	// printf("!!!!! eventTriggered %i",eventTriggered);
205	// if(!eventTriggered) return;
206	Double_t besttimeC=99999.;
207	Double_t besttimeA=99999.;
208
209	// Float_t coefA = 0.891466;
210	// Float_t coefC = 0.922772;
211	Float_t coefA = 1;
212	Float_t coefC = 1;
213	Float_t sumampA=0;
214	Float_t sumampC=0;
215	Int_t highA=0, highC=0;
216	const Double_t *amp = esd->GetT0amplitude();
217	const Double_t *time = esd->GetT0time();
218	Float_t vertex = esd->GetT0zVertex();
219	if(vertex<999) fVertex->Fill(vertex);
220	Float_t mean = esd->GetT0();
221	fMean->Fill(mean);
222
223	for (Int_t i=0; i<12; i++) {
224	sumampC += amp[i];
225	fTime->Fill(amp[i],time[i]);
226	fAmp->Fill(amp[i]);
227	if(time[i]<besttimeC && time[i]>0) besttimeC=time[i]; //timeC
228	}
229	fOrC->Fill(besttimeC);
230
231	for (Int_t i=12; i<24; i++){
232	sumampA += amp[i];
233	fAmp->Fill(amp[i]);
234	fTime->Fill(time[i]);
235	if(time[i]<besttimeA && time[i]>0) besttimeA=time[i]; //timeC
236	}
237
238	fOrA->Fill(besttimeA);
239
240	fMultRecRealA ->Fill(primaryA, sumampA*coefA);
241	fMultRecRealC ->Fill(primaryC, sumampC*coefC);
242	fMultRecTotal->Fill(numPrim, sumampAcoefA + sumampCcoefC);
243	fRec->Fill(Int_t(sumampA+sumampC));
244
245	// Post output data.
246	PostData(1, fListOfHistos);
247	}
248
249
250	void AliT0MultiplicityTask::Terminate(Option_t *)
251	{
252	// Draw result to the screen
253	// Called once at the end of the query
254	// printf(" AliT0MultiplicityTask::Terminate ");
255
256	fListOfHistos = dynamic_cast<TList*>(GetOutputData(1));
257	if (!fListOfHistos) {
258	Printf("ERROR: fListOfHistos not available");
259	return;
260	}
261	// printf(" before cast ");
262	fOrA = dynamic_cast<TH1F*>(fListOfHistos->At(0));
263	fOrC = dynamic_cast<TH1F*>(fListOfHistos->At(1));
264	fMean = dynamic_cast<TH1F*>(fListOfHistos->At(2));
265	fVertex = dynamic_cast<TH1F*>(fListOfHistos->At(3));
266	fAmp = dynamic_cast<TH1F*>(fListOfHistos->At(4));
267	fTime = dynamic_cast<TH1F*>(fListOfHistos->At(5));
268	fTotalMult = dynamic_cast<TH1F*>(fListOfHistos->At(6));
269	fMultRecTotal = dynamic_cast<TH2F*>(fListOfHistos->At(7));
270	fMultRecRealA = dynamic_cast<TH2F*>(fListOfHistos->At(8));
271	fMultRecRealC= dynamic_cast<TH2F*>(fListOfHistos->At(9));
272	fPrim= dynamic_cast<TH1F*>(fListOfHistos->At(10));
273	fRef= dynamic_cast<TH1F*>(fListOfHistos->At(11));
274	fRec= dynamic_cast<TH1F*>(fListOfHistos->At(12));
275
276	TFile fc("MultHist.root\n", "RECREATE");
277	// printf(" File MultHist.root recreated\n");
278	fOrC->Write();
279	// printf("1\n");
280	fOrA->Write();
281	// printf("2\n");
282	fMean->Write();
283	// printf("3\n");
284	fVertex->Write();
285	// printf("4\n");
286	fAmp->Write();
287	// printf("5\n");
288	fTime->Write();
289	fTotalMult->Write();
290	// printf("6\n");
291	fMultRecTotal->Write();
292	// printf("7\n");
293	fMultRecRealA->Write();
294	// printf("8\n");
295	fMultRecRealC->Write();
296	// printf("9\n");
297	fPrim->Write();
298	fRef->Write();
299	fRec->Write();
300
301	fc.Close();
302	printf(" fc.Close()\n");
303
304	}