Changes to compile with Root6 on macosx64
[u/mrichter/AliRoot.git] / PWGLF / SPECTRA / PiKaPr / TPCTOF / AliAnalysisCombinedHadronSpectra.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, proviyaded 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 purapose. It is         *
13  * provided "as is" without express or implied warranty.                  *
14  **************************************************************************/
15
16 ///////////////////////////////////////////////////////////////////////////
17 //                                                                       //
18 //                                                                       //
19 // Analysis for identified charged hadron spectra.                       //
20 //                                                                       //
21 //                                                                       //
22 ///////////////////////////////////////////////////////////////////////////
23
24 #include "Riostream.h"
25 #include "TChain.h"
26 #include "TTree.h"
27 #include "TH1F.h"
28 #include "TH2D.h"
29 #include "TH3F.h"
30 #include "TList.h"
31 #include "TMath.h"
32 #include "TCanvas.h"
33 #include "TObjArray.h"
34 #include "TF1.h"
35 #include "TFile.h"
36
37 #include "AliAnalysisTaskSE.h"
38 #include "AliAnalysisManager.h"
39
40 #include "AliHeader.h"
41 #include "AliGenPythiaEventHeader.h"
42 #include "AliGenCocktailEventHeader.h"
43
44 #include "AliPID.h"
45 #include "AliESDtrackCuts.h"
46 #include "AliESDVertex.h"
47 #include "AliESDEvent.h"
48 #include "AliESDInputHandler.h"
49 #include "AliESDtrack.h"
50 #include "AliESDpid.h"
51 #include "AliCentrality.h"
52 #include "AliESDUtils.h"
53 #include "AliMultiplicity.h"
54
55 #include "AliMCEventHandler.h"
56 #include "AliMCEvent.h"
57 #include "AliStack.h"
58
59 #include "AliLog.h"
60
61 #include "AliAnalysisCombinedHadronSpectra.h"
62
63
64 ClassImp(AliAnalysisCombinedHadronSpectra)
65
66 //________________________________________________________________________
67 AliAnalysisCombinedHadronSpectra::AliAnalysisCombinedHadronSpectra() 
68   : AliAnalysisTaskSE("TaskChargedHadron"), fESD(0), fListHist(0), fESDtrackCuts(0),fESDTrackCutsMult(0),fESDpid(0),
69     fMCtrue(0),
70     fOnlyQA(0),
71     fUseHBTmultiplicity(0),
72     fUseTPConlyTracks(0),
73     fSaveMotherPDG(0),
74   fSmallTHnSparse(0),
75   fTPCnSigmaCutLow(0),
76   fTPCnSigmaCutHigh(0),
77   fRapidityCutLow(0),
78   fRapidityCutHigh(0),
79   fEvenDCAbinning(0),
80     fAlephParameters(),
81     fHistRealTracks(0),
82     fHistMCparticles(0),
83     fHistPidQA(0),
84     fHistMult(0),
85     fHistCentrality(0)
86 {
87   // default Constructor
88   /* fast compilation test
89      gSystem->Load("libANALYSIS");
90      gSystem->Load("libANALYSISalice");
91      .L /d/alice09/akalweit/train/trunk/akalweit_hadronSpectra/AliAnalysisCombinedHadronSpectra.cxx++
92    */
93 }
94
95
96 //________________________________________________________________________
97 AliAnalysisCombinedHadronSpectra::AliAnalysisCombinedHadronSpectra(const char *name) 
98   : AliAnalysisTaskSE(name), fESD(0), fListHist(0), fESDtrackCuts(0),fESDTrackCutsMult(0),fESDpid(0),
99     fMCtrue(0),
100     fOnlyQA(0),
101     fUseHBTmultiplicity(0),
102     fUseTPConlyTracks(0),
103     fSaveMotherPDG(0),
104   fSmallTHnSparse(0),
105   fTPCnSigmaCutLow(0),
106   fTPCnSigmaCutHigh(0),
107   fRapidityCutLow(0),
108   fRapidityCutHigh(0),
109   fEvenDCAbinning(0),
110     fAlephParameters(),
111     fHistRealTracks(0),
112     fHistMCparticles(0),
113     fHistPidQA(0),
114     fHistMult(0),
115     fHistCentrality(0)
116     {
117   //
118   // standard constructur which should be used
119   //
120   Printf("*** CONSTRUCTOR CALLED ****");
121   //
122   fMCtrue = kTRUE; 
123   fOnlyQA = kFALSE;
124   fUseHBTmultiplicity = kTRUE;
125   fUseTPConlyTracks = kFALSE;
126
127   fSmallTHnSparse = kFALSE;
128   fTPCnSigmaCutLow = -3.;
129   fTPCnSigmaCutHigh = 3.;
130   fRapidityCutLow = -0.2;
131   fRapidityCutHigh = 0.2;
132   fEvenDCAbinning = kFALSE;
133
134   /* real */
135   fAlephParameters[0] = 0.0283086;
136   fAlephParameters[1] = 2.63394e+01;
137   fAlephParameters[2] = 5.04114e-11;
138   fAlephParameters[3] = 2.12543e+00;
139   fAlephParameters[4] = 4.88663e+00;
140   //
141   // initialize PID object
142   //
143   //fESDpid = new AliESDpid();
144   //
145   // create track cuts
146   //
147   fESDtrackCuts = new AliESDtrackCuts("AliESDtrackCuts","AliESDtrackCuts");
148   //
149   //Initialize();
150   // Output slot #0 writes into a TList container
151   DefineOutput(1, TList::Class());
152
153 }
154
155
156 //________________________________________________________________________
157 void AliAnalysisCombinedHadronSpectra::Initialize()
158 {
159   //
160   // updating parameters in case of changes
161   //
162   // 1. pid parameters
163   //
164   //fESDpid->GetTPCResponse().SetBetheBlochParameters(fAlephParameters[0],fAlephParameters[1],fAlephParameters[2],fAlephParameters[3],fAlephParameters[4]);
165   //
166   // 2. track cuts
167   //
168   /*
169   fESDtrackCuts->SetMaxCovDiagonalElements(2, 2, 0.5, 0.5, 2);  // BEWARE STANDARD VALUES ARE: 2, 2, 0.5, 0.5, 2
170   fESDtrackCuts->SetMaxNsigmaToVertex(3);
171   fESDtrackCuts->SetRequireSigmaToVertex(kTRUE);
172   fESDtrackCuts->SetAcceptKinkDaughters(kFALSE);
173   fESDtrackCuts->SetMinNClustersTPC(70);
174   fESDtrackCuts->SetMaxChi2PerClusterTPC(4);
175   fESDtrackCuts->SetMaxDCAToVertexXY(3);
176   fESDtrackCuts->SetMaxDCAToVertexZ(3);
177   fESDtrackCuts->SetRequireTPCRefit(kTRUE);
178   fESDtrackCuts->SetRequireITSRefit(kTRUE);
179   fESDtrackCuts->SetClusterRequirementITS(AliESDtrackCuts::kSPD, AliESDtrackCuts::kAny); //TEMPORARY <-> REMOVE
180   fESDtrackCuts->SetMinNClustersITS(3);
181   */
182   //fESDtrackCuts = AliESDtrackCuts::GetStandardITSTPCTrackCuts2010(kTRUE); // kTRUE = sel. primaries --> patch for the moment, do TFractionFitter later
183
184
185   if (!fUseTPConlyTracks) {
186           fESDtrackCuts = AliESDtrackCuts::GetStandardITSTPCTrackCuts2010(kFALSE);
187           fESDtrackCuts->SetMaxDCAToVertexXY(3);
188           fESDtrackCuts->SetMaxDCAToVertexZ(2);
189           fESDtrackCuts->SetEtaRange(-0.9,0.9);
190   }
191   else {
192           //fESDtrackCuts = AliESDtrackCuts::GetStandardTPCOnlyTrackCuts();
193           fESDtrackCuts->SetMinNClustersTPC(70);
194           fESDtrackCuts->SetMaxChi2PerClusterTPC(4);
195           fESDtrackCuts->SetAcceptKinkDaughters(kFALSE);
196           fESDtrackCuts->SetRequireTPCRefit(kFALSE);
197
198           fESDtrackCuts->SetMaxDCAToVertexXY(15);
199           fESDtrackCuts->SetMaxDCAToVertexZ(6);
200           fESDtrackCuts->SetDCAToVertex2D(kFALSE);
201           fESDtrackCuts->SetRequireSigmaToVertex(kFALSE);
202
203           fESDtrackCuts->SetEtaRange(-0.9,0.9);
204   }
205   //
206   //
207   //
208   //
209   fESDTrackCutsMult = AliESDtrackCuts::GetStandardITSTPCTrackCuts2010(kTRUE);
210   fESDTrackCutsMult->SetEtaRange(-1.2,+1.2);
211   fESDTrackCutsMult->SetPtRange(0.15,1e10);
212
213 }
214
215
216 //________________________________________________________________________
217 void AliAnalysisCombinedHadronSpectra::UserCreateOutputObjects() 
218 {
219   // Create histograms
220   // Called once
221   fListHist = new TList();
222   fListHist->SetOwner(kTRUE);
223   //
224   const Int_t kPtBins = 35;
225   const Int_t kMultBins = 11;
226
227   Int_t kDcaBinsTemp = 76;
228   if (fEvenDCAbinning) kDcaBinsTemp = 150;
229   const Int_t kDcaBins = (const Int_t) kDcaBinsTemp;
230
231   const Float_t kDcaBinsTPConlyFactor = 5; //need to change binning of DCA plot for tpconly
232   // sort pT-bins ..
233   Double_t binsPt[77] = {0., 0.05, 0.1, 0.15, 0.2, 0.25, 0.3, 0.35, 0.4, 0.45, 0.5, 0.55, 0.6, 0.65, 0.7, 0.75, 0.8, 0.85, 0.9, 0.95, 1.0, 1.1, 1.2, 1.3, 1.4, 1.5, 1.6, 1.7, 1.8, 1.9, 2.0, 2.2, 2.4, 2.6, 2.8, 3.0};
234
235   Double_t binsDca[77] =  {-3,-2.85,-2.7,-2.55,-2.4,-2.25,-2.1,-1.95,-1.8,-1.65,-1.5,-1.35,-1.2,-1.05,-0.9,-0.75,-0.6,-0.45,-0.3,-0.285,-0.27,-0.255,-0.24,-0.225,-0.21,-0.195,-0.18,-0.165,-0.15,-0.135,-0.12,-0.105,-0.09,-0.075,-0.06,-0.045,-0.03,-0.015,0,0.015,0.03,0.045,0.06,0.075,0.09,0.105,0.12,0.135,0.15,0.165,0.18,0.195,0.21,0.225,0.24,0.255,0.27,0.285,0.3,0.45,0.6,0.75,0.9,1.05,1.2,1.35,1.5,1.65,1.8,1.95,2.1,2.25,2.4,2.55,2.7,2.85,3};
236
237   // DCA bins borders get multiplied by constant factor for TPConlyTracks
238   Double_t binsDcaTPConly[77];
239   for (Int_t i = 0; i< 77; i++) {
240         binsDcaTPConly[i] = kDcaBinsTPConlyFactor * binsDca[i];
241   }
242
243   //
244   // create the histograms with all necessary information --> it is filled 4x for each particle assumption
245   //
246   // (0.) assumed particle: 0. pion, 1. kaon, 2. proton, 3. deuteron
247   // (1.) multiplicity or centrality -- number of accepted ESD tracks per events (deprecated), but now classes from 1 to 10, 0: Min. Bias
248   // (2.) pT
249   // (3.) sign
250   // (4.) rapidity --> filled 4x
251   // (5.)  pull TPC dEx --> filled 4x
252   // (6.) has valid TOF pid signal
253   // (7.) nsigma TOF --> filled 4x
254   // (8..) dca_xy
255   // (9.) CODE -- only MC 0-generated, 1-true rec. primaries, 2-misident prim, 3-second weak, 4-second material, 5-misident sec, 6-sec. K0, 7-sec. lambda, 8-sec sigma+
256   //
257   
258   //dimensions of standard THnSparse
259   //                              0,           1,           2,  3,       4,    5,     6,   7,     8
260   Int_t    binsHistReal[9] = {   3,   kMultBins,     kPtBins,  2,       10,    50,    2,  80,    kDcaBins};
261   Double_t xminHistReal[9] = {-0.5,        -0.5,           0, -2,      -0.5,   -5, -0.5,  -8,       -3};
262   Double_t xmaxHistReal[9] = { 2.5,        10.5,           3,  2,       0.5,    5,  1.5,   8,        3};
263
264   //dimensions of small THnSparse
265   //                              0,           1,           2,  3,                        4,   5,       6
266   Int_t    binsHistRealSm[7] = {   3,   kMultBins,     kPtBins,  2,  /*    10,   50,*/    2,  80, kDcaBins};
267   Double_t xminHistRealSm[7] = {-0.5,        -0.5,           0, -2,  /*  -0.5,   -5,*/ -0.5,  -8,       -3};
268   Double_t xmaxHistRealSm[7] = { 2.5,        10.5,           3,  2,  /*   0.5,    5,*/  1.5,   8,        3};
269
270   if (!fSmallTHnSparse) fHistRealTracks = new THnSparseF("fHistRealTracks","real tracks",9,binsHistReal,xminHistReal,xmaxHistReal);
271   else                  fHistRealTracks = new THnSparseF("fHistRealTracks","real tracks",7,binsHistRealSm,xminHistRealSm,xmaxHistRealSm);
272   //
273   fHistRealTracks->GetAxis(2)->Set(kPtBins, binsPt);
274
275   //different DCAxy binning for TPConlyTracks
276
277   Int_t dcaAxisNumber = 8;
278   if (fSmallTHnSparse) dcaAxisNumber = 6;
279
280   if (!fUseTPConlyTracks){
281     if (fEvenDCAbinning) fHistRealTracks->GetAxis(dcaAxisNumber)->Set(kDcaBins,-3.,3.);
282     else fHistRealTracks->GetAxis(dcaAxisNumber)->Set(kDcaBins, binsDca);
283   }
284   else {
285     if (fEvenDCAbinning) fHistRealTracks->GetAxis(dcaAxisNumber)->Set(kDcaBins,-15.,15.);
286     else fHistRealTracks->GetAxis(dcaAxisNumber)->Set(kDcaBins, binsDcaTPConly);
287    }
288   fListHist->Add(fHistRealTracks);
289   //
290   //                      0.ptot,1.tpcSig,2.hasTOF, 3. assumed part., 4. nclDedx, 5. nSigmaTPC (4x), 6. nSigmaTOF (4x), 7. centrality
291   fHistPidQA = new TH3D("fHistPidQA","PID QA",500,0.1,10,1000,0,1000,2,-2,2);
292   BinLogAxis(fHistPidQA);
293   fListHist->Add(fHistPidQA);
294   
295   // dimensions of standard THnSparse
296   //                            0,            1,           2,  3,      4,   5,    6,   7,       8.,    9.
297   Int_t    binsHistMC[10] = {   3,    kMultBins,     kPtBins,  2,     10,  50,    2,  80, kDcaBins,    6};
298   Double_t xminHistMC[10] = {-0.5,         -0.5,           0, -2,   -0.5,  -5, -0.5,  -8,       -3, -0.5};
299   Double_t xmaxHistMC[10] = { 2.5,         10.5,           3,  2,    0.5,   5,  1.5,   8,        3,  5.5};
300
301   //dimensions of small THnSparse
302   //                              0,           1,           2,  3,                     4,   5,       6.,    7.
303   Int_t    binsHistMCSm[8] = {   3,    kMultBins,     kPtBins,  2,  /*   10,  50,*/    2,  80, kDcaBins,    6};
304   Double_t xminHistMCSm[8] = {-0.5,         -0.5,           0, -2,  /* -0.5,  -5,*/ -0.5,  -8,       -3, -0.5};
305   Double_t xmaxHistMCSm[8] = { 2.5,         10.5,           3,  2,  /*  0.5,   5,*/  1.5,   8,        3,  5.5};
306
307   //different binning for CODE axis, if we want to save motherPDG
308   if (fSaveMotherPDG) {
309     binsHistMC[9] = 9;
310     xmaxHistMC[9] = 8.5;
311     binsHistMCSm[7] = 9;
312     xmaxHistMCSm[7] = 8.5;
313   }
314
315
316   if (!fSmallTHnSparse) fHistMCparticles = new THnSparseF("fHistMCparticles","MC histogram",10,binsHistMC,xminHistMC,xmaxHistMC);
317   else                  fHistMCparticles = new THnSparseF("fHistMCparticles","MC histogram",8,binsHistMCSm,xminHistMCSm,xmaxHistMCSm);
318
319
320   fHistMCparticles->GetAxis(2)->Set(kPtBins, binsPt);
321
322   /*
323   //different DCAxy binning for TPConlyTracks
324   if (!fEvenDCAbinning){
325     if (!fUseTPConlyTracks) fHistMCparticles->GetAxis(dcaAxisNumber)->Set(kDcaBins, binsDca);
326     else fHistMCparticles->GetAxis(dcaAxisNumber)->Set(kDcaBins, binsDcaTPConly);
327   }
328   */
329
330   if (!fUseTPConlyTracks){
331     if (fEvenDCAbinning) fHistMCparticles->GetAxis(dcaAxisNumber)->Set(kDcaBins,-3.,3.);
332     else fHistMCparticles->GetAxis(dcaAxisNumber)->Set(kDcaBins, binsDca);
333   }
334   else {
335     if (fEvenDCAbinning) fHistMCparticles->GetAxis(dcaAxisNumber)->Set(kDcaBins,-15.,15.);
336     else fHistMCparticles->GetAxis(dcaAxisNumber)->Set(kDcaBins, binsDcaTPConly);
337    }
338
339
340
341   fListHist->Add(fHistMCparticles);
342   //
343   fHistMult = new TH2D("fHistMult", "control histogram to count number of events", 502, -2.5, 499.5,4,-0.5,3.5);
344   fHistCentrality = new TH1D("fHistCentrality", "control histogram to count number of events", 22, -1.5, 20.5);
345   fListHist->Add(fHistMult);
346   fListHist->Add(fHistCentrality);
347   
348 }
349
350 //________________________________________________________________________
351 void AliAnalysisCombinedHadronSpectra::UserExec(Option_t *) 
352 {
353   //
354   // main event loop
355   //
356   if (!fESDpid) fESDpid = ((AliESDInputHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->GetESDpid();
357   if (!fESDpid) {
358     fESDpid = new AliESDpid(); // HACK FOR MC PBPB --> PLEASE REMOVE AS SOON AS POSSIBLE
359     fESDpid->GetTPCResponse().SetBetheBlochParameters(1.28778e+00/50., 3.13539e+01, TMath::Exp(-3.16327e+01), 1.87901e+00, 6.41583e+00);
360   }
361   //AliLog::SetGlobalLogLevel(AliLog::kError);
362   //
363   // Check Monte Carlo information and other access first:
364   //
365   AliMCEventHandler* eventHandler = dynamic_cast<AliMCEventHandler*> (AliAnalysisManager::GetAnalysisManager()->GetMCtruthEventHandler());
366   if (!eventHandler) {
367     //Printf("ERROR: Could not retrieve MC event handler");
368     fMCtrue = kFALSE;
369   }
370   //
371   AliMCEvent* mcEvent = 0x0;
372   AliStack* stack = 0x0;
373   if (eventHandler) mcEvent = eventHandler->MCEvent();
374   if (!mcEvent) {
375     //Printf("ERROR: Could not retrieve MC event");
376     if (fMCtrue) return;
377   }
378   if (fMCtrue) {
379     stack = mcEvent->Stack();
380     if (!stack) return;
381   }
382   //
383   fESD = dynamic_cast<AliESDEvent*>( InputEvent() );
384   if (!fESD) {
385     //Printf("ERROR: fESD not available");
386     return;
387   }
388   
389   if (!fESDtrackCuts) {
390     Printf("ERROR: fESDtrackCuts not available");
391     return;
392   }
393   //
394   // check if event is selected by physics selection class
395   //
396   Bool_t isSelected = kTRUE; // for reasons of backward compatibility --> check is now in AddTask macro
397   /*
398   Bool_t isSelected = kFALSE;
399   isSelected = (((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->IsEventSelected() & AliVEvent::kMB) == AliVEvent::kMB;
400   */
401   //
402   // monitor vertex position
403   //
404   const AliESDVertex *vertex = fESD->GetPrimaryVertexTracks();
405   if(vertex->GetNContributors()<1) {
406     // SPD vertex
407     vertex = fESD->GetPrimaryVertexSPD();
408     if(vertex->GetNContributors()<1) vertex = 0x0;
409   }  
410   //
411   // small track loop to determine trigger Pt, multiplicity or centrality
412   //
413   Double_t triggerPt = 0;
414   for (Int_t i=0;i<fESD->GetNumberOfTracks();++i) {
415     AliESDtrack *track =fESD->GetTrack(i);
416     if (!fESDTrackCutsMult->AcceptTrack(track)) continue;
417     if (track->Pt() > triggerPt) triggerPt = track->Pt();
418   }
419   Int_t trackCounter = fESDTrackCutsMult->CountAcceptedTracks(fESD);
420   //
421   // 2nd multiplicity estimator SPD hits; replaces multiplicity from HBT
422   //
423   Float_t spdCorr = -1;
424   const AliMultiplicity *mult = fESD->GetMultiplicity();
425   Float_t nClusters[6]={0.0,0.0,0.0,0.0,0.0,0.0};
426   for(Int_t ilay=0; ilay<6; ilay++) nClusters[ilay] = (Float_t)mult->GetNumberOfITSClusters(ilay);
427   if (vertex) spdCorr = AliESDUtils::GetCorrSPD2(nClusters[1],vertex->GetZ());
428   //
429   Float_t centrality = -1;
430   //
431   // IMPORTANT CENTRALITY DEFINITION FOR pp
432   //
433   if (!fUseHBTmultiplicity) {
434     // test binning to check if multiplicity-dependence is due to diffractive events being s
435     if (trackCounter >= 0  && trackCounter <= 0)  centrality = 0;
436     if (trackCounter >= 1  && trackCounter <= 1)  centrality = 1;
437     if (trackCounter >= 2  && trackCounter <= 2)  centrality = 2;
438     if (trackCounter >= 3  && trackCounter <= 3)  centrality = 3;
439     if (trackCounter >= 4  && trackCounter <= 4)  centrality = 4;
440     if (trackCounter >= 5  && trackCounter <= 5)  centrality = 5;
441     // this was all the original bin 1 being [0..5]
442     if (trackCounter >= 6  && trackCounter <= 9)  centrality = 6;
443     if (trackCounter >= 10 && trackCounter <= 14) centrality = 7;
444     if (trackCounter >= 15 && trackCounter <= 22) centrality = 8;
445     if (trackCounter >= 23 && trackCounter <= 32) centrality = 9;
446     if (trackCounter >= 33 && trackCounter <= 42) centrality = 10;
447     /*
448     if (trackCounter >= 43 && trackCounter <= 52) centrality = 7
449     if (trackCounter >= 53 && trackCounter <= 62) centrality = 8;
450     if (trackCounter >= 63 && trackCounter <= 72) centrality = 9;
451     if (trackCounter >= 73 && trackCounter <= 82) centrality = 10;
452     */
453   } else {
454     if (spdCorr >= 0  && spdCorr <=  2)  centrality  = 0;
455     if (spdCorr >= 3  && spdCorr <=  5)  centrality  = 1;
456     if (spdCorr >= 6  && spdCorr <=  8)  centrality  = 2;
457     if (spdCorr >= 9  && spdCorr <= 11)  centrality  = 3;
458     if (spdCorr >= 12 && spdCorr <= 14)  centrality  = 4;
459     if (spdCorr >= 15 && spdCorr <= 16)  centrality  = 5;
460     // this was all the original bin 1 being [0..16]
461     if (spdCorr >= 17 && spdCorr <= 30)  centrality =  6;
462     if (spdCorr >= 31 && spdCorr <= 45)  centrality =  7;
463     if (spdCorr >= 46 && spdCorr <= 68)  centrality =  8;
464     if (spdCorr >= 69 && spdCorr <= 97)  centrality =  9;
465     if (spdCorr >= 98)                   centrality = 10;
466     /*
467     if (spdCorr >= 17 && spdCorr <= 30)  centrality = 2;
468     if (spdCorr >= 31 && spdCorr <= 45)  centrality = 3;
469     if (spdCorr >= 46 && spdCorr <= 68)  centrality = 4;
470     if (spdCorr >= 69 && spdCorr <= 97)  centrality = 5;
471     if (spdCorr >= 98)                   centrality = 6;
472     */
473   }
474   //
475   Int_t rootS = fESD->GetBeamEnergy() < 1000 ? 0 : 1;
476   if (fESD->GetEventSpecie() == 4) { // PbPb
477     rootS = 2;
478     AliCentrality *esdCentrality = fESD->GetCentrality();
479     centrality = esdCentrality->GetCentralityClass10("V0M") + 1; // centrality percentile determined with V0
480     if (TMath::Abs(centrality - 1) < 1e-5) {
481       centrality = esdCentrality->GetCentralityClass5("V0M");
482     }
483   }
484   Int_t nContributors = 0;
485   if (fESD->GetPrimaryVertexTPC()) nContributors = fESD->GetPrimaryVertexTPC()->GetNContributors();
486   //
487   Int_t processtype = 0;
488   Int_t processCode = 0;
489   //
490   // important change: fill generated only after vertex cut in case of heavy-ions
491   //
492   if (!vertex && fESD->GetEventSpecie() == 4) {
493     fHistMult->Fill(-1, processCode);
494     PostData(1, fListHist);
495     return;
496   } else {
497     if (TMath::Abs(vertex->GetZv()) > 10 && fESD->GetEventSpecie() == 4) {
498       fHistMult->Fill(-1, processCode);
499       PostData(1, fListHist);
500       return;
501     }
502   }
503   //
504   if (fMCtrue) {
505     //
506     //
507     //
508     AliHeader * header = mcEvent->Header();
509     processtype = GetPythiaEventProcessType(header);
510     // non diffractive
511     if (processtype !=92 && processtype !=93 && processtype != 94) processCode = 1;
512     // single diffractive
513     if ((processtype == 92 || processtype == 93)) processCode = 2;
514     // double diffractive
515     if (processtype == 94) processCode = 3;
516     //
517     for(Int_t i = 0; i < stack->GetNtrack(); i++) {
518       TParticle * trackMC = stack->Particle(i);
519       Int_t pdg = trackMC->GetPdgCode();
520       //
521       Double_t xv = trackMC->Vx();
522       Double_t yv = trackMC->Vy();
523       Double_t zv = trackMC->Vz();
524       Double_t dxy = 0;
525       dxy = TMath::Sqrt(xv*xv + yv*yv); // so stupid to avoid warnings
526       Double_t dz = 0;
527       dz = TMath::Abs(zv); // so stupid to avoid warnings
528       //
529       // vertex cut - selection of primaries
530       //
531       //if (dxy > 3 || dz > 10) continue; // fixed cut at 3cm in r
532       //
533       if (!stack->IsPhysicalPrimary(i)) continue;
534       //
535       // fill MC histograms here...
536       // 
537       Double_t rap = trackMC->Y();
538       Double_t pT  = trackMC->Pt();
539       Int_t sign = pdg < 0 ? -1 : 1; // only works for charged pi,K,p !!
540 //      Double_t transMass = TMath::Sqrt(trackMC->Pt()*trackMC->Pt() + trackMC->GetMass()*trackMC->GetMass()) - trackMC->GetMass();
541       //
542       Int_t iPart = -1;
543       if (TMath::Abs(pdg) == 211)  iPart = 0; // select Pi+/Pi- only
544       if (TMath::Abs(pdg) == 321)  iPart = 1; // select K+/K- only
545       if (TMath::Abs(pdg) == 2212) iPart = 2; // select p+/p- only
546       if (iPart == -1) continue;
547       //
548
549       if (!fSmallTHnSparse){
550         Double_t vecHistMC[10] = {static_cast<Double_t>(iPart), centrality,  pT, static_cast<Double_t>(sign), rap, 0, 1, 0, dxy, 0};
551         if (!fOnlyQA) fHistMCparticles->Fill(vecHistMC);
552       }
553       else{
554         if (rap>fRapidityCutLow && rap<fRapidityCutHigh){
555           Double_t vecHistMC[8] = {static_cast<Double_t>(iPart), centrality,  pT, static_cast<Double_t>(sign), 1, 0, dxy, 0};
556           if (!fOnlyQA) fHistMCparticles->Fill(vecHistMC);
557         }
558       }
559     }
560   }
561   //
562   if (!isSelected && !fOnlyQA) {
563     PostData(1, fListHist);
564     return;
565   }
566   //
567   if (!vertex) {
568     fHistMult->Fill(-1, processCode);
569     PostData(1, fListHist);
570     return;
571   } else {
572     if (TMath::Abs(vertex->GetZv()) > 10) {
573       fHistMult->Fill(-1, processCode);
574       PostData(1, fListHist);
575       return;
576     }
577   }
578   //
579   // count events after physics selection and after vertex selection
580   //
581   //cout << "MULTIPLICITY " << trackCounter << " " << fESD->GetEventNumberInFile() <<endl;
582   fHistMult->Fill(trackCounter, processCode);
583   fHistCentrality->Fill(centrality);
584
585
586   //***************************************************
587   // track loop
588   //***************************************************
589   //const Float_t kNsigmaCut = 3;
590   //const Float_t k2sigmaCorr = 1/(0.5*(TMath::Erf(kNsigmaCut/sqrt(2))-TMath::Erf(-kNsigmaCut/sqrt(2))))/*1/0.9545*/;
591   //
592   Float_t dca[2], cov[3]; // dca_xy, dca_z, sigma_xy, sigma_xy_z, sigma_z for the vertex cut
593   //
594   for (Int_t i=0;i<fESD->GetNumberOfTracks();++i) {
595         
596         AliESDtrack *track = 0;
597         AliESDtrack *trackForTOF = 0; //normal track for all TOF information needed when using tpconly-tracks
598         
599         //normal tracks, if tpconly flag is set, use tpconlytracks
600         if (!fUseTPConlyTracks){
601         track =fESD->GetTrack(i); 
602         }
603         else {
604                 track = fESDtrackCuts->GetTPCOnlyTrack(fESD,i);
605                 if (!track) continue;
606                 trackForTOF = fESD->GetTrack(i);
607         }
608     //
609     if (!track->GetInnerParam()) {
610                 if (fUseTPConlyTracks) {delete track; track = 0;} //need to delete tpconlytrack
611                 continue;
612         }
613     Double_t ptot = track->GetInnerParam()->GetP(); // momentum for dEdx determination
614     Double_t pT = track->Pt();
615     track->GetImpactParameters(dca, cov);
616     //
617     //
618     // cut for dead regions in the detector
619     // if (track->Eta() > 0.1 && (track->Eta() < 0.2 && track->Phi() > 0.1 && track->Phi() < 0.1) continue;
620     //
621     // 2.a) apply some standard track cuts according to general recommendations
622     //
623     if (!fESDtrackCuts->AcceptTrack(track)) {
624                 if (fUseTPConlyTracks) {delete track; track = 0;} //need to delete tpconlytrack
625                 continue;
626         }
627
628         UInt_t status = 0;
629     if (!fUseTPConlyTracks) status = track->GetStatus();
630         else status = trackForTOF->GetStatus();
631     Bool_t hasTOFout  = status&AliESDtrack::kTOFout; 
632     Bool_t hasTOFtime = status&AliESDtrack::kTIME;
633     Bool_t hasTOFpid  = status&AliESDtrack::kTOFpid;
634     Bool_t hasTOF     = kFALSE;
635     if (hasTOFout && hasTOFtime && hasTOFpid) hasTOF = kTRUE;
636         Float_t length = 0.;
637     if (!fUseTPConlyTracks) length = track->GetIntegratedLength(); 
638     else length = trackForTOF->GetIntegratedLength();
639
640     if (length < 350.) hasTOF = kFALSE;
641     //
642     // calculate rapidities and kinematics
643     // 
644     //
645     Double_t pvec[3];
646     track->GetPxPyPz(pvec);
647     Double_t energyPion = TMath::Sqrt(track->GetP()*track->GetP() + AliPID::ParticleMass(AliPID::kPion)*AliPID::ParticleMass(AliPID::kPion));
648     Double_t energyKaon = TMath::Sqrt(track->GetP()*track->GetP() + AliPID::ParticleMass(AliPID::kKaon)*AliPID::ParticleMass(AliPID::kKaon));
649     Double_t energyProton = TMath::Sqrt(track->GetP()*track->GetP() + AliPID::ParticleMass(AliPID::kProton)*AliPID::ParticleMass(AliPID::kProton));
650     Double_t energyDeuteron = TMath::Sqrt(track->GetP()*track->GetP() + 4*AliPID::ParticleMass(AliPID::kProton)*AliPID::ParticleMass(AliPID::kProton));
651     //
652     Double_t rapPion = 0.5*TMath::Log((energyPion + pvec[2])/(energyPion - pvec[2]));
653     Double_t rapKaon = 0.5*TMath::Log((energyKaon + pvec[2])/(energyKaon - pvec[2]));
654     Double_t rapProton = 0.5*TMath::Log((energyProton + pvec[2])/(energyProton - pvec[2]));
655     Double_t rapDeuteron = 0.5*TMath::Log((energyDeuteron + pvec[2])/(energyDeuteron - pvec[2]));
656     //
657 //    Double_t transMassPion = TMath::Sqrt(track->Pt()*track->Pt() + AliPID::ParticleMass(AliPID::kPion)*AliPID::ParticleMass(AliPID::kPion))      -  AliPID::ParticleMass(AliPID::kPion);
658 //    Double_t transMassKaon = TMath::Sqrt(track->Pt()*track->Pt() + AliPID::ParticleMass(AliPID::kKaon)*AliPID::ParticleMass(AliPID::kKaon))     -  AliPID::ParticleMass(AliPID::kKaon);
659  //   Double_t transMassProton = TMath::Sqrt(track->Pt()*track->Pt() + AliPID::ParticleMass(AliPID::kProton)*AliPID::ParticleMass(AliPID::kProton))    -  AliPID::ParticleMass(AliPID::kProton);
660 //    Double_t transMassDeuteron = TMath::Sqrt(track->Pt()*track->Pt() + 4*AliPID::ParticleMass(AliPID::kProton)*AliPID::ParticleMass(AliPID::kProton))    -  2*AliPID::ParticleMass(AliPID::kProton);
661     //
662     // 3. make the PID
663     //
664     Double_t sign = track->GetSign();   
665     Double_t tpcSignal = track->GetTPCsignal();
666     //
667     //
668     // 3.a. calculate expected signals in nsigma
669     //
670     //  
671     // (0.) assumed particle: 0. pion, 1. kaon, 2. proton, 3. deuteron
672     // (1.) multiplicity or centrality -- number of accepted ESD tracks per events (deprecated), but now classes from 1 to 10, 0: Min. Bias
673     // (2.) pT
674     // (3.) sign
675     // (4.) rapidity --> filled 4x
676     // (5.)  pull TPC dEx --> filled 4x
677     // (6.) has valid TOF pid signal
678     // (7.) nsigma TOF --> filled 4x
679     // (8..) dca_xy
680     // (9.) CODE -- only MC 0-generated, 1-true rec. primaries, 2-misident, 3-second weak, 4-second material
681     //
682 //    Double_t transMass[4] = {transMassPion,transMassKaon,transMassProton,transMassDeuteron};
683     Double_t rap[4] = {rapPion,rapKaon,rapProton,rapDeuteron};
684     Double_t pullsTPC[4] = {fESDpid->NumberOfSigmasTPC(track,AliPID::kPion),
685                             fESDpid->NumberOfSigmasTPC(track,AliPID::kKaon),
686                             fESDpid->NumberOfSigmasTPC(track,AliPID::kProton),
687                             0}; // ASK FOR PUTTING THE DEUTERON TO AliPID !!!!!!!!!!!!!!
688     Float_t time0 = fESDpid->GetTOFResponse().GetTimeZero();
689     //fESDpid->GetTOFResponse().SetTimeResolution(130.);
690     Double_t pullsTOF[4] ={0.,0.,0.,0.};
691     if (!fUseTPConlyTracks) {
692       pullsTOF[0] = fESDpid->NumberOfSigmasTOF(track,AliPID::kPion, time0);
693       pullsTOF[1] = fESDpid->NumberOfSigmasTOF(track,AliPID::kKaon, time0);
694       pullsTOF[2] = fESDpid->NumberOfSigmasTOF(track,AliPID::kProton, time0);
695       pullsTOF[3] = 0; // ASK FOR PUTTING THE DEUTERON TO AliPID !!!!!!!!!!!!!!;
696     }
697     else {
698       pullsTOF[0] = fESDpid->NumberOfSigmasTOF(trackForTOF,AliPID::kPion, time0);
699       pullsTOF[1] = fESDpid->NumberOfSigmasTOF(trackForTOF,AliPID::kKaon, time0);
700       pullsTOF[2] = fESDpid->NumberOfSigmasTOF(trackForTOF,AliPID::kProton, time0);
701       pullsTOF[3] = 0; // ASK FOR PUTTING THE DEUTERON TO AliPID !!!!!!!!!!!!!!;
702     }
703
704     //
705 //    Double_t tpcQA[4] = {fESDpid->NumberOfSigmasTPC(track,AliPID::kElectron),
706 //                       fESDpid->NumberOfSigmasTPC(track,AliPID::kPion),
707 //                       fESDpid->NumberOfSigmasTPC(track,AliPID::kKaon),
708 //                       fESDpid->NumberOfSigmasTPC(track,AliPID::kProton)};
709
710     Double_t tofQA[4] = {0.,0.,0.,0.}; 
711     if (!fUseTPConlyTracks) {
712       tofQA[0] = fESDpid->NumberOfSigmasTOF(track,AliPID::kElectron, time0);
713       tofQA[0] = fESDpid->NumberOfSigmasTOF(track,AliPID::kPion, time0);
714       tofQA[0] = fESDpid->NumberOfSigmasTOF(track,AliPID::kKaon, time0);
715       tofQA[0] = fESDpid->NumberOfSigmasTOF(track,AliPID::kProton, time0);
716     }
717     else{
718       tofQA[0] = fESDpid->NumberOfSigmasTOF(trackForTOF,AliPID::kElectron, time0);
719       tofQA[0] = fESDpid->NumberOfSigmasTOF(trackForTOF,AliPID::kPion, time0);
720       tofQA[0] = fESDpid->NumberOfSigmasTOF(trackForTOF,AliPID::kKaon, time0);
721       tofQA[0] = fESDpid->NumberOfSigmasTOF(trackForTOF,AliPID::kProton, time0);
722     }
723
724
725     //save information for every particle type  // loop over assumed particle type
726     for(Int_t iPart = 0; iPart < 3; iPart++) {
727
728       if (!fSmallTHnSparse) {
729         Double_t vecHistReal[9]  = {static_cast<Double_t>(iPart),  centrality,   pT, static_cast<Double_t>(sign),  rap[iPart], pullsTPC[iPart], static_cast<Double_t>(hasTOF), pullsTOF[iPart], dca[0]};
730         if (!fOnlyQA) fHistRealTracks->Fill(vecHistReal);
731       }
732       else {
733         if (pullsTPC[iPart]>fTPCnSigmaCutLow && pullsTPC[iPart]<fTPCnSigmaCutHigh && rap[iPart]>fRapidityCutLow && rap[iPart]<fRapidityCutHigh) {
734           Double_t vecHistReal[7]  = {static_cast<Double_t>(iPart),  centrality,   pT, static_cast<Double_t>(sign), static_cast<Double_t>(hasTOF), pullsTOF[iPart], dca[0]};
735           if (!fOnlyQA) fHistRealTracks->Fill(vecHistReal);
736         }
737       }
738
739
740       // using MC truth for precise efficiencies...
741       //
742       if (fMCtrue && !fOnlyQA) {
743         Int_t code = 9; // code: 0-generated, 1-true rec. primaries, 2-misident, 3-second weak, 4-second material
744         Int_t assumedPdg = 0;//2212(proton); 321(Kaon); 211(pion);
745         Int_t motherCode = -1;
746         if (iPart == 0) assumedPdg = 211;
747         if (iPart == 1) assumedPdg = 321;
748         if (iPart == 2) assumedPdg = 2212;
749         //
750         //
751         TParticle *trackMC = stack->Particle(TMath::Abs(track->GetLabel()));
752         Int_t pdg = TMath::Abs(trackMC->GetPdgCode());
753         //
754         if (pdg != assumedPdg && stack->IsPhysicalPrimary(TMath::Abs(track->GetLabel()))) code = 2;
755         if (pdg != assumedPdg && stack->IsSecondaryFromWeakDecay(TMath::Abs(track->GetLabel()))) code = 5;
756         if (pdg == assumedPdg && stack->IsPhysicalPrimary(TMath::Abs(track->GetLabel()))) code = 1;
757         if (pdg == assumedPdg && stack->IsSecondaryFromWeakDecay(TMath::Abs(track->GetLabel()))) {
758           code = 3;
759           if (fSaveMotherPDG){
760             TParticle *trackMother =  stack->Particle(TMath::Abs(trackMC->GetFirstMother()));
761             if (trackMother->GetPdgCode() == 310) motherCode = 6; //K0
762             if (trackMother->GetPdgCode() == 3122) motherCode = 7; //Lambda
763             if (trackMother->GetPdgCode() == 3222) motherCode = 8; //Sigma+
764           }
765         }
766         if (pdg == assumedPdg && stack->IsSecondaryFromMaterial(TMath::Abs(track->GetLabel()))) code = 4;
767         
768         //
769         // muons need special treatment, because they are indistinguishable from pions
770         //
771         if (iPart == 0 && pdg == 13  && stack->IsPhysicalPrimary(TMath::Abs(track->GetLabel()))) code = 1;
772         if (iPart == 0 && pdg == 13  && stack->IsSecondaryFromWeakDecay(TMath::Abs(track->GetLabel()))) code = 3;
773         //
774         // check TOF mismatch on MC basis with TOF label
775         //
776         Int_t tofLabel[3];
777         if (!fUseTPConlyTracks) track->GetTOFLabel(tofLabel);
778         else trackForTOF->GetTOFLabel(tofLabel);
779         if (TMath::Abs(track->GetLabel()) != TMath::Abs(tofLabel[0]) || tofLabel[1] > 0) hasTOF = kFALSE;
780         //
781         // IMPORTANT BIG PROBLEM HERE THE PROBABLILITY TO HAVE A PID SIGNAL MUST BE IN !!!!!!!!!!!!
782         //
783         if (!fSmallTHnSparse){
784           Double_t vectorHistMC[10] = {static_cast<Double_t>(iPart),  centrality,  pT, static_cast<Double_t>(sign),  rap[iPart], pullsTPC[iPart], static_cast<Double_t>(hasTOF), pullsTOF[iPart], dca[0], static_cast<Double_t>(code)};
785           if (!fOnlyQA) { 
786             fHistMCparticles->Fill(vectorHistMC);
787             if (motherCode != -1 && fSaveMotherPDG) { //if mother of weak decay is K0, lambda or sigma+ add track again with this information
788               Double_t vectorHistMCmother[10] = {static_cast<Double_t>(iPart),  centrality,  pT, static_cast<Double_t>(sign),  rap[iPart], pullsTPC[iPart], static_cast<Double_t>(hasTOF), pullsTOF[iPart], dca[0], static_cast<Double_t>(motherCode)};
789               fHistMCparticles->Fill(vectorHistMCmother);
790             }
791           }
792         }
793         else{
794           if (pullsTPC[iPart]>fTPCnSigmaCutLow && pullsTPC[iPart]<fTPCnSigmaCutHigh && rap[iPart]>fRapidityCutLow && rap[iPart]<fRapidityCutHigh) {
795             //                              0,           1,   2,    3,           4,               5,      6,               7,      8,   9
796             Double_t vectorHistMC[8] = {static_cast<Double_t>(iPart),  centrality,  pT, static_cast<Double_t>(sign), static_cast<Double_t>(hasTOF), pullsTOF[iPart], dca[0], static_cast<Double_t>(code)};
797             if (!fOnlyQA) { 
798               fHistMCparticles->Fill(vectorHistMC);
799               if (motherCode != -1 && fSaveMotherPDG) { //if mother of weak decay is K0, lambda or sigma+ add track again with this information
800                 Double_t vectorHistMCmother[8] = {static_cast<Double_t>(iPart),  centrality,  pT, static_cast<Double_t>(sign), static_cast<Double_t>(hasTOF), pullsTOF[iPart], dca[0], static_cast<Double_t>(motherCode)};
801                 fHistMCparticles->Fill(vectorHistMCmother);
802               }
803             }
804           }
805         }
806       }
807       //
808       //
809       Int_t tpcShared = track->GetTPCnclsS();
810       if (TMath::Abs(track->Eta()) < 0.8 && iPart == 0 && tpcShared < 4) fHistPidQA->Fill(ptot,tpcSignal,sign);
811     } // end loop over assumed particle type
812
813           //need to delete tpconlytrack
814           if (fUseTPConlyTracks){
815                 delete track; 
816                 track = 0;     
817           }
818
819   } // end of track loop
820   
821   // Post output data  
822   PostData(1, fListHist);
823   
824 }      
825
826
827 //________________________________________________________________________
828 void AliAnalysisCombinedHadronSpectra::Terminate(Option_t *) 
829 {
830   // Draw result to the screen
831   // Called once at the end of the query
832   Printf("*** CONSTRUCTOR CALLED ****");
833
834 }
835
836
837 //________________________________________________________________________
838 Bool_t AliAnalysisCombinedHadronSpectra::SelectOnImpPar(AliESDtrack* t) {
839   //
840   // cut on transverse impact parameter // DEPRECATED
841   //
842   Float_t d0z0[2],covd0z0[3];
843   t->GetImpactParameters(d0z0,covd0z0);
844   Float_t sigma= 0.0050+0.0060/TMath::Power(t->Pt(),0.9);
845   Float_t d0max = 7.*sigma;
846   //
847   Float_t sigmaZ = 0.0146+0.0070/TMath::Power(t->Pt(),1.114758);
848   if (t->Pt() > 1) sigmaZ = 0.0216;
849   Float_t d0maxZ = 5.*sigmaZ;
850   //
851   if(TMath::Abs(d0z0[0]) < d0max && TMath::Abs(d0z0[1]) < d0maxZ) return kTRUE;
852   return kFALSE;
853 }
854
855
856 //________________________________________________________________________
857 void AliAnalysisCombinedHadronSpectra::BinLogAxis(const TH1 *h) {
858   //
859   // Method for the correct logarithmic binning of histograms
860   //
861   TAxis *axis = h->GetXaxis();
862   int bins = axis->GetNbins();
863
864   Double_t from = axis->GetXmin();
865   Double_t to = axis->GetXmax();
866   Double_t *newBins = new Double_t[bins + 1];
867    
868   newBins[0] = from;
869   Double_t factor = pow(to/from, 1./bins);
870   
871   for (int i = 1; i <= bins; i++) {
872    newBins[i] = factor * newBins[i-1];
873   }
874   axis->Set(bins, newBins);
875   delete [] newBins;
876   
877 }
878
879
880 //________________________________________________________________________
881 Int_t AliAnalysisCombinedHadronSpectra::GetPythiaEventProcessType(const AliHeader* aHeader, Bool_t adebug) const {
882   //
883   // get the process type of the event.
884   //
885
886   // can only read pythia headers, either directly or from cocktalil header
887   AliGenPythiaEventHeader* pythiaGenHeader = dynamic_cast<AliGenPythiaEventHeader*>(aHeader->GenEventHeader());
888
889   if (!pythiaGenHeader) {
890
891     AliGenCocktailEventHeader* genCocktailHeader = dynamic_cast<AliGenCocktailEventHeader*>(aHeader->GenEventHeader());
892     if (!genCocktailHeader) {
893       //printf("AliAnalysisCombinedHadronSpectra::GetProcessType : Unknown header type (not Pythia or Cocktail). \n");
894       return -1;
895     }
896
897     TList* headerList = genCocktailHeader->GetHeaders();
898     if (!headerList) {
899       return -1;
900     }
901
902     for (Int_t i=0; i<headerList->GetEntries(); i++) {
903       pythiaGenHeader = dynamic_cast<AliGenPythiaEventHeader*>(headerList->At(i));
904       if (pythiaGenHeader)
905         break;
906     }
907
908     if (!pythiaGenHeader) {
909       //printf("AliAnalysisCombinedHadronSpectra::GetProcessType : Could not find Pythia header. \n");
910       return -1;
911     }
912   }
913
914   if (adebug) {
915     //printf("AliAnalysisCombinedHadronSpectra::GetProcessType : Pythia process type found: %d \n",pythiaGenHeader->ProcessType());
916   }
917
918   return pythiaGenHeader->ProcessType();
919 }