DA find only coorect peak
[u/mrichter/AliRoot.git] / T0 / AliT0MultiplicityTask.cxx
CommitLineData
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
27ClassImp(AliT0MultiplicityTask)
28
29
30AliT0MultiplicityTask::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
57AliT0MultiplicityTask::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
83void 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
127void 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, sumampA*coefA + sumampC*coefC);
243 fRec->Fill(Int_t(sumampA+sumampC));
244
245 // Post output data.
246 PostData(1, fListOfHistos);
247}
248
249
250void 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}