AliPhysics  323987d (323987d)
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
AliAnalysisTaskJetsEvshape.cxx
Go to the documentation of this file.
1 // ROOT
2 #include "TFile.h"
3 #include "TList.h"
4 #include "TH1.h"
5 #include "TH2.h"
6 #include "TH3.h"
7 #include "TF1.h"
8 #include "TFormula.h"
9 #include "TRandom.h"
10 #include "TSpline.h"
11 
12 // analysis framework
13 #include "AliAnalysisManager.h"
14 #include "AliInputEventHandler.h"
15 #include "AliVEvent.h"
16 #include "AliVTrack.h"
17 #include "AliVVertex.h"
18 #include "AliVMultiplicity.h"
19 
20 // MC stuff
21 #include "AliMCEvent.h"
22 #include "AliGenPythiaEventHeader.h"
23 #include "AliStack.h"
24 
25 // ESD stuff
26 #include "AliESDEvent.h"
27 #include "AliESDInputHandler.h"
28 #include "AliESDtrack.h"
29 #include "AliESDtrackCuts.h"
30 
31 // AOD stuff
32 #include "AliAODEvent.h"
33 #include "AliAODJet.h"
34 #include "AliAODTrack.h"
35 
36 // EMCAL framework and jet tasks
37 #include "AliJetContainer.h"
38 #include "AliParticleContainer.h"
39 #include "AliClusterContainer.h"
40 #include "AliEmcalJet.h"
41 
44 
45 #include <iostream>
46 #include <cmath>
47 
49  AliAnalysisTaskEmcalJet(name, kTRUE),
50  fMCEvent(0x0),
51  fESDEvent(0x0),
52  fAODEvent(0x0),
53  fRunNumber(-1),
54  fJetsCont(0),
55  fTracksCont(0),
56  fCaloClustersCont(0),
57  fThrust(0.),
58  fSphericityT(0.),
59  fUseMC(kTRUE),
60  fEtaMaxForEvshape(.9),
61  fOutputList(0x0),
62  fShortTaskId("jets_evshape")
63 {
64  // default ctor
65 
67 
68  DefineOutput(kOutputEmcal, TList::Class());
69  AliInfo(Form("creating output slot #%i", kOutputTask));
70  DefineOutput(kOutputTask, TList::Class());
71 }
72 
74 {
75  // dtor
76 }
77 
79 {
80  // create user output objects
81  AliInfo("creating output objects");
82 
83  // common EMCAL framework
85  PostData(kOutputEmcal, fOutput);
87  printf("setup jet container: %p, from available:\n", fJetsCont);
88  fJetCollArray.Print();
89  printf("end\n");
90  if(fJetsCont) { //get particles and clusters connected to jets
94  } else { //no jets, just analysis tracks and clusters
97  }
98  if(fTracksCont)
99  fTracksCont->SetClassName("AliVParticle");
101  fCaloClustersCont->SetClassName("AliVCluster");
102 
103  // setup list
105  fOutputList = new TList();
106  fOutputList->SetOwner();
107 
108  // setup histograms
109  TH1 *hist;
110  TH1 *histStat = AddHistogram(ID(kHistStat), "event statistics;;counts",
111  kStatLast-1, .5, kStatLast-.5);
112  histStat->GetXaxis()->SetBinLabel(ID(kStatSeen));
113  histStat->GetXaxis()->SetBinLabel(ID(kStatTrg));
114  histStat->GetXaxis()->SetBinLabel(ID(kStatUsed));
115  histStat->GetXaxis()->SetBinLabel(ID(kStatEvCuts));
116 
117  AddHistogram(ID(kHistJetPt), "jet spectrum;p_{T}^{jet,ch} (GeV/#it{c});counts",
118  40, 0., 40.);
119  AddHistogram(ID(kHistMult), "multiplicity;multiplicity;counts",
120  100, 0., 4000.);
121  AddHistogram(ID(kHistJetPtVsMult), "jet p_{T} vs mult;multiplicity;p_{T}^{jet,ch} (GeV/#it{c})",
122  100, 0., 4000.,
123  40, 0., 40.);
124  AddHistogram(ID(kHistJetPtMultEvshape), "jet spectrum;p_{T}^{jet,ch} (GeV/#it{c});multiplicity;S_{T}",
125  40, 0., 40., 100, 0., 500., 50, 0., 1.);
126  AddHistogram(ID(kHistThrust), "thrust;T;counts",
127  50, 0., 1.);
128  AddHistogram(ID(kHistSphericityT), "transverse sphericity;S_{T};counts",
129  50, 0., 1.);
130  AddHistogram(ID(kHistSphericityTVsMult), "transverse sphericity vs mult;multiplicity;S_{T};counts",
131  100, 0., 4000.,
132  50, 0., 1.);
133 
134  PostData(kOutputTask, fOutputList);
135 }
136 
138 {
139  // actions to be taken upon notification about input file change
140 
141  return AliAnalysisTaskEmcalJet::Notify();
142 }
143 
145 {
146  // actions at task termination
147 
148  AliAnalysisTaskEmcalJet::Terminate(option);
149 }
150 
152 {
153  AliAnalysisTaskEmcalJet::PrintTask(option, indent);
154 
155  std::cout << std::setw(indent) << " " << "nothing to say: " << std::endl;
156 }
157 
159 {
160  Bool_t eventGood = kTRUE;
161 
162  // setup pointers to input data (null if unavailable)
163  // mcEvent: MC input
164  // esdEvent: ESD input
165  // outEvent: AOD output
166  // aodEvent: AOD input if available, otherwise AOD output
167 
168  fMCEvent = this->MCEvent();
169  fESDEvent = dynamic_cast<AliESDEvent*>(this->InputEvent()); // could also be AOD input
170  AliAODEvent* outEvent = this->AODEvent();
171  fAODEvent = dynamic_cast<AliAODEvent*> (this->InputEvent());
172  if (!fAODEvent)
173  fAODEvent = outEvent;
174 
175  if ((fDebug > 0) && fESDEvent)
176  printf("event: %s-%06i\n", CurrentFileName(), fESDEvent->GetEventNumberInFile());
177 
178  // check for run change
179  if (fRunNumber != InputEvent()->GetRunNumber()) {
180  fRunNumber = InputEvent()->GetRunNumber();
181  }
182 
183  return eventGood;
184 }
185 
187 {
188  return kTRUE;
189 }
190 
192 {
194 }
195 
197 {
198  // record number of sampled events and detect trigger contributions
200 
201  // so far, no trigger selection, we accept all
203 
204  // prepare the event
205  if (PrepareEvent()) {
206  // InputEvent()->GetList()->ls();
207 
208  // multiplicity selection
209 
210  // event shape selection
211  if (!CalculateThrust())
212  return kFALSE;
213 
214  if (!CalculateSphericityT())
215  return kFALSE;
216 
217  // here we have passed the event cuts
219 
220  return kTRUE;
221  } else {
222  return kFALSE;
223  }
224 }
225 
227 {
229 
230  const AliStack *stack = fMCEvent ? fMCEvent->Stack() : 0x0;
231  const AliVMultiplicity *mult = InputEvent()->GetMultiplicity();
232  const Int_t nTracklets =
233  stack ? stack->GetNprimary() :
234  mult ? mult->GetNumberOfTracklets() :
235  -1;
236  FillH1(kHistMult, nTracklets);
240 
241  if (fJetsCont) {
242  // const Int_t nJets = fJetsCont->GetNJets();
243  // const Int_t nAcceptedJets = fJetsCont->GetNAcceptedJets();
244 
245  fJetsCont->ResetCurrentID();
246  while (AliEmcalJet *jet = fJetsCont->GetNextAcceptJet()) {
247  FillH1(kHistJetPt, jet->Pt());
248  FillH2(kHistJetPtVsMult, nTracklets, jet->Pt());
249  FillH3(kHistJetPtMultEvshape, jet->Pt(), nTracklets, fSphericityT);
250  }
251  }
252 
253  return kTRUE;
254 }
255 
257 {
258  // actual work
259 
261 
262  CleanUpEvent();
263 
264  PostData(kOutputEmcal, fOutput);
265  PostData(kOutputTask, fOutputList);
266 }
267 
269 {
270  const AliVEvent *event = fUseMC ? MCEvent() : InputEvent();
271 
272  // to be added: track cuts for event shape observables
273  const Int_t nTracksAll = event ? event->GetNumberOfTracks() : 0;
274  if (nTracksAll == 0)
275  return kFALSE;
276 
277  std::vector<TVector3> p;
278  for (Int_t iTrack = 0; iTrack < nTracksAll; ++iTrack) {
279  AliVParticle *track = event->GetTrack(iTrack);
280  if (fUseMC && !MCEvent()->IsPhysicalPrimary(iTrack))
281  continue;
282  if (TMath::Abs(track->Eta()) > fEtaMaxForEvshape)
283  continue;
284  Double_t pi[3];
285  track->PxPyPz(pi);
286  TVector3 p3(pi);
287  p.push_back(p3);
288  }
289 
290  // thrust calculation following H. Yamamoto, J. Comp. Physics 52, 1983
291  const Int_t nTracks = p.size();
292  Char_t *signset1 = new Char_t[nTracks];
293  Char_t *signset2 = new Char_t[nTracks];
294  Char_t *sign = signset1;
295  Char_t *signmax = signset2;
296  Double_t thrustmax = 0;
297  Double_t norm = 0;
298 
299  // looping over pairs of tracks (i,j) w. j > i
300  TVector3 Pn;
301  TVector3 Pthrust;
302  TVector3 Ptemp[4];
303  for (Int_t i = 0; i < nTracks; i++) {
304  norm += p[i].Mag();
305  for (Int_t j = i+1; j < nTracks; j++) {
306  Pn = p[i].Cross(p[j]);
307 
308  // calculate hypothetical thrust
309  Double_t thrust;
310  Pthrust.SetXYZ(0., 0., 0.);
311  for (Int_t iTrack = 0; iTrack < nTracks; iTrack++) {
312  if ((iTrack == i) || (iTrack == j))
313  continue;
314  sign[iTrack] = (Pn.Dot(p[iTrack]) > 0) ? 1 : -1;
315  if (sign[iTrack] > 0)
316  Pthrust += p[iTrack];
317  else
318  Pthrust -= p[iTrack];
319  }
320 
321  Ptemp[0] = Ptemp[1] = Ptemp[2] = Ptemp[3] = Pthrust;
322  Ptemp[0] += p[i]; Ptemp[0] -= p[j];
323  Ptemp[1] += p[i]; Ptemp[0] += p[j];
324  Ptemp[2] -= p[i]; Ptemp[0] -= p[j];
325  Ptemp[3] -= p[i]; Ptemp[0] += p[j];
326 
327  Double_t thrustmaxtemp = 0;
328  Int_t kmax = -1;
329  for (Int_t k = 0; k < 4; k++) {
330  if (Ptemp[k].Mag() > thrustmaxtemp) {
331  kmax = k;
332  thrustmaxtemp = Ptemp[k].Mag();
333  }
334  }
335  if (thrustmaxtemp > thrustmax) {
336  Char_t *temp = signmax;
337  signmax = sign;
338  sign = temp;
339  thrustmax = thrustmaxtemp;
340  }
341  }
342  }
343 
344  if (norm > 0) {
345  fThrust = thrustmax / norm;
346  return kTRUE;
347  } else {
348  return kFALSE;
349  }
350 }
351 
353 {
354  const AliVEvent *event = fUseMC ? MCEvent() : InputEvent();
355 
356  const Int_t nTracks = event->GetNumberOfTracks();
357  if (nTracks == 0)
358  return kFALSE;
359 
360  // calculate transverse sphericity
361  Float_t smatrix[3] = { 0. };
362  Float_t sumpt = 0.;
363  for (Int_t iTrack = 0; iTrack < nTracks; iTrack++) {
364  AliVParticle *part= event->GetTrack(iTrack);
365  if (fUseMC && !MCEvent()->IsPhysicalPrimary(iTrack))
366  continue;
367  if (TMath::Abs(part->Eta()) > fEtaMaxForEvshape)
368  continue;
369  const Float_t pt = part->Pt();
370  const Float_t px = part->Px();
371  const Float_t py = TMath::Sqrt(pt * pt - px * px);
372  smatrix[0] += px * px / pt;
373  smatrix[1] += px * py / pt;
374  smatrix[2] += py * py / pt;
375  sumpt += part->Pt();
376  }
377  if (sumpt > 0) {
378  smatrix[0] /= sumpt;
379  smatrix[1] /= sumpt;
380  smatrix[2] /= sumpt;
381  }
382 
383  // calculate eigenvalues
384  Float_t lambda1 = ((smatrix[0] + smatrix[2]) + TMath::Sqrt((smatrix[0] + smatrix[2])*(smatrix[0] + smatrix[2]) - 4 * (smatrix[0] * smatrix[2] - smatrix[1] * smatrix[1]))) / 2;
385  Float_t lambda2 = ((smatrix[0] + smatrix[2]) - TMath::Sqrt((smatrix[0] + smatrix[2])*(smatrix[0] + smatrix[2]) - 4 * (smatrix[0] * smatrix[2] - smatrix[1] * smatrix[1]))) / 2;
386 
387  if (TMath::Abs(lambda1 + lambda2) > 1.e-9) {
388  fSphericityT = 2. * TMath::Min(lambda1, lambda2) / (lambda1 + lambda2);
389  return kTRUE;
390  } else {
391  return kFALSE;
392  }
393 }
394 
395 // ----- histogram management -----
397  Int_t xbins, Float_t xmin, Float_t xmax,
398  Int_t binType)
399 {
400  TString hName;
401  hName.Form("%s_%s", fShortTaskId, hid);
402  hName.ToLower();
403  TH1 *h = 0x0;
404  if (binType == 0)
405  h = new TH1I(hName.Data(), title,
406  xbins, xmin, xmax);
407  else
408  h = new TH1F(hName.Data(), title,
409  xbins, xmin, xmax);
410  GetHistogram(hist) = h;
411  fOutputList->Add(h);
412  return h;
413 }
414 
416  Int_t xbins, Float_t xmin, Float_t xmax,
417  Int_t ybins, Float_t ymin, Float_t ymax,
418  Int_t binType)
419 {
420  TString hName;
421  hName.Form("%s_%s", fShortTaskId, hid);
422  hName.ToLower();
423  TH1 *h = 0x0;
424  if (binType == 0)
425  h = GetHistogram(hist) = new TH2I(hName.Data(), title,
426  xbins, xmin, xmax,
427  ybins, ymin, ymax);
428  else
429  h = GetHistogram(hist) = new TH2F(hName.Data(), title,
430  xbins, xmin, xmax,
431  ybins, ymin, ymax);
432  fOutputList->Add(h);
433  return (TH2*) h;
434 }
435 
437  Int_t xbins, Float_t xmin, Float_t xmax,
438  Int_t ybins, Float_t ymin, Float_t ymax,
439  Int_t zbins, Float_t zmin, Float_t zmax,
440  Int_t binType)
441 {
442  TString hName;
443  hName.Form("%s_%s", fShortTaskId, hid);
444  hName.ToLower();
445  TH1 *h = 0x0;
446  if (binType == 0)
447  h = GetHistogram(hist) = new TH3I(hName.Data(), title,
448  xbins, xmin, xmax,
449  ybins, ymin, ymax,
450  zbins, zmin, zmax);
451  else
452  h = GetHistogram(hist) = new TH3F(hName.Data(), title,
453  xbins, xmin, xmax,
454  ybins, ymin, ymax,
455  zbins, zmin, zmax);
456  fOutputList->Add(h);
457  return (TH3*) h;
458 }
const Double_t ymax
Definition: AddTaskCFDStar.C:7
TH1 * AddHistogram(Hist_t hist, const char *hid, TString title, Int_t xbins, Float_t xmin, Float_t xmax, Int_t binType=1)
double Double_t
Definition: External.C:58
Definition: External.C:260
Definition: External.C:236
const char * title
Definition: MakeQAPdf.C:26
AliJetContainer * GetJetContainer(Int_t i=0) const
Definition: External.C:244
char Char_t
Definition: External.C:18
AliClusterContainer * GetClusterContainer() const
void FillH1(Hist_t hist, Float_t x, Float_t weight=1., Int_t idx=0)
virtual void PrintTask(Option_t *option, Int_t indent) const
const char * fShortTaskId
pointers to histogram
virtual void Terminate(const Option_t *option)
AliParticleContainer * GetParticleContainer(Int_t i=0) const
AliParticleContainer * GetParticleContainer() const
int Int_t
Definition: External.C:63
Definition: External.C:204
float Float_t
Definition: External.C:68
void FillH3(Hist_t hist, Float_t x, Float_t y, Float_t z, Float_t weight=1., Int_t idx=0)
const Double_t zmin
#define ID(x)
AliClusterContainer * GetClusterContainer(Int_t i=0) const
AliClusterContainer * fCaloClustersCont
Tracks.
Double_t Mag(const AliVParticle &trk)
AliEmcalJet * GetNextAcceptJet()
TObjArray fJetCollArray
jet collection array
void UserExec(Option_t *option)
AliParticleContainer * fTracksCont
Jets.
TH1 *& GetHistogram(Hist_t hist, const Int_t idx=0)
short identifier for the task
AliAnalysisTaskJetsEvshape(const char *name="jets_trg_trd")
AliEmcalList * fOutput
!output list
void FillH2(Hist_t hist, Float_t x, Float_t y, Float_t weight=1., Int_t idx=0)
Definition: External.C:220
void SetMakeGeneralHistograms(Bool_t g)
const Double_t ymin
Definition: AddTaskCFDStar.C:6
Base task in the EMCAL jet framework.
Represent a jet reconstructed using the EMCal jet framework.
Definition: AliEmcalJet.h:44
const char Option_t
Definition: External.C:48
const Double_t zmax
const Double_t pi
bool Bool_t
Definition: External.C:53
virtual void UserExec(Option_t *option)
Definition: External.C:196
Bool_t PrepareEvent()
current run number
TList * OpenFile(const char *fname)
Definition: DrawAnaELoss.C:65