README Update
[u/mrichter/AliRoot.git] / PWGPP / AliFilteredTreeEventCuts.cxx
CommitLineData
1059f583 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
31#include "AliFilteredTreeEventCuts.h"
32
33using namespace std;
34
35ClassImp(AliFilteredTreeEventCuts)
36
37Int_t AliFilteredTreeEventCuts::fgLastProcessType = -1;
38
39//_____________________________________________________________________________
40AliFilteredTreeEventCuts::AliFilteredTreeEventCuts(const Char_t* name,const Char_t *title) :
41AliAnalysisCuts(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//_____________________________________________________________________________
67AliFilteredTreeEventCuts::~AliFilteredTreeEventCuts()
68{
69 // destructor
70}
71
72//_____________________________________________________________________________
73void AliFilteredTreeEventCuts::Init()
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//_____________________________________________________________________________
89Bool_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->GetZv() < fMinZv) return kFALSE;
121 if(vtx->GetZv() > fMaxZv) return kFALSE;
122
123return retValue;
124}
125
126//_____________________________________________________________________________
127Bool_t AliFilteredTreeEventCuts::AcceptMCEvent(AliMCEvent *mcEvent)
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
162return retValue;
163}
164
165//_____________________________________________________________________________
166Long64_t AliFilteredTreeEventCuts::Merge(TCollection* list)
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
188return count;
189}
190
191//_____________________________________________________________________________
192AliFilteredTreeEventCuts::MCProcessType AliFilteredTreeEventCuts::GetEventProcessType(AliHeader* aHeader, Bool_t adebug) {
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//____________________________________________________________________
240AliFilteredTreeEventCuts::MCProcessType AliFilteredTreeEventCuts::GetEventProcessType(AliESDEvent* esd, AliHeader* header, AliStack* stack, DiffTreatment diffTreatment)
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
291AliFilteredTreeEventCuts::MCProcessType AliFilteredTreeEventCuts::GetPythiaEventProcessType(AliGenEventHeader* aHeader, Bool_t adebug) {
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
325AliFilteredTreeEventCuts::MCProcessType AliFilteredTreeEventCuts::GetDPMjetEventProcessType(AliGenEventHeader* aHeader, Bool_t adebug) {
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//____________________________________________________________________
361Bool_t AliFilteredTreeEventCuts::IsHadronLevelSingleDiffractive(AliStack* stack, Float_t cms, Float_t xiMin, Float_t xiMax)
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//____________________________________________________________________
462Double_t AliFilteredTreeEventCuts::Rapidity(Double_t pt, Double_t pz, Double_t m)
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