[u/mrichter/AliRoot.git] / PWGPP / AliFilteredTreeEventCuts.cxx

/**************************************************************************
* Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
*                                                                        *
* Author: The ALICE Off-line Project.                                    *
* Contributors are mentioned in the code where appropriate.              *
*                                                                        *
* Permission to use, copy, modify and distribute this software and its   *
* documentation strictly for non-commercial purposes is hereby granted   *
* without fee, provided that the above copyright notice appears in all   *
* copies and that both the copyright notice and this permission notice   *
* appear in the supporting documentation. The authors make no claims     *
* about the suitability of this software for any purpose. It is          *
* provided "as is" without express or implied warranty.                  *
**************************************************************************/

#include <iostream>
#include <TList.h>

#include "AliLog.h"
#include "AliESDEvent.h"
#include "AliESDVertex.h"
#include "AliMCEvent.h"
#include "AliHeader.h"
#include "AliGenEventHeader.h"
#include "AliGenPythiaEventHeader.h"
#include "AliGenCocktailEventHeader.h"
#include "AliGenDPMjetEventHeader.h"
#include "AliStack.h"


#include "AliFilteredTreeEventCuts.h"

using namespace std;

ClassImp(AliFilteredTreeEventCuts)

Int_t AliFilteredTreeEventCuts::fgLastProcessType = -1;

//_____________________________________________________________________________
AliFilteredTreeEventCuts::AliFilteredTreeEventCuts(const Char_t* name,const Char_t *title) : 
AliAnalysisCuts(name, title)
, fTriggerRequired(kTRUE)
, fRecVertexRequired(kTRUE)
, fEventProcessType(kInvalidProcess)
, fMinNContributors(0)
, fMaxNContributors(0)
, fMaxR(0)
, fMinZv(0)
, fMaxZv(0)
, fMeanXv(0)
, fMeanYv(0)
, fMeanZv(0)
, fSigmaMeanXv(0)
, fSigmaMeanYv(0)
, fSigmaMeanZv(0)
, fRedoTPCVertex(kTRUE)
, fUseBeamSpotConstraint(kTRUE)
, fEventSelectedRequired(kFALSE)
{
  // default constructor 
  
  // init data members with defaults
  Init();
}

//_____________________________________________________________________________
AliFilteredTreeEventCuts::~AliFilteredTreeEventCuts()  
{
  // destructor
}

//_____________________________________________________________________________
void AliFilteredTreeEventCuts::Init()  
{
  // set default values
  SetTriggerRequired();
  SetRecVertexRequired();
  SetEventProcessType();
  SetNContributorsRange();
  SetMaxR();
  SetZvRange();
  SetMeanXYZv();
  SetSigmaMeanXYZv();
  SetRedoTPCVertex();
  SetUseBeamSpotConstraint();
}

//_____________________________________________________________________________
Bool_t AliFilteredTreeEventCuts::AcceptEvent(AliESDEvent *esdEvent,AliMCEvent *mcEvent, const AliESDVertex *vtx)
{
  // Check event selection cuts
  Bool_t retValue=kTRUE;

  if(!esdEvent) return kFALSE;
  if(!IsRecVertexRequired()) return kTRUE;
  if(!vtx) return kFALSE;
  if(!vtx->GetStatus()) return kFALSE;

  if(mcEvent) {
   // check MC event conditions
   AliHeader* header = mcEvent->Header();
   if(!header) return kFALSE;
  
    // select event type (ND-non diffractive, SD-single diffractive, DD-double diffractive)
    if(fEventProcessType == kInvalidProcess) { 
      retValue=kTRUE;
    } 
    else if(fEventProcessType == kSD || fEventProcessType == kDD) {
      MCProcessType processType = GetEventProcessType(header);
      if(processType == kND) retValue=kFALSE;
      else retValue=kTRUE;
    }
    else if(fEventProcessType == GetEventProcessType(header)) { 
      retValue=kTRUE;
    }
    else 
      retValue=kFALSE;
  }

  if(vtx->GetZ() < fMinZv) return kFALSE; 
  if(vtx->GetZ() > fMaxZv) return kFALSE; 

return retValue;  
}

