1 /**************************************************************************
2 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
4 * Author: The ALICE Off-line Project. *
5 * Contributors are mentioned in the code where appropriate. *
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 **************************************************************************/
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"
31 #include "AliFilteredTreeEventCuts.h"
35 ClassImp(AliFilteredTreeEventCuts)
37 Int_t AliFilteredTreeEventCuts::fgLastProcessType = -1;
39 //_____________________________________________________________________________
40 AliFilteredTreeEventCuts::AliFilteredTreeEventCuts(const Char_t* name,const Char_t *title) :
41 AliAnalysisCuts(name, title)
42 , fTriggerRequired(kTRUE)
43 , fRecVertexRequired(kTRUE)
44 , fEventProcessType(kInvalidProcess)
45 , fMinNContributors(0)
46 , fMaxNContributors(0)
56 , fRedoTPCVertex(kTRUE)
57 , fUseBeamSpotConstraint(kTRUE)
58 , fEventSelectedRequired(kFALSE)
60 // default constructor
62 // init data members with defaults
66 //_____________________________________________________________________________
67 AliFilteredTreeEventCuts::~AliFilteredTreeEventCuts()
72 //_____________________________________________________________________________
73 void AliFilteredTreeEventCuts::Init()
77 SetRecVertexRequired();
78 SetEventProcessType();
79 SetNContributorsRange();
85 SetUseBeamSpotConstraint();
88 //_____________________________________________________________________________
89 Bool_t AliFilteredTreeEventCuts::AcceptEvent(AliESDEvent *esdEvent,AliMCEvent *mcEvent, const AliESDVertex *vtx)
91 // Check event selection cuts
92 Bool_t retValue=kTRUE;
94 if(!esdEvent) return kFALSE;
95 if(!IsRecVertexRequired()) return kTRUE;
96 if(!vtx) return kFALSE;
97 if(!vtx->GetStatus()) return kFALSE;
100 // check MC event conditions
101 AliHeader* header = mcEvent->Header();
102 if(!header) return kFALSE;
104 // select event type (ND-non diffractive, SD-single diffractive, DD-double diffractive)
105 if(fEventProcessType == kInvalidProcess) {
108 else if(fEventProcessType == kSD || fEventProcessType == kDD) {
109 MCProcessType processType = GetEventProcessType(header);
110 if(processType == kND) retValue=kFALSE;
113 else if(fEventProcessType == GetEventProcessType(header)) {
120 if(vtx->GetZv() < fMinZv) return kFALSE;
121 if(vtx->GetZv() > fMaxZv) return kFALSE;
126 //_____________________________________________________________________________
127 Bool_t AliFilteredTreeEventCuts::AcceptMCEvent(AliMCEvent *mcEvent)
129 // Check event selection cuts
130 if(!mcEvent) return kFALSE;
132 Bool_t retValue=kTRUE;
134 // check MC event conditions
135 AliHeader* header = mcEvent->Header();
136 if(!header) return kFALSE;
138 AliGenEventHeader* genHeader = header->GenEventHeader();
140 AliDebug(AliLog::kError, "Could not retrieve genHeader from Header");
144 genHeader->PrimaryVertex(vtxMC);
146 // select event type (ND-non diffractive, SD-single diffractive, DD-double diffractive)
147 if(fEventProcessType == kInvalidProcess) {
150 if(fEventProcessType == GetEventProcessType(header)) retValue=kTRUE;
151 else retValue=kFALSE;
155 Float_t R = TMath::Sqrt(vtxMC[0]*vtxMC[0]+vtxMC[1]*vtxMC[1]);
156 if(R > fMaxR) return kFALSE;
159 if(vtxMC[2] < fMinZv) return kFALSE;
160 if(vtxMC[2] > fMaxZv) return kFALSE;
165 //_____________________________________________________________________________
166 Long64_t AliFilteredTreeEventCuts::Merge(TCollection* list)
168 // Merge list of objects (needed by PROOF)
175 TIterator* iter = list->MakeIterator();
179 while((obj = iter->Next()) != 0)
181 AliFilteredTreeEventCuts* entry = dynamic_cast<AliFilteredTreeEventCuts*>(obj);
191 //_____________________________________________________________________________
192 AliFilteredTreeEventCuts::MCProcessType AliFilteredTreeEventCuts::GetEventProcessType(AliHeader* aHeader, Bool_t adebug) {
194 // get the process type of the event.
198 // Check for simple headers first
200 AliGenPythiaEventHeader* pythiaGenHeader = dynamic_cast<AliGenPythiaEventHeader*>(aHeader->GenEventHeader());
201 if (pythiaGenHeader) {
202 return GetPythiaEventProcessType(pythiaGenHeader,adebug);
205 AliGenDPMjetEventHeader* dpmJetGenHeader = dynamic_cast<AliGenDPMjetEventHeader*>(aHeader->GenEventHeader());
206 if (dpmJetGenHeader) {
207 return GetDPMjetEventProcessType(dpmJetGenHeader,adebug);
211 // check for cocktail
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;
219 TList* headerList = genCocktailHeader->GetHeaders();
221 return kInvalidProcess;
224 for (Int_t i=0; i<headerList->GetEntries(); i++) {
226 pythiaGenHeader = dynamic_cast<AliGenPythiaEventHeader*>(headerList->At(i));
227 if (pythiaGenHeader) {
228 return GetPythiaEventProcessType(pythiaGenHeader,adebug);
231 dpmJetGenHeader = dynamic_cast<AliGenDPMjetEventHeader*>(headerList->At(i));
232 if (dpmJetGenHeader) {
233 return GetDPMjetEventProcessType(dpmJetGenHeader,adebug);
236 return kInvalidProcess;
239 //____________________________________________________________________
240 AliFilteredTreeEventCuts::MCProcessType AliFilteredTreeEventCuts::GetEventProcessType(AliESDEvent* esd, AliHeader* header, AliStack* stack, DiffTreatment 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
252 MCProcessType mcProcessType = GetEventProcessType(header);
254 if (diffTreatment == kMCFlags)
255 return mcProcessType;
259 Printf("ERROR: AliFilteredTreeEventCuts::GetEventProcessType: diffTreatment != kMCFlags and esd == 0");
260 return kInvalidProcess;
263 Float_t cms = esd->GetESDRun()->GetBeamEnergy();
264 if (esd->GetESDRun()->IsBeamEnergyIsSqrtSHalfGeV())
266 //Printf("cms = %f", cms);
268 if (diffTreatment == kUA5Cuts && mcProcessType == kSD)
270 if (IsHadronLevelSingleDiffractive(stack, cms, 0, 0.05))
273 else if (diffTreatment == kE710Cuts && mcProcessType == kSD)
275 if (IsHadronLevelSingleDiffractive(stack, cms, 2. / cms / cms, 0.05))
278 else if (diffTreatment == kALICEHadronLevel)
280 if (IsHadronLevelSingleDiffractive(stack, cms, 0, 0.05))
284 if (mcProcessType == kSD)
287 return mcProcessType;
291 AliFilteredTreeEventCuts::MCProcessType AliFilteredTreeEventCuts::GetPythiaEventProcessType(AliGenEventHeader* aHeader, Bool_t adebug) {
293 AliGenPythiaEventHeader* pythiaGenHeader = dynamic_cast<AliGenPythiaEventHeader*>(aHeader);
295 if (!pythiaGenHeader) {
296 printf("AliFilteredTreeEventCuts::GetProcessType : Unknown gen Header type). \n");
297 return kInvalidProcess;
301 Int_t pythiaType = pythiaGenHeader->ProcessType();
302 fgLastProcessType = pythiaType;
303 MCProcessType globalType = kInvalidProcess;
307 printf("AliFilteredTreeEventCuts::GetProcessType : Pythia process type found: %d \n",pythiaType);
311 if(pythiaType==92||pythiaType==93){
314 else if(pythiaType==94){
317 //else if(pythiaType != 91){ // also exclude elastic to be sure... CKB??}
325 AliFilteredTreeEventCuts::MCProcessType AliFilteredTreeEventCuts::GetDPMjetEventProcessType(AliGenEventHeader* aHeader, Bool_t adebug) {
327 // get the process type of the event.
330 // can only read pythia headers, either directly or from cocktalil header
331 AliGenDPMjetEventHeader* dpmJetGenHeader = dynamic_cast<AliGenDPMjetEventHeader*>(aHeader);
333 if (!dpmJetGenHeader) {
334 printf("AliFilteredTreeEventCuts::GetDPMjetProcessType : Unknown header type (not DPMjet or). \n");
335 return kInvalidProcess;
338 Int_t dpmJetType = dpmJetGenHeader->ProcessType();
339 fgLastProcessType = dpmJetType;
340 MCProcessType globalType = kInvalidProcess;
344 printf("AliFilteredTreeEventCuts::GetDPMJetProcessType : DPMJet process type found: %d \n",dpmJetType);
348 if (dpmJetType == 1 || dpmJetType == 4) { // explicitly inelastic plus central diffraction
351 else if (dpmJetType==5 || dpmJetType==6) {
354 else if (dpmJetType==7) {
360 //____________________________________________________________________
361 Bool_t AliFilteredTreeEventCuts::IsHadronLevelSingleDiffractive(AliStack* stack, Float_t cms, Float_t xiMin, Float_t xiMax)
364 // return if process is single diffractive on hadron level
366 // xiMax and xiMin cut on M^2/s
368 // Based on code from Martin Poghoysan
371 TParticle* part1 = 0;
372 TParticle* part2 = 0;
374 Double_t smallestY = 1e10;
375 Double_t largestY = -1e10;
377 for (Int_t iParticle = 0; iParticle < stack->GetNprimary(); iParticle++)
379 TParticle* part = stack->Particle(iParticle);
383 Int_t pdg = TMath::Abs(part->GetPdgCode());
385 Int_t child1 = part->GetFirstDaughter();
386 Int_t ist = part->GetStatusCode();
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))
393 Int_t imother = part->GetFirstMother();
396 TParticle *partM = stack->Particle(imother);
397 Int_t pdgM=TMath::Abs(partM->GetPdgCode());
398 if (pdgM==111 || pdgM==3124 || pdgM==3212)
404 // fix for problem with getting mass of particle 3124
406 y = Rapidity(part->Pt(), part->Pz(), part->GetMass());
408 y = Rapidity(part->Pt(), part->Pz(), 1.5195);
423 if (part1 == 0 || part2 == 0)
426 Int_t pdg1 = part1->GetPdgCode();
427 Int_t pdg2 = part2->GetPdgCode();
429 Double_t pt1 = part1->Pt();
430 Double_t pt2 = part2->Pt();
431 Double_t pz1 = part1->Pz();
432 Double_t pz2 = part2->Pz();
434 Double_t y1 = TMath::Abs(Rapidity(pt1, pz1, 0.938));
435 Double_t y2 = TMath::Abs(Rapidity(pt2, pz2, 0.938));
438 if (pdg1 == 2212 && pdg2 == 2212)
445 else if (pdg1 == 2212)
447 else if (pdg2 == 2212)
450 Double_t M02s = 1. - 2 * part1->Energy() / cms;
451 Double_t M12s = 1. - 2 * part2->Energy() / cms;
453 if (arm == 0 && M02s > xiMin && M02s < xiMax)
455 else if (arm == 1 && M12s > xiMin && M12s < xiMax)
461 //____________________________________________________________________
462 Double_t AliFilteredTreeEventCuts::Rapidity(Double_t pt, Double_t pz, Double_t m)
465 // calculates rapidity keeping the sign in case E == pz
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));
473 return TMath::Sign(1.e30,pz);