Stupid bug fix in new superlight mode (from Zurich airport)
[u/mrichter/AliRoot.git] / PWGPP / AliFilteredTreeEventCuts.cxx
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
120   if(vtx->GetZ() < fMinZv) return kFALSE; 
121   if(vtx->GetZ() > fMaxZv) return kFALSE; 
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(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