//_____________________________________________________________________________
Bool_t AliFilteredTreeEventCuts::AcceptMCEvent(AliMCEvent *mcEvent)
{
  // Check event selection cuts
  if(!mcEvent) return kFALSE;

  Bool_t retValue=kTRUE;

  // check MC event conditions
  AliHeader* header = mcEvent->Header();
  if(!header) return kFALSE;

  AliGenEventHeader* genHeader = header->GenEventHeader();
  if (!genHeader) {
    AliDebug(AliLog::kError, "Could not retrieve genHeader from Header");
    return kFALSE;
  }
  TArrayF vtxMC(3);
  genHeader->PrimaryVertex(vtxMC);
  
  // select event type (ND-non diffractive, SD-single diffractive, DD-double diffractive)
  if(fEventProcessType == kInvalidProcess) { 
     retValue=kTRUE;
  } else {
     if(fEventProcessType == GetEventProcessType(header)) retValue=kTRUE;
     else retValue=kFALSE;
  }

  /*
  Float_t R = TMath::Sqrt(vtxMC[0]*vtxMC[0]+vtxMC[1]*vtxMC[1]);
  if(R > fMaxR) return kFALSE; 
  */

  if(vtxMC[2] < fMinZv) return kFALSE; 
  if(vtxMC[2] > fMaxZv) return kFALSE; 

return retValue;  
}

//_____________________________________________________________________________
Long64_t AliFilteredTreeEventCuts::Merge(TCollection* list) 
{
  // Merge list of objects (needed by PROOF)
  if (!list)
  return 0;

  if (list->IsEmpty())
  return 1;

  TIterator* iter = list->MakeIterator();
  TObject* obj = 0;

  Int_t count=0;
  while((obj = iter->Next()) != 0) 
  {
    AliFilteredTreeEventCuts* entry = dynamic_cast<AliFilteredTreeEventCuts*>(obj);
    if (entry == 0)  
      continue; 

  count++;
  }

return count;
}

//_____________________________________________________________________________
AliFilteredTreeEventCuts::MCProcessType AliFilteredTreeEventCuts::GetEventProcessType(AliHeader* aHeader, Bool_t adebug) {
  //
  // get the process type of the event.
  //


  // Check for simple headers first

  AliGenPythiaEventHeader* pythiaGenHeader = dynamic_cast<AliGenPythiaEventHeader*>(aHeader->GenEventHeader());
  if (pythiaGenHeader) {
    return GetPythiaEventProcessType(pythiaGenHeader,adebug);
  }

  AliGenDPMjetEventHeader* dpmJetGenHeader = dynamic_cast<AliGenDPMjetEventHeader*>(aHeader->GenEventHeader());
  if (dpmJetGenHeader) {
    return GetDPMjetEventProcessType(dpmJetGenHeader,adebug);
  }
  

  // check for cocktail

  AliGenCocktailEventHeader* genCocktailHeader = dynamic_cast<AliGenCocktailEventHeader*>(aHeader->GenEventHeader());
  if (!genCocktailHeader) {
    printf("AliFilteredTreeEventCuts::GetProcessType : Unknown header type (not Pythia or Cocktail). \n");
    return kInvalidProcess;
  }

  TList* headerList = genCocktailHeader->GetHeaders();
  if (!headerList) {
    return kInvalidProcess;
  }

  for (Int_t i=0; i<headerList->GetEntries(); i++) {

    pythiaGenHeader = dynamic_cast<AliGenPythiaEventHeader*>(headerList->At(i));
    if (pythiaGenHeader) {
      return GetPythiaEventProcessType(pythiaGenHeader,adebug);
    }

    dpmJetGenHeader = dynamic_cast<AliGenDPMjetEventHeader*>(headerList->At(i));
    if (dpmJetGenHeader) {
      return GetDPMjetEventProcessType(dpmJetGenHeader,adebug);
    }
  }
  return kInvalidProcess;
}

