AliPhysics  2c6b7ad (2c6b7ad)
AliFilteredTreeEventCuts.cxx
Go to the documentation of this file.
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 #include <iostream>
17 #include <TList.h>
18 
19 #include "AliLog.h"
20 #include "AliESDEvent.h"
21 #include "AliESDVertex.h"
22 #include "AliMCEvent.h"
23 #include "AliHeader.h"
24 #include "AliGenEventHeader.h"
25 #include "AliGenPythiaEventHeader.h"
26 #include "AliGenCocktailEventHeader.h"
27 #include "AliGenDPMjetEventHeader.h"
28 #include "AliStack.h"
29 
30 
32 
33 using namespace std;
34 
36 
37 Int_t AliFilteredTreeEventCuts::fgLastProcessType = -1;
38 
39 //_____________________________________________________________________________
41 AliAnalysisCuts(name, title)
42 , fTriggerRequired(kTRUE)
43 , fRecVertexRequired(kTRUE)
44 , fEventProcessType(kInvalidProcess)
45 , fMinNContributors(0)
46 , fMaxNContributors(0)
47 , fMaxR(0)
48 , fMinZv(0)
49 , fMaxZv(0)
50 , fMeanXv(0)
51 , fMeanYv(0)
52 , fMeanZv(0)
53 , fSigmaMeanXv(0)
54 , fSigmaMeanYv(0)
55 , fSigmaMeanZv(0)
56 , fRedoTPCVertex(kTRUE)
57 , fUseBeamSpotConstraint(kTRUE)
58 , fEventSelectedRequired(kFALSE)
59 {
60  // default constructor
61 
62  // init data members with defaults
63  Init();
64 }
65 
66 //_____________________________________________________________________________
68 {
69  // destructor
70 }
71 
72 //_____________________________________________________________________________
74 {
75  // set default values
76  SetTriggerRequired();
77  SetRecVertexRequired();
78  SetEventProcessType();
79  SetNContributorsRange();
80  SetMaxR();
81  SetZvRange();
82  SetMeanXYZv();
83  SetSigmaMeanXYZv();
84  SetRedoTPCVertex();
85  SetUseBeamSpotConstraint();
86 }
87 
88 //_____________________________________________________________________________
89 Bool_t AliFilteredTreeEventCuts::AcceptEvent(AliESDEvent *esdEvent,AliMCEvent *mcEvent, const AliESDVertex *vtx)
90 {
91  // Check event selection cuts
92  Bool_t retValue=kTRUE;
93 
94  if(!esdEvent) return kFALSE;
95  if(!IsRecVertexRequired()) return kTRUE;
96  if(!vtx) return kFALSE;
97  if(!vtx->GetStatus()) return kFALSE;
98 
99  if(mcEvent) {
100  // check MC event conditions
101  AliHeader* header = mcEvent->Header();
102  if(!header) return kFALSE;
103 
104  // select event type (ND-non diffractive, SD-single diffractive, DD-double diffractive)
105  if(fEventProcessType == kInvalidProcess) {
106  retValue=kTRUE;
107  }
108  else if(fEventProcessType == kSD || fEventProcessType == kDD) {
109  MCProcessType processType = GetEventProcessType(header);
110  if(processType == kND) retValue=kFALSE;
111  else retValue=kTRUE;
112  }
113  else if(fEventProcessType == GetEventProcessType(header)) {
114  retValue=kTRUE;
115  }
116  else
117  retValue=kFALSE;
118  }
119 
120  if(vtx->GetZ() < fMinZv) return kFALSE;
121  if(vtx->GetZ() > fMaxZv) return kFALSE;
122 
123 return retValue;
124 }
125 
126 //_____________________________________________________________________________
128 {
129  // Check event selection cuts
130  if(!mcEvent) return kFALSE;
131 
132  Bool_t retValue=kTRUE;
133 
134  // check MC event conditions
135  AliHeader* header = mcEvent->Header();
136  if(!header) return kFALSE;
137 
138  AliGenEventHeader* genHeader = header->GenEventHeader();
139  if (!genHeader) {
140  AliDebug(AliLog::kError, "Could not retrieve genHeader from Header");
141  return kFALSE;
142  }
143  TArrayF vtxMC(3);
144  genHeader->PrimaryVertex(vtxMC);
145 
146  // select event type (ND-non diffractive, SD-single diffractive, DD-double diffractive)
147  if(fEventProcessType == kInvalidProcess) {
148  retValue=kTRUE;
149  } else {
150  if(fEventProcessType == GetEventProcessType(header)) retValue=kTRUE;
151  else retValue=kFALSE;
152  }
153 
154  /*
155  Float_t R = TMath::Sqrt(vtxMC[0]*vtxMC[0]+vtxMC[1]*vtxMC[1]);
156  if(R > fMaxR) return kFALSE;
157  */
158 
159  if(vtxMC[2] < fMinZv) return kFALSE;
160  if(vtxMC[2] > fMaxZv) return kFALSE;
161 
162 return retValue;
163 }
164 
165 //_____________________________________________________________________________
167 {
168  // Merge list of objects (needed by PROOF)
169  if (!list)
170  return 0;
171 
172  if (list->IsEmpty())
173  return 1;
174 
175  TIterator* iter = list->MakeIterator();
176  TObject* obj = 0;
177 
178  Int_t count=0;
179  while((obj = iter->Next()) != 0)
180  {
181  AliFilteredTreeEventCuts* entry = dynamic_cast<AliFilteredTreeEventCuts*>(obj);
182  if (entry == 0)
183  continue;
184 
185  count++;
186  }
187 
188 return count;
189 }
190 
191 //_____________________________________________________________________________
193  //
194  // get the process type of the event.
195  //
196 
197 
198  // Check for simple headers first
199 
200  AliGenPythiaEventHeader* pythiaGenHeader = dynamic_cast<AliGenPythiaEventHeader*>(aHeader->GenEventHeader());
201  if (pythiaGenHeader) {
202  return GetPythiaEventProcessType(pythiaGenHeader,adebug);
203  }
204 
205  AliGenDPMjetEventHeader* dpmJetGenHeader = dynamic_cast<AliGenDPMjetEventHeader*>(aHeader->GenEventHeader());
206  if (dpmJetGenHeader) {
207  return GetDPMjetEventProcessType(dpmJetGenHeader,adebug);
208  }
209 
210 
211  // check for cocktail
212 
213  AliGenCocktailEventHeader* genCocktailHeader = dynamic_cast<AliGenCocktailEventHeader*>(aHeader->GenEventHeader());
214  if (!genCocktailHeader) {
215  printf("AliFilteredTreeEventCuts::GetProcessType : Unknown header type (not Pythia or Cocktail). \n");
216  return kInvalidProcess;
217  }
218 
219  TList* headerList = genCocktailHeader->GetHeaders();
220  if (!headerList) {
221  return kInvalidProcess;
222  }
223 
224  for (Int_t i=0; i<headerList->GetEntries(); i++) {
225 
226  pythiaGenHeader = dynamic_cast<AliGenPythiaEventHeader*>(headerList->At(i));
227  if (pythiaGenHeader) {
228  return GetPythiaEventProcessType(pythiaGenHeader,adebug);
229  }
230 
231  dpmJetGenHeader = dynamic_cast<AliGenDPMjetEventHeader*>(headerList->At(i));
232  if (dpmJetGenHeader) {
233  return GetDPMjetEventProcessType(dpmJetGenHeader,adebug);
234  }
235  }
236  return kInvalidProcess;
237 }
238 
239 //____________________________________________________________________
241 {
242  //
243  // get process type
244  //
245  // diffTreatment:
246  // kMCFlags: use MC flags
247  // kUA5Cuts: SD events are those that have the MC SD flag and fulfill M^2/s < 0.05; DD from MC flag; Remainder is ND
248  // kE710Cuts: SD events are those that have the MC SD flag and fulfill 2 < M^2 < 0.05s; DD from MC flag; Remainder is ND
249  // kALICEHadronLevel: SD events are those that fulfill M^2/s < 0.05; DD from MC flag; Remainder is ND
250  //
251 
252  MCProcessType mcProcessType = GetEventProcessType(header);
253 
254  if (diffTreatment == kMCFlags)
255  return mcProcessType;
256 
257  if (!esd)
258  {
259  Printf("ERROR: AliFilteredTreeEventCuts::GetEventProcessType: diffTreatment != kMCFlags and esd == 0");
260  return kInvalidProcess;
261  }
262 
263  Float_t cms = esd->GetESDRun()->GetBeamEnergy();
264  if (esd->GetESDRun()->IsBeamEnergyIsSqrtSHalfGeV())
265  cms *= 2;
266  //Printf("cms = %f", cms);
267 
268  if (diffTreatment == kUA5Cuts && mcProcessType == kSD)
269  {
270  if (IsHadronLevelSingleDiffractive(stack, cms, 0, 0.05))
271  return kSD;
272  }
273  else if (diffTreatment == kE710Cuts && mcProcessType == kSD)
274  {
275  if (IsHadronLevelSingleDiffractive(stack, cms, 2. / cms / cms, 0.05))
276  return kSD;
277  }
278  else if (diffTreatment == kALICEHadronLevel)
279  {
280  if (IsHadronLevelSingleDiffractive(stack, cms, 0, 0.05))
281  return kSD;
282  }
283 
284  if (mcProcessType == kSD)
285  return kND;
286 
287  return mcProcessType;
288 }
289 
290 
292 
293  AliGenPythiaEventHeader* pythiaGenHeader = dynamic_cast<AliGenPythiaEventHeader*>(aHeader);
294 
295  if (!pythiaGenHeader) {
296  printf("AliFilteredTreeEventCuts::GetProcessType : Unknown gen Header type). \n");
297  return kInvalidProcess;
298  }
299 
300 
301  Int_t pythiaType = pythiaGenHeader->ProcessType();
302  fgLastProcessType = pythiaType;
303  MCProcessType globalType = kInvalidProcess;
304 
305 
306  if (adebug) {
307  printf("AliFilteredTreeEventCuts::GetProcessType : Pythia process type found: %d \n",pythiaType);
308  }
309 
310 
311  if(pythiaType==92||pythiaType==93){
312  globalType = kSD;
313  }
314  else if(pythiaType==94){
315  globalType = kDD;
316  }
317  //else if(pythiaType != 91){ // also exclude elastic to be sure... CKB??}
318  else {
319  globalType = kND;
320  }
321  return globalType;
322 }
323 
324 
326  //
327  // get the process type of the event.
328  //
329 
330  // can only read pythia headers, either directly or from cocktalil header
331  AliGenDPMjetEventHeader* dpmJetGenHeader = dynamic_cast<AliGenDPMjetEventHeader*>(aHeader);
332 
333  if (!dpmJetGenHeader) {
334  printf("AliFilteredTreeEventCuts::GetDPMjetProcessType : Unknown header type (not DPMjet or). \n");
335  return kInvalidProcess;
336  }
337 
338  Int_t dpmJetType = dpmJetGenHeader->ProcessType();
339  fgLastProcessType = dpmJetType;
340  MCProcessType globalType = kInvalidProcess;
341 
342 
343  if (adebug) {
344  printf("AliFilteredTreeEventCuts::GetDPMJetProcessType : DPMJet process type found: %d \n",dpmJetType);
345  }
346 
347 
348  if (dpmJetType == 1 || dpmJetType == 4) { // explicitly inelastic plus central diffraction
349  globalType = kND;
350  }
351  else if (dpmJetType==5 || dpmJetType==6) {
352  globalType = kSD;
353  }
354  else if (dpmJetType==7) {
355  globalType = kDD;
356  }
357  return globalType;
358 }
359 
360 //____________________________________________________________________
362 {
363  //
364  // return if process is single diffractive on hadron level
365  //
366  // xiMax and xiMin cut on M^2/s
367  //
368  // Based on code from Martin Poghoysan
369  //
370 
371  TParticle* part1 = 0;
372  TParticle* part2 = 0;
373 
374  Double_t smallestY = 1e10;
375  Double_t largestY = -1e10;
376 
377  for (Int_t iParticle = 0; iParticle < stack->GetNprimary(); iParticle++)
378  {
379  TParticle* part = stack->Particle(iParticle);
380  if (!part)
381  continue;
382 
383  Int_t pdg = TMath::Abs(part->GetPdgCode());
384 
385  Int_t child1 = part->GetFirstDaughter();
386  Int_t ist = part->GetStatusCode();
387 
388  Int_t mfl = Int_t (pdg / TMath::Power(10, Int_t(TMath::Log10(pdg))));
389  if (child1 > -1 || ist != 1)
390  mfl = 0; // select final state charm and beauty
391  if (!(stack->IsPhysicalPrimary(iParticle) || pdg == 111 || pdg == 3212 || pdg==3124 || mfl >= 4))
392  continue;
393  Int_t imother = part->GetFirstMother();
394  if (imother>0)
395  {
396  TParticle *partM = stack->Particle(imother);
397  Int_t pdgM=TMath::Abs(partM->GetPdgCode());
398  if (pdgM==111 || pdgM==3124 || pdgM==3212)
399  continue;
400  }
401 
402  Double_t y = 0;
403 
404  // fix for problem with getting mass of particle 3124
405  if (pdg != 3124)
406  y = Rapidity(part->Pt(), part->Pz(), part->GetMass());
407  else
408  y = Rapidity(part->Pt(), part->Pz(), 1.5195);
409 
410  if (y < smallestY)
411  {
412  smallestY = y;
413  part1 = part;
414  }
415 
416  if (y > largestY)
417  {
418  largestY = y;
419  part2 = part;
420  }
421  }
422 
423  if (part1 == 0 || part2 == 0)
424  return kFALSE;
425 
426  Int_t pdg1 = part1->GetPdgCode();
427  Int_t pdg2 = part2->GetPdgCode();
428 
429  Double_t pt1 = part1->Pt();
430  Double_t pt2 = part2->Pt();
431  Double_t pz1 = part1->Pz();
432  Double_t pz2 = part2->Pz();
433 
434  Double_t y1 = TMath::Abs(Rapidity(pt1, pz1, 0.938));
435  Double_t y2 = TMath::Abs(Rapidity(pt2, pz2, 0.938));
436 
437  Int_t arm = -99999;
438  if (pdg1 == 2212 && pdg2 == 2212)
439  {
440  if (y1 > y2)
441  arm = 0;
442  else
443  arm = 1;
444  }
445  else if (pdg1 == 2212)
446  arm = 0;
447  else if (pdg2 == 2212)
448  arm = 1;
449 
450  Double_t M02s = 1. - 2 * part1->Energy() / cms;
451  Double_t M12s = 1. - 2 * part2->Energy() / cms;
452 
453  if (arm == 0 && M02s > xiMin && M02s < xiMax)
454  return kTRUE;
455  else if (arm == 1 && M12s > xiMin && M12s < xiMax)
456  return kTRUE;
457 
458  return kFALSE;
459 }
460 
461 //____________________________________________________________________
463 {
464  //
465  // calculates rapidity keeping the sign in case E == pz
466  //
467 
468  Double_t energy = TMath::Sqrt(pt*pt+pz*pz+m*m);
469  if (energy != TMath::Abs(pz))
470  return 0.5*TMath::Log((energy+pz)/(energy-pz));
471 
472  Printf("W- mt=0");
473  return TMath::Sign(1.e30,pz);
474 }
475 
Int_t pdg
double Double_t
Definition: External.C:58
const char * title
Definition: MakeQAPdf.C:27
long long Long64_t
Definition: External.C:43
energy
Definition: HFPtSpectrum.C:44
char Char_t
Definition: External.C:18
static MCProcessType GetPythiaEventProcessType(AliGenEventHeader *aHeader, Bool_t adebug=kFALSE)
static Bool_t IsHadronLevelSingleDiffractive(AliStack *stack, Float_t cms, Float_t xiMin, Float_t xiMax)
AliStack * stack
int Int_t
Definition: External.C:63
float Float_t
Definition: External.C:68
static MCProcessType GetDPMjetEventProcessType(AliGenEventHeader *aHeader, Bool_t adebug=kFALSE)
virtual Long64_t Merge(TCollection *list)
Bool_t AcceptMCEvent(AliMCEvent *mcEvent=0)
bool Bool_t
Definition: External.C:53
Bool_t AcceptEvent(AliESDEvent *event=0, AliMCEvent *mcEvent=0, const AliESDVertex *vtx=0)
static Double_t Rapidity(Double_t pt, Double_t pz, Double_t m)