//____________________________________________________________________
AliFilteredTreeEventCuts::MCProcessType AliFilteredTreeEventCuts::GetEventProcessType(AliESDEvent* esd, AliHeader* header, AliStack* stack, DiffTreatment diffTreatment)
{
  // 
  // get process type
  //
  // diffTreatment:
  //   kMCFlags: use MC flags
  //   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
  //   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
  //   kALICEHadronLevel: SD events are those that fulfill M^2/s < 0.05; DD from MC flag; Remainder is ND
  //

  MCProcessType mcProcessType = GetEventProcessType(header);
  
  if (diffTreatment == kMCFlags)
    return mcProcessType;
    
  if (!esd)
  {
    Printf("ERROR: AliFilteredTreeEventCuts::GetEventProcessType: diffTreatment != kMCFlags and esd == 0");
    return kInvalidProcess;
  }
    
  Float_t cms = esd->GetESDRun()->GetBeamEnergy();
  if (esd->GetESDRun()->IsBeamEnergyIsSqrtSHalfGeV())
    cms *= 2;
  //Printf("cms = %f", cms);

  if (diffTreatment == kUA5Cuts && mcProcessType == kSD)
  {
    if (IsHadronLevelSingleDiffractive(stack, cms, 0, 0.05))
      return kSD;
  }
  else if (diffTreatment == kE710Cuts && mcProcessType == kSD)
  {
    if (IsHadronLevelSingleDiffractive(stack, cms, 2. / cms / cms, 0.05))
      return kSD;
  }
  else if (diffTreatment == kALICEHadronLevel)
  {
    if (IsHadronLevelSingleDiffractive(stack, cms, 0, 0.05))
      return kSD;
  }
  
  if (mcProcessType == kSD)
    return kND;
    
  return mcProcessType;
}


AliFilteredTreeEventCuts::MCProcessType AliFilteredTreeEventCuts::GetPythiaEventProcessType(AliGenEventHeader* aHeader, Bool_t adebug) {

  AliGenPythiaEventHeader* pythiaGenHeader = dynamic_cast<AliGenPythiaEventHeader*>(aHeader);

  if (!pythiaGenHeader) {
    printf("AliFilteredTreeEventCuts::GetProcessType : Unknown gen Header type). \n");
    return kInvalidProcess;
  }


  Int_t pythiaType = pythiaGenHeader->ProcessType();
  fgLastProcessType = pythiaType;
  MCProcessType globalType = kInvalidProcess;  


  if (adebug) {
    printf("AliFilteredTreeEventCuts::GetProcessType : Pythia process type found: %d \n",pythiaType);
  }


  if(pythiaType==92||pythiaType==93){
    globalType = kSD;
  }
  else if(pythiaType==94){
    globalType = kDD;
  }
  //else if(pythiaType != 91){ // also exclude elastic to be sure... CKB??}
  else {
    globalType = kND;
  }
  return globalType;
}


AliFilteredTreeEventCuts::MCProcessType AliFilteredTreeEventCuts::GetDPMjetEventProcessType(AliGenEventHeader* aHeader, Bool_t adebug) {
  //
  // get the process type of the event.
  //

  // can only read pythia headers, either directly or from cocktalil header
  AliGenDPMjetEventHeader* dpmJetGenHeader = dynamic_cast<AliGenDPMjetEventHeader*>(aHeader);

  if (!dpmJetGenHeader) {
    printf("AliFilteredTreeEventCuts::GetDPMjetProcessType : Unknown header type (not DPMjet or). \n");
    return kInvalidProcess;
  }

  Int_t dpmJetType = dpmJetGenHeader->ProcessType();
  fgLastProcessType = dpmJetType;
  MCProcessType globalType = kInvalidProcess;  


  if (adebug) {
    printf("AliFilteredTreeEventCuts::GetDPMJetProcessType : DPMJet process type found: %d \n",dpmJetType);
  }


  if (dpmJetType == 1 || dpmJetType == 4) { // explicitly inelastic plus central diffraction
    globalType = kND;
  }  
  else if (dpmJetType==5 || dpmJetType==6) {
    globalType = kSD;
  }
  else if (dpmJetType==7) {
    globalType = kDD;
  }
  return globalType;
}

//____________________________________________________________________
Bool_t AliFilteredTreeEventCuts::IsHadronLevelSingleDiffractive(AliStack* stack, Float_t cms, Float_t xiMin, Float_t xiMax)
{
  //
  // return if process is single diffractive on hadron level
  // 
  // xiMax and xiMin cut on M^2/s
  //
  // Based on code from Martin Poghoysan
  //
  
  TParticle* part1 = 0;
  TParticle* part2 = 0;
  
  Double_t smallestY = 1e10;
  Double_t largestY  = -1e10;
  
  for (Int_t iParticle = 0; iParticle < stack->GetNprimary(); iParticle++)
  {
    TParticle* part = stack->Particle(iParticle);
    if (!part)
      continue;

    Int_t pdg = TMath::Abs(part->GetPdgCode());

    Int_t child1 = part->GetFirstDaughter();
    Int_t ist    = part->GetStatusCode();

    Int_t mfl  = Int_t (pdg / TMath::Power(10, Int_t(TMath::Log10(pdg))));
    if (child1 > -1 || ist != 1)
      mfl = 0; // select final state charm and beauty
    if (!(stack->IsPhysicalPrimary(iParticle) || pdg == 111 || pdg == 3212 || pdg==3124 || mfl >= 4)) 
      continue;
    Int_t imother = part->GetFirstMother();
    if (imother>0)
    {
      TParticle *partM = stack->Particle(imother);
      Int_t pdgM=TMath::Abs(partM->GetPdgCode());
      if (pdgM==111 || pdgM==3124 || pdgM==3212) 
        continue;
    }
    
    Double_t y = 0;

    // fix for problem with getting mass of particle 3124
    if (pdg != 3124)
      y = Rapidity(part->Pt(), part->Pz(), part->GetMass());
    else 
      y = Rapidity(part->Pt(), part->Pz(), 1.5195);
      
    if (y < smallestY)
    {
      smallestY = y;
      part1 = part;
    }
    
    if (y > largestY)
    {
      largestY = y;
      part2 = part;
    }
  }
  
  if (part1 == 0 || part2 == 0)
    return kFALSE;

  Int_t pdg1 = part1->GetPdgCode();
  Int_t pdg2 = part2->GetPdgCode();

  Double_t pt1 = part1->Pt();
  Double_t pt2 = part2->Pt();
  Double_t pz1 = part1->Pz();
  Double_t pz2 = part2->Pz();
  
  Double_t y1 = TMath::Abs(Rapidity(pt1, pz1, 0.938));
  Double_t y2 = TMath::Abs(Rapidity(pt2, pz2, 0.938));
  
  Int_t arm = -99999;
  if (pdg1 == 2212 && pdg2 == 2212)
  {
    if (y1 > y2) 
      arm = 0;
    else
      arm = 1;
  }
  else if (pdg1 == 2212) 
    arm = 0;
  else if (pdg2 == 2212) 
    arm = 1;

  Double_t M02s = 1. - 2 * part1->Energy() / cms;
  Double_t M12s = 1. - 2 * part2->Energy() / cms;

  if (arm == 0 && M02s > xiMin && M02s < xiMax)
    return kTRUE;
  else if (arm == 1 && M12s > xiMin && M12s < xiMax)
    return kTRUE;

  return kFALSE;
}

//____________________________________________________________________
Double_t AliFilteredTreeEventCuts::Rapidity(Double_t pt, Double_t pz, Double_t m)
{
  //
  // calculates rapidity keeping the sign in case E == pz
  //

  Double_t energy = TMath::Sqrt(pt*pt+pz*pz+m*m);
  if (energy != TMath::Abs(pz))
    return 0.5*TMath::Log((energy+pz)/(energy-pz));

  Printf("W- mt=0");
  return TMath::Sign(1.e30,pz);
}
Commit	Line	Data
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
	33	using namespace std;
	34
	35	ClassImp(AliFilteredTreeEventCuts)
	36
	37	Int_t AliFilteredTreeEventCuts::fgLastProcessType = -1;
	38
	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)
	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	//_____________________________________________________________________________
67	AliFilteredTreeEventCuts::~AliFilteredTreeEventCuts()
68	{
69	// destructor
70	}
71
72	//_____________________________________________________________________________
73	void 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	//_____________________________________________________________________________
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
e690d4d0	120	if(vtx->GetZ() < fMinZv) return kFALSE;
e690d4d0	121	if(vtx->GetZ() > fMaxZv) return kFALSE;
1059f583	122
	123	return retValue;
	124	}
	125
	126	//_____________________________________________________________________________
	127	Bool_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
	162	return retValue;
	163	}
	164
	165	//_____________________________________________________________________________
	166	Long64_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
188	return count;
189	}
190
191	//_____________________________________________________________________________
192	AliFilteredTreeEventCuts::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	//____________________________________________________________________
240	AliFilteredTreeEventCuts::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
291	AliFilteredTreeEventCuts::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
325	AliFilteredTreeEventCuts::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	//____________________________________________________________________
361	Bool_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	//____________________________________________________________________
462	Double_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(ptpt+pzpz+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