new functionality and data memebers added
[u/mrichter/AliRoot.git] / PWG0 / dNdPt / AlidNdPtHelper.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 // static dNdPt helper functions
17 //
18 // basic functionality to select events and tracks 
19 // for dNdPt analysis
20 //
21 // Origin: Jan Fiete Grosse-Oetringhaus
22 // Modified and Extended: Jacek Otwinowski 19/11/2009
23 //
24 // 
25
26 #include <TROOT.h>
27 #include <TCanvas.h>
28 #include <TF1.h>
29 #include <TH1.h>
30 #include <TH2.h>
31 #include <TH3.h>
32
33 #include <AliHeader.h>
34 #include <AliStack.h>
35 #include <AliLog.h>
36 #include <AliESD.h>
37 #include <AliESDEvent.h>
38 #include <AliMCEvent.h>
39 #include <AliESDVertex.h>
40 #include <AliVertexerTracks.h>
41 #include <AliMathBase.h>
42 #include <AliESDtrackCuts.h>
43 #include <AliTracker.h>
44 #include "dNdPt/AlidNdPtEventCuts.h"
45 #include "dNdPt/AlidNdPtAcceptanceCuts.h"
46 #include "dNdPt/AlidNdPtHelper.h"
47
48 //____________________________________________________________________
49 ClassImp(AlidNdPtHelper)
50
51 //____________________________________________________________________
52 const AliESDVertex* AlidNdPtHelper::GetVertex(AliESDEvent* const aEsd, AlidNdPtEventCuts *const evtCuts, AlidNdPtAcceptanceCuts *const accCuts, AliESDtrackCuts *const trackCuts, AnalysisMode analysisMode, Bool_t debug, Bool_t bRedoTPC, Bool_t bUseMeanVertex)
53 {
54   // Get the vertex from the ESD and returns it if the vertex is valid
55   //
56   // Second argument decides which vertex is used (this selects
57   // also the quality criteria that are applied)
58
59   if(!aEsd) 
60   { 
61     ::Error("AlidNdPtHelper::GetVertex()","esd event is NULL");
62     return NULL;  
63   }
64  
65   if(!evtCuts || !accCuts || !trackCuts) 
66   { 
67     ::Error("AlidNdPtHelper::GetVertex()","cuts not available");
68     return NULL;  
69   }
70
71   const AliESDVertex* vertex = 0;
72   AliESDVertex *initVertex = 0;
73   if (analysisMode == kSPD || analysisMode == kTPCITS || 
74       analysisMode == kTPCSPDvtx || analysisMode == kTPCSPDvtxUpdate || analysisMode == kTPCITSHybrid)
75   {
76     vertex = aEsd->GetPrimaryVertexSPD();
77     if (debug)
78       Printf("AlidNdPtHelper::GetVertex: Returning SPD vertex");
79   }
80   else if (analysisMode == kTPC)
81   {
82     if(bRedoTPC) {
83
84       Double_t kBz = aEsd->GetMagneticField();
85       AliVertexerTracks vertexer(kBz);
86
87       if(bUseMeanVertex) {
88          Double_t pos[3]={evtCuts->GetMeanXv(),evtCuts->GetMeanYv(),evtCuts->GetMeanZv()};
89          Double_t err[3]={evtCuts->GetSigmaMeanXv(),evtCuts->GetSigmaMeanYv(),evtCuts->GetSigmaMeanZv()};
90          initVertex = new AliESDVertex(pos,err);
91          vertexer.SetVtxStart(initVertex);
92          vertexer.SetConstraintOn();
93       }
94
95       Double_t maxDCAr = accCuts->GetMaxDCAr();
96       Double_t maxDCAz = accCuts->GetMaxDCAz();
97       Int_t minTPCClust = trackCuts->GetMinNClusterTPC();
98
99       //vertexer.SetTPCMode(Double_t dcacut=0.1, Double_t dcacutIter0=1.0, Double_t maxd0z0=5.0, Int_t minCls=10, Int_t mintrks=1, Double_t nsigma=3., Double_t mindetfitter=0.1, Double_t maxtgl=1.5, Double_t fidR=3., Double_t fidZ=30., Int_t finderAlgo=1, Int_t finderAlgoIter0=4);
100       vertexer.SetTPCMode(0.1,1.0,5.0,minTPCClust,1,3.,0.1,2.0,maxDCAr,maxDCAz,1,4);
101
102       // TPC track preselection
103       Int_t ntracks = aEsd->GetNumberOfTracks();
104       TObjArray array(ntracks);
105       UShort_t *id = new UShort_t[ntracks];
106
107       Int_t count=0;
108       for (Int_t i=0;i <ntracks; i++) {
109         AliESDtrack *t = aEsd->GetTrack(i);
110         if (!t) continue;
111         if (t->Charge() == 0) continue;
112         if (!t->GetTPCInnerParam()) continue;
113         if (t->GetTPCNcls()<vertexer.GetMinClusters()) continue;
114         AliExternalTrackParam  *tpcTrack  = new AliExternalTrackParam(*(t->GetTPCInnerParam()));
115         if(tpcTrack) { 
116           array.AddLast(tpcTrack);
117           id[count] = (UShort_t)t->GetID();
118           count++;
119         }
120       } 
121       AliESDVertex *vTPC = vertexer.VertexForSelectedTracks(&array,id, kTRUE, kTRUE, bUseMeanVertex);
122       
123       // set recreated TPC vertex
124       aEsd->SetPrimaryVertexTPC(vTPC);
125
126       for (Int_t i=0; i<aEsd->GetNumberOfTracks(); i++) {
127         AliESDtrack *t = aEsd->GetTrack(i);
128         if(!t) continue;
129
130         Double_t x[3]; t->GetXYZ(x);
131         Double_t b[3]; AliTracker::GetBxByBz(x,b);
132         t->RelateToVertexTPCBxByBz(vTPC, b, kVeryBig);
133       }
134       
135       delete vTPC;
136       array.Delete();
137       delete [] id; id=NULL;
138
139     }
140     vertex = aEsd->GetPrimaryVertexTPC();
141     if (debug)
142      Printf("AlidNdPtHelper::GetVertex: Returning vertex from tracks");
143    }
144    else
145      Printf("AlidNdPtHelper::GetVertex: ERROR: Invalid second argument %d", analysisMode);
146
147     if (!vertex) {
148      if (debug)
149       Printf("AlidNdPtHelper::GetVertex: No vertex found in ESD");
150       return 0;
151     }
152
153   if (debug)
154   {
155     Printf("AlidNdPtHelper::GetVertex: Returning valid vertex: %s", vertex->GetTitle());
156     vertex->Print();
157   }
158   
159   if(initVertex) delete initVertex; initVertex=NULL;
160   return vertex;
161 }
162
163 //____________________________________________________________________
164 Bool_t AlidNdPtHelper::TestRecVertex(const AliESDVertex* vertex, AnalysisMode analysisMode, Bool_t debug)
165 {
166   // Checks if a vertex meets the needed quality criteria
167   if(!vertex) return kFALSE;
168   if(!vertex->GetStatus()) return kFALSE;
169
170   Float_t requiredZResolution = -1;
171   if (analysisMode == kSPD || analysisMode == kTPCITS || 
172       analysisMode == kTPCSPDvtx || analysisMode == kTPCSPDvtxUpdate || analysisMode == kTPCITSHybrid)
173   {
174     requiredZResolution = 1000;
175   }
176   else if (analysisMode == kTPC)
177     requiredZResolution = 10.;
178
179   // check resolution
180   Double_t zRes = vertex->GetZRes();
181
182   if (zRes > requiredZResolution) {
183     if (debug)
184       Printf("AlidNdPtHelper::TestVertex: Resolution too poor %f (required: %f", zRes, requiredZResolution);
185     return kFALSE;
186   }
187  
188   if (vertex->IsFromVertexerZ())
189   {
190     if (vertex->GetDispersion() > 0.02) 
191     {
192       if (debug)
193         Printf("AliPWG0Helper::TestVertex: Delta Phi too large in Vertexer Z: %f (required: %f", vertex->GetDispersion(), 0.02);
194       return kFALSE;
195     }
196   }
197
198   /*
199   // check Ncontributors
200   if (vertex->GetNContributors() <= 0) {
201     if (debug){
202       Printf("AlidNdPtHelper::GetVertex: NContributors() <= 0: %d",vertex->GetNContributors());
203       Printf("AlidNdPtHelper::GetVertex: NIndices(): %d",vertex->GetNIndices());
204       vertex->Print();
205     }
206     return kFALSE;
207   }
208   */
209
210   return kTRUE;
211 }
212
213 //____________________________________________________________________
214 Bool_t AlidNdPtHelper::IsPrimaryParticle(AliStack* const stack, Int_t idx, ParticleMode particleMode)
215 {
216 //
217 // check primary particles 
218 // depending on the particle mode
219 //
220   if(!stack) return kFALSE;
221
222   TParticle* particle = stack->Particle(idx);
223   if (!particle) return  kFALSE;
224
225   // only charged particles
226   Double_t charge = particle->GetPDG()->Charge()/3.;
227   if (TMath::Abs(charge) < 0.001) return kFALSE;
228
229   Int_t pdg = TMath::Abs(particle->GetPdgCode());
230
231   // physical primary
232   Bool_t prim = stack->IsPhysicalPrimary(idx);
233
234   if(particleMode==kMCPion) {
235     if(prim && pdg==kPiPlus) return kTRUE;
236     else return kFALSE;
237   } 
238
239   if (particleMode==kMCKaon) {
240     if(prim && pdg==kKPlus) return kTRUE;
241     else return kFALSE;
242   }
243     
244   if (particleMode==kMCProton) {
245     if(prim && pdg==kProton) return kTRUE;
246     else return kFALSE;
247   }
248
249 return prim;
250 }
251
252 //____________________________________________________________________
253 Bool_t AlidNdPtHelper::IsCosmicTrack(AliESDtrack *const track1, AliESDtrack *const track2)
254 {
255 //
256 // check cosmic tracks
257 //
258   if(!track1) return kFALSE;
259   if(!track2) return kFALSE;
260
261   //
262   // cosmic tracks in TPC
263   //
264   //if( TMath::Abs( track1->Theta() - track2->Theta() ) < 0.004  && 
265   //  ((TMath::Abs(track1->Phi()-track2->Phi()) - TMath::Pi() )<0.004) && 
266   //if( track1->Pt() > 4.0 && (TMath::Abs(track1->Phi()-track2->Phi())-TMath::Pi())<0.1 
267   //    && (TMath::Abs(track1->Eta()+track2->Eta())-0.1) < 0.0 && (track1->Charge()+track2->Charge()) == 0)
268
269   //Float_t scaleF= 6.0;
270   if ( track1->Pt() > 4 && track2->Pt() > 4 && 
271        //(TMath::Abs(track1->GetSnp()-track2->GetSnp())-2.) < scaleF * TMath::Sqrt(track1->GetSigmaSnp2()+track2->GetSigmaSnp2()) &&
272        //TMath::Abs(track1->GetTgl()-track2->GetTgl())   < scaleF * TMath::Sqrt(track1->GetSigmaTgl2()+track2->GetSigmaTgl2()) &&
273        //TMath::Abs(track1->OneOverPt()-track2->OneOverPt()) < scaleF * TMath::Sqrt(track1->GetSigma1Pt2()+track2->GetSigma1Pt2()) && 
274        (track1->Charge()+track2->Charge()) == 0 && 
275        track1->Eta()*track2->Eta()<0.0 && TMath::Abs(track1->Eta()+track2->Eta())<0.03 &&
276        TMath::Abs(TMath::Abs(track1->Phi()-track2->Phi())-TMath::Pi())<0.1
277      )  
278   {
279     printf("COSMIC  candidate \n");
280
281     printf("track1->Pt() %f, track1->Theta() %f, track1->Eta() %f, track1->Phi() %f, track1->Charge() %d  \n", track1->Pt(), track1->Theta(), track1->Eta(), track1->Phi(), track1->Charge());
282     printf("track2->Pt() %f, track2->Theta() %f, track2->Eta() %f, track2->Phi() %f, track2->Charge() %d  \n", track2->Pt(), track2->Theta(), track2->Eta(), track2->Phi(), track2->Charge());
283     printf("dtheta %f, deta %f, dphi %f, dq %d  \n", track1->Theta()-track2->Theta(),  track1->Eta()-track2->Eta(), track1->Phi()-track2->Phi(), track1->Charge()+track2->Charge()); 
284     printf("dsphi %f, errsphi %f, dtanl %f, errtanl %f  \n", TMath::Abs(track1->GetSnp()-track2->GetSnp()), TMath::Sqrt(track1->GetSigmaSnp2()+track2->GetSigmaSnp2()), TMath::Abs(track1->GetTgl()-track2->GetTgl()), TMath::Sqrt(track1->GetSigmaTgl2()+track2->GetSigmaTgl2())); 
285     return kTRUE;
286   }
287      
288 return kFALSE; 
289 }
290
291 //____________________________________________________________________
292 void AlidNdPtHelper::PrintConf(AnalysisMode analysisMode, AliTriggerAnalysis::Trigger trigger)
293 {
294   //
295   // Prints the given configuration
296   //
297
298   TString str(">>>> Running with ");
299
300   switch (analysisMode)
301   {
302     case kInvalid: str += "invalid setting"; break;
303     case kSPD : str += "SPD-only"; break;
304     case kTPC : str += "TPC-only"; break;
305     case kTPCITS : str += "Global tracking"; break;
306     case kTPCSPDvtx : str += "TPC tracking + SPD event vertex"; break;
307     case kTPCSPDvtxUpdate : str += "TPC tracks updated with SPD event vertex point"; break;
308     case kTPCITSHybrid : str += "TPC tracking + ITS refit + >1 SPD cluster"; break;
309     case kMCRec : str += "TPC tracking + Replace rec. with MC values"; break;
310   }
311   str += " and trigger ";
312
313   str += AliTriggerAnalysis::GetTriggerName(trigger);
314
315   str += " <<<<";
316
317   Printf("%s", str.Data());
318 }
319
320 //____________________________________________________________________
321 Int_t AlidNdPtHelper::ConvertPdgToPid(TParticle *particle) {
322 //
323 // Convert Abs(pdg) to pid 
324 // (0 - e, 1 - muons, 2 - pions, 3 - kaons, 4 - protons, 5 -all rest)
325 //
326 Int_t pid=-1;
327
328   if (TMath::Abs(particle->GetPdgCode()) == kElectron)         { pid = 0; }
329   else if (TMath::Abs(particle->GetPdgCode()) == kMuonMinus) { pid = 1; }
330   else if (TMath::Abs(particle->GetPdgCode()) == kPiPlus)    { pid = 2; }
331   else if (TMath::Abs(particle->GetPdgCode()) == kKPlus)     { pid = 3; }
332   else if (TMath::Abs(particle->GetPdgCode()) == kProton)    { pid = 4; }
333   else                                                       { pid = 5; }
334
335 return pid;
336 }
337
338 //_____________________________________________________________________________
339 TH1F* AlidNdPtHelper::CreateResHisto(TH2F* const hRes2, TH1F **phMean, Int_t integ,  Bool_t drawBinFits, Int_t minHistEntries)
340 {
341 //
342 // Create mean and resolution 
343 // histograms
344 //
345   TVirtualPad* currentPad = gPad;
346   TAxis* axis = hRes2->GetXaxis();
347   Int_t nBins = axis->GetNbins();
348   Bool_t overflowBinFits = kFALSE;
349   TH1F* hRes, *hMean;
350   if (axis->GetXbins()->GetSize()){
351     hRes = new TH1F("hRes", "", nBins, axis->GetXbins()->GetArray());
352     hMean = new TH1F("hMean", "", nBins, axis->GetXbins()->GetArray());
353   }
354   else{
355     hRes = new TH1F("hRes", "", nBins, axis->GetXmin(), axis->GetXmax());
356     hMean = new TH1F("hMean", "", nBins, axis->GetXmin(), axis->GetXmax());
357
358   }
359   hRes->SetStats(false);
360   hRes->SetOption("E");
361   hRes->SetMinimum(0.);
362   //
363   hMean->SetStats(false);
364   hMean->SetOption("E");
365  
366   // create the fit function
367   TF1 * fitFunc = new TF1("G","[0]*exp(-(x-[1])*(x-[1])/(2.*[2]*[2]))",-3,3);
368   
369   fitFunc->SetLineWidth(2);
370   fitFunc->SetFillStyle(0);
371   // create canvas for fits
372   TCanvas* canBinFits = NULL;
373   Int_t nPads = (overflowBinFits) ? nBins+2 : nBins;
374   Int_t nx = Int_t(sqrt(nPads-1.));// + 1;
375   Int_t ny = (nPads-1) / nx + 1;
376   if (drawBinFits) {
377     canBinFits = (TCanvas*)gROOT->FindObject("canBinFits");
378     if (canBinFits) delete canBinFits;
379     canBinFits = new TCanvas("canBinFits", "fits of bins", 200, 100, 500, 700);
380     canBinFits->Divide(nx, ny);
381   }
382
383   // loop over x bins and fit projection
384   Int_t dBin = ((overflowBinFits) ? 1 : 0);
385   for (Int_t bin = 1-dBin; bin <= nBins+dBin; bin++) {
386     if (drawBinFits) canBinFits->cd(bin + dBin);
387     Int_t bin0=TMath::Max(bin-integ,0);
388     Int_t bin1=TMath::Min(bin+integ,nBins);
389     TH1D* hBin = hRes2->ProjectionY("hBin", bin0, bin1);
390     //    
391     if (hBin->GetEntries() > minHistEntries) {
392       fitFunc->SetParameters(hBin->GetMaximum(),hBin->GetMean(),hBin->GetRMS());
393       hBin->Fit(fitFunc,"s");
394       Double_t sigma = TMath::Abs(fitFunc->GetParameter(2));
395
396       if (sigma > 0.){
397         hRes->SetBinContent(bin, TMath::Abs(fitFunc->GetParameter(2)));
398         hMean->SetBinContent(bin, fitFunc->GetParameter(1));    
399       }
400       else{
401         hRes->SetBinContent(bin, 0.);
402         hMean->SetBinContent(bin,0);
403       }
404       hRes->SetBinError(bin, fitFunc->GetParError(2));
405       hMean->SetBinError(bin, fitFunc->GetParError(1));
406       
407       //
408       //
409
410     } else {
411       hRes->SetBinContent(bin, 0.);
412       hRes->SetBinError(bin, 0.);
413       hMean->SetBinContent(bin, 0.);
414       hMean->SetBinError(bin, 0.);
415     }
416     
417
418     if (drawBinFits) {
419       char name[256];
420       if (bin == 0) {
421         sprintf(name, "%s < %.4g", axis->GetTitle(), axis->GetBinUpEdge(bin));
422       } else if (bin == nBins+1) {
423         sprintf(name, "%.4g < %s", axis->GetBinLowEdge(bin), axis->GetTitle());
424       } else {
425         sprintf(name, "%.4g < %s < %.4g", axis->GetBinLowEdge(bin),
426                 axis->GetTitle(), axis->GetBinUpEdge(bin));
427       }
428       canBinFits->cd(bin + dBin);
429       hBin->SetTitle(name);
430       hBin->SetStats(kTRUE);
431       hBin->DrawCopy("E");
432       canBinFits->Update();
433       canBinFits->Modified();
434       canBinFits->Update();
435     }
436     
437     delete hBin;
438   }
439
440   delete fitFunc;
441   currentPad->cd();
442   *phMean = hMean;
443   return hRes;
444 }
445
446 //_____________________________________________________________________________
447 TH1F* AlidNdPtHelper::MakeResol(TH2F * his, Int_t integ, Bool_t type, Bool_t drawBins, Int_t minHistEntries){
448 // Create resolution histograms
449   
450      TH1F *hisr=0, *hism=0;
451      if (!gPad) new TCanvas;
452          hisr = CreateResHisto(his,&hism,integ,drawBins,minHistEntries);
453          if (type) return hism;
454          else return hisr;
455
456 return hisr;     
457 }
458
459 //_____________________________________________________________________________
460 TObjArray* AlidNdPtHelper::GetAllChargedTracks(AliESDEvent *esdEvent, AnalysisMode analysisMode)
461 {
462   //
463   // all charged TPC particles 
464   //
465   TObjArray *allTracks = new TObjArray();
466   if(!allTracks) return allTracks;
467
468   AliESDtrack *track=0;
469   for (Int_t iTrack = 0; iTrack < esdEvent->GetNumberOfTracks(); iTrack++) 
470   { 
471     if(analysisMode == AlidNdPtHelper::kTPC) { 
472       //
473       // track must be deleted by user 
474       // esd track parameters are replaced by TPCinner
475       //
476       track = AliESDtrackCuts::GetTPCOnlyTrack(esdEvent,iTrack);
477       if(!track) continue;
478     } 
479     else if (analysisMode == AlidNdPtHelper::kTPCSPDvtx || analysisMode == AlidNdPtHelper::kTPCSPDvtxUpdate)
480     {
481       //
482       // track must be deleted by the user 
483       // esd track parameters are replaced by TPCinner
484       //
485       track = AlidNdPtHelper::GetTPCOnlyTrackSPDvtx(esdEvent,iTrack,kFALSE);
486       if(!track) continue;
487     }
488     else if( analysisMode == AlidNdPtHelper::kTPCITSHybrid )
489     {
490       track = AlidNdPtHelper::GetTrackSPDvtx(esdEvent,iTrack,kFALSE);
491     }
492     else 
493     {
494       track = esdEvent->GetTrack(iTrack);
495     }
496
497     if(!track) continue;
498
499     if(track->Charge()==0) { 
500       if(analysisMode == AlidNdPtHelper::kTPC || analysisMode == AlidNdPtHelper::kTPCSPDvtx || 
501          analysisMode == AlidNdPtHelper::kTPCSPDvtxUpdate) 
502       {
503         delete track; continue; 
504       } else {
505         continue;
506       } 
507     }
508
509     allTracks->Add(track);
510   }
511
512   if(analysisMode == AlidNdPtHelper::kTPC || analysisMode == AlidNdPtHelper::kTPCSPDvtx || 
513      analysisMode == AlidNdPtHelper::kTPCSPDvtxUpdate) {
514      
515      allTracks->SetOwner(kTRUE);
516   }
517
518 return allTracks;
519 }
520
521 //_____________________________________________________________________________
522 AliESDtrack *AlidNdPtHelper::GetTPCOnlyTrackSPDvtx(AliESDEvent* esdEvent, Int_t iTrack, Bool_t bUpdate)
523 {
524 //
525 // Create ESD tracks from TPCinner parameters.
526 // Propagte to DCA to SPD vertex.
527 // Update using SPD vertex point (parameter)
528 //
529 // It is user responsibility to delete these tracks
530 //
531
532   if (!esdEvent) return NULL;
533   if (!esdEvent->GetPrimaryVertexSPD() ) { return NULL; }
534   if (!esdEvent->GetPrimaryVertexSPD()->GetStatus() ) { return  NULL; }
535    
536   // 
537   AliESDtrack* track = esdEvent->GetTrack(iTrack);
538   if (!track)
539     return NULL;
540
541   Bool_t isOK = kFALSE;
542   Double_t x[3]; track->GetXYZ(x);
543   Double_t b[3]; AliTracker::GetBxByBz(x,b);
544
545   // create new ESD track
546   AliESDtrack *tpcTrack = new AliESDtrack();
547  
548   // relate TPC-only tracks (TPCinner) to SPD vertex
549   AliExternalTrackParam cParam;
550   if(bUpdate) {  
551     isOK = track->RelateToVertexTPCBxByBz(esdEvent->GetPrimaryVertexSPD(),b,kVeryBig,&cParam);
552     track->Set(cParam.GetX(),cParam.GetAlpha(),cParam.GetParameter(),cParam.GetCovariance());
553
554     // reject fake tracks
555     if(track->Pt() > 10000.)  {
556       ::Error("Exclude no physical tracks","pt>10000. GeV");
557       delete tpcTrack; 
558       return NULL;
559     }
560   }
561   else {
562     isOK = track->RelateToVertexTPCBxByBz(esdEvent->GetPrimaryVertexSPD(), b, kVeryBig);
563   }
564
565   // only true if we have a tpc track
566   if (!track->FillTPCOnlyTrack(*tpcTrack))
567   {
568     delete tpcTrack;
569     return NULL;
570   }
571   
572   if(!isOK) return NULL;
573
574 return tpcTrack;
575
576
577 //_____________________________________________________________________________
578 AliESDtrack *AlidNdPtHelper::GetTrackSPDvtx(AliESDEvent* esdEvent, Int_t iTrack, Bool_t bUpdate)
579 {
580 //
581 // Propagte track to DCA to SPD vertex.
582 // Update using SPD vertex point (parameter)
583 //
584   if (!esdEvent) return NULL;
585   if (!esdEvent->GetPrimaryVertexSPD() ) { return NULL; }
586   if (!esdEvent->GetPrimaryVertexSPD()->GetStatus() ) { return  NULL; }
587    
588   // 
589   AliESDtrack* track = esdEvent->GetTrack(iTrack);
590   if (!track)
591     return NULL;
592
593   Bool_t isOK = kFALSE;
594   Double_t x[3]; track->GetXYZ(x);
595   Double_t b[3]; AliTracker::GetBxByBz(x,b);
596
597   // relate tracks to SPD vertex
598   AliExternalTrackParam cParam;
599   if(bUpdate) {  
600     isOK = track->RelateToVertexBxByBz(esdEvent->GetPrimaryVertexSPD(),b,kVeryBig,&cParam);
601     track->Set(cParam.GetX(),cParam.GetAlpha(),cParam.GetParameter(),cParam.GetCovariance());
602
603     // reject fake tracks
604     if(track->Pt() > 10000.)  {
605       ::Error("Exclude no physical tracks","pt>10000. GeV");
606       return NULL;
607     }
608   }
609   else {
610     isOK = track->RelateToVertexBxByBz(esdEvent->GetPrimaryVertexSPD(), b, kVeryBig);
611   }
612  
613   if(!isOK) return NULL;
614
615 return track;
616
617
618 //_____________________________________________________________________________
619 Int_t AlidNdPtHelper::GetTPCMBTrackMult(AliESDEvent *const esdEvent, AlidNdPtEventCuts *const evtCuts, AlidNdPtAcceptanceCuts *const accCuts, AliESDtrackCuts *const trackCuts)
620 {
621   //
622   // get MB event track multiplicity
623   //
624   if(!esdEvent) 
625   { 
626     ::Error("AlidNdPtHelper::GetTPCMBTrackMult()","esd event is NULL");
627     return 0;  
628   }
629  
630   if(!evtCuts || !accCuts || !trackCuts) 
631   { 
632     ::Error("AlidNdPtHelper::GetTPCMBTrackMult()","cuts not available");
633     return 0;  
634   }
635
636   //
637   Double_t pos[3]={evtCuts->GetMeanXv(),evtCuts->GetMeanYv(),evtCuts->GetMeanZv()};
638   Double_t err[3]={evtCuts->GetSigmaMeanXv(),evtCuts->GetSigmaMeanYv(),evtCuts->GetSigmaMeanZv()};
639   AliESDVertex vtx0(pos,err);
640
641   //
642   Float_t maxDCAr = accCuts->GetMaxDCAr();
643   Float_t maxDCAz = accCuts->GetMaxDCAz();
644   Float_t minTPCClust = trackCuts->GetMinNClusterTPC();
645   //
646   Int_t ntracks = esdEvent->GetNumberOfTracks();
647   Double_t dca[2],cov[3];
648   Int_t mult=0;
649   for (Int_t i=0;i <ntracks; i++){
650     AliESDtrack *t = esdEvent->GetTrack(i);
651     if (!t) continue;
652     if (t->Charge() == 0) continue;
653     if (!t->GetTPCInnerParam()) continue;
654     if (t->GetTPCNcls()<minTPCClust) continue;
655     //
656     Double_t x[3]; t->GetXYZ(x);
657     Double_t b[3]; AliTracker::GetBxByBz(x,b);
658     const Double_t kMaxStep = 1;   //max step over the material
659
660     AliExternalTrackParam  *tpcTrack  = new AliExternalTrackParam(*(t->GetTPCInnerParam()));
661     if (!tpcTrack->PropagateToDCABxByBz(&vtx0,b,kMaxStep,dca,cov)) 
662     {
663       if(tpcTrack) delete tpcTrack; 
664       continue;
665     }
666     //
667     if (TMath::Abs(dca[0])>maxDCAr || TMath::Abs(dca[1])>maxDCAz) {
668       if(tpcTrack) delete tpcTrack; 
669       continue;
670     }
671
672     mult++;    
673
674     if(tpcTrack) delete tpcTrack; 
675   }
676
677 return mult;  
678 }
679
680 //_____________________________________________________________________________
681 Int_t AlidNdPtHelper::GetTPCMBPrimTrackMult(AliESDEvent *const esdEvent, AliStack *const  stack, AlidNdPtEventCuts *const evtCuts, AlidNdPtAcceptanceCuts *const accCuts, AliESDtrackCuts *const trackCuts)
682 {
683   //
684   // get MB primary event track multiplicity
685   //
686   if(!esdEvent) 
687   { 
688     ::Error("AlidNdPtHelper::GetTPCMBPrimTrackMult()","esd event is NULL");
689     return 0;  
690   }
691
692   if(!stack) 
693   { 
694     ::Error("AlidNdPtHelper::GetTPCMBPrimTrackMult()","esd event is NULL");
695     return 0;  
696   }
697  
698   if(!evtCuts || !accCuts || !trackCuts) 
699   { 
700     ::Error("AlidNdPtHelper::GetTPCMBPrimTrackMult()","cuts not available");
701     return 0;  
702   }
703
704   //
705   Double_t pos[3]={evtCuts->GetMeanXv(),evtCuts->GetMeanYv(),evtCuts->GetMeanZv()};
706   Double_t err[3]={evtCuts->GetSigmaMeanXv(),evtCuts->GetSigmaMeanYv(),evtCuts->GetSigmaMeanZv()};
707   AliESDVertex vtx0(pos,err);
708
709   //
710   Float_t maxDCAr = accCuts->GetMaxDCAr();
711   Float_t maxDCAz = accCuts->GetMaxDCAz();
712   Float_t minTPCClust = trackCuts->GetMinNClusterTPC();
713
714   //
715   Int_t ntracks = esdEvent->GetNumberOfTracks();
716   Double_t dca[2],cov[3];
717   Int_t mult=0;
718   for (Int_t i=0;i <ntracks; i++){
719     AliESDtrack *t = esdEvent->GetTrack(i);
720     if (!t) continue;
721     if (t->Charge() == 0) continue;
722     if (!t->GetTPCInnerParam()) continue;
723     if (t->GetTPCNcls()<minTPCClust) continue;
724     //
725     Double_t x[3]; t->GetXYZ(x);
726     Double_t b[3]; AliTracker::GetBxByBz(x,b);
727     const Double_t kMaxStep = 1;   //max step over the material
728
729     AliExternalTrackParam  *tpcTrack  = new AliExternalTrackParam(*(t->GetTPCInnerParam()));
730     if (!tpcTrack->PropagateToDCABxByBz(&vtx0,b,kMaxStep,dca,cov)) 
731     {
732       if(tpcTrack) delete tpcTrack; 
733       continue;
734     }
735     //
736     if (TMath::Abs(dca[0])>maxDCAr || TMath::Abs(dca[1])>maxDCAz) {
737       if(tpcTrack) delete tpcTrack; 
738       continue;
739     }
740
741     Int_t label = TMath::Abs(t->GetLabel());
742     TParticle *part = stack->Particle(label);
743     if(!part) { 
744       if(tpcTrack) delete tpcTrack; 
745       continue;
746     }
747     if(!stack->IsPhysicalPrimary(label)) 
748     { 
749       if(tpcTrack) delete tpcTrack; 
750       continue;
751     }
752
753     mult++;    
754
755     if(tpcTrack) delete tpcTrack; 
756   }
757
758 return mult;  
759 }
760
761
762
763
764
765 //_____________________________________________________________________________
766 Int_t AlidNdPtHelper::GetMCTrueTrackMult(AliMCEvent *const mcEvent, AlidNdPtEventCuts *const evtCuts, AlidNdPtAcceptanceCuts *const accCuts)
767 {
768   //
769   // calculate mc event true track multiplicity
770   //
771   if(!mcEvent) return 0;
772
773   AliStack* stack = 0;
774   Int_t mult = 0;
775
776   // MC particle stack
777   stack = mcEvent->Stack();
778   if (!stack) return 0;
779
780   Bool_t isEventOK = evtCuts->AcceptMCEvent(mcEvent);
781   if(!isEventOK) return 0; 
782
783   Int_t nPart  = stack->GetNtrack();
784   for (Int_t iMc = 0; iMc < nPart; ++iMc) 
785   {
786      TParticle* particle = stack->Particle(iMc);
787      if (!particle)
788      continue;
789
790      // only charged particles
791      Double_t charge = particle->GetPDG()->Charge()/3.;
792      if (TMath::Abs(charge) < 0.001)
793      continue;
794       
795      // physical primary
796      Bool_t prim = stack->IsPhysicalPrimary(iMc);
797      if(!prim) continue;
798
799      // checked accepted without pt cut
800      //if(accCuts->AcceptTrack(particle)) 
801      if( particle->Eta() > accCuts->GetMinEta() && particle->Eta() < accCuts->GetMaxEta() ) 
802      {
803        mult++;
804      }
805   }
806
807 return mult;  
808 }
809
810 //_______________________________________________________________________
811 void  AlidNdPtHelper::PrintMCInfo(AliStack *const pStack,Int_t label)
812 {
813 // print information about particles in the stack
814
815   if(!pStack)return;
816   label = TMath::Abs(label);
817   TParticle *part = pStack->Particle(label);
818   Printf("########################");
819   Printf("%s:%d %d UniqueID %d PDG %d P %3.3f",(char*)__FILE__,__LINE__,label,part->GetUniqueID(),part->GetPdgCode(),part->P());
820   part->Print();
821   TParticle* mother = part;
822   Int_t imo = part->GetFirstMother();
823   Int_t nprim = pStack->GetNprimary();
824
825   while((imo >= nprim)) {
826       mother =  pStack->Particle(imo);
827       Printf("Mother %s:%d Label %d UniqueID %d PDG %d P %3.3f",(char*)__FILE__,__LINE__,imo,mother->GetUniqueID(),mother->GetPdgCode(),mother->P());
828       mother->Print();
829       imo =  mother->GetFirstMother();
830  }
831
832  Printf("########################");
833 }
834
835
836 //_____________________________________________________________________________
837 TH1* AlidNdPtHelper::GetContCorrHisto(TH1 *const hist) 
838 {
839 //
840 // get contamination histogram
841 //
842  if(!hist) return 0;
843
844  Int_t nbins = hist->GetNbinsX();
845  TH1 *hCont = (TH1D *)hist->Clone();
846
847  for(Int_t i=0; i<=nbins+1; i++) {
848    Double_t binContent = hist->GetBinContent(i);
849    Double_t binError = hist->GetBinError(i);
850
851    hCont->SetBinContent(i,1.-binContent);
852    hCont->SetBinError(i,binError);
853  }
854
855 return hCont;
856 }
857
858
859 //_____________________________________________________________________________
860 TH1* AlidNdPtHelper::ScaleByBinWidth(TH1 *const hist) 
861 {
862 //
863 // scale by bin width
864 //
865  if(!hist) return 0;
866
867  TH1 *hScale = (TH1D *)hist->Clone();
868  hScale->Scale(1.,"width");
869
870 return hScale;
871 }
872
873 //_____________________________________________________________________________
874 TH1* AlidNdPtHelper::CalcRelativeDifference(TH1 *const hist1, TH1 *const hist2) 
875 {
876 //
877 // calculate rel. difference 
878 //
879
880  if(!hist1) return 0;
881  if(!hist2) return 0;
882
883  TH1 *h1Clone = (TH1D *)hist1->Clone();
884  h1Clone->Sumw2();
885
886  // (rec-mc)/mc
887  h1Clone->Add(hist2,-1);
888  h1Clone->Divide(hist2);
889
890 return h1Clone;
891 }
892
893 //_____________________________________________________________________________
894 TH1* AlidNdPtHelper::CalcRelativeDifferenceFun(TH1 *const hist1, TF1 *const fun) 
895 {
896 //
897 // calculate rel. difference
898 // between histogram and function
899 //
900  if(!hist1) return 0;
901  if(!fun) return 0;
902
903  TH1 *h1Clone = (TH1D *)hist1->Clone();
904  h1Clone->Sumw2();
905
906  // 
907  h1Clone->Add(fun,-1);
908  h1Clone->Divide(hist1);
909
910 return h1Clone;
911 }
912
913 //_____________________________________________________________________________
914 TH1* AlidNdPtHelper::NormalizeToEvent(TH2 *const hist1, TH1 *const hist2) 
915 {
916 // normalise to event for a given multiplicity bin
917 // return pt histogram 
918
919  if(!hist1) return 0;
920  if(!hist2) return 0;
921  char name[256];
922
923  Int_t nbinsX = hist1->GetNbinsX();
924  //Int_t nbinsY = hist1->GetNbinsY();
925
926  TH1D *histNorm = 0;
927  for(Int_t i=0; i<=nbinsX+1; i++) {
928    sprintf(name,"mom_%d",i);
929    TH1D *hist = (TH1D*)hist1->ProjectionY(name,i+1,i+1);
930
931    sprintf(name,"mom_norm");
932    if(i==0) { 
933      histNorm = (TH1D *)hist->Clone(name);
934      histNorm->Reset();
935    }
936
937    Double_t nbEvents = hist2->GetBinContent(i);
938    if(!nbEvents) { nbEvents = 1.; };
939
940    hist->Scale(1./nbEvents);
941    histNorm->Add(hist);
942  }
943
944 return histNorm;
945 }
946
947 //_____________________________________________________________________________
948 THnSparse* AlidNdPtHelper::GenerateCorrMatrix(THnSparse *const hist1, THnSparse *const hist2, char *const name) {
949 // generate correction matrix
950 if(!hist1 || !hist2) return 0; 
951
952 THnSparse *h =(THnSparse*)hist1->Clone(name);;
953 h->Divide(hist1,hist2,1,1,"B");
954
955 return h;
956 }
957
958 //_____________________________________________________________________________
959 TH2* AlidNdPtHelper::GenerateCorrMatrix(TH2 *const hist1, TH2 *const hist2, char *const name) {
960 // generate correction matrix
961 if(!hist1 || !hist2) return 0; 
962
963 TH2D *h =(TH2D*)hist1->Clone(name);;
964 h->Divide(hist1,hist2,1,1,"B");
965
966 return h;
967 }
968
969 //_____________________________________________________________________________
970 TH1* AlidNdPtHelper::GenerateCorrMatrix(TH1 *const hist1, TH1 *const hist2, char *const name) {
971 // generate correction matrix
972 if(!hist1 || !hist2) return 0; 
973
974 TH1D *h =(TH1D*)hist1->Clone(name);;
975 h->Divide(hist1,hist2,1,1,"B");
976
977 return h;
978 }
979
980 //_____________________________________________________________________________
981 THnSparse* AlidNdPtHelper::GenerateContCorrMatrix(THnSparse *const hist1, THnSparse *const hist2, char *const name) {
982 // generate contamination correction matrix
983 if(!hist1 || !hist2) return 0; 
984
985 THnSparse *hist =  GenerateCorrMatrix(hist1, hist2, name);
986 if(!hist) return 0;
987
988 // only for non ZERO bins!!!!
989
990 Int_t* coord = new Int_t[hist->GetNdimensions()];
991 memset(coord, 0, sizeof(Int_t) * hist->GetNdimensions());
992
993   for (Long64_t i = 0; i < hist->GetNbins(); ++i) {
994     Double_t v = hist->GetBinContent(i, coord);
995     hist->SetBinContent(coord, 1.0-v);
996     //printf("v %f, hist->GetBinContent(i, coord) %f \n",v,hist->GetBinContent(i, coord));
997     Double_t err = hist->GetBinError(coord);
998     hist->SetBinError(coord, err);
999   }
1000
1001 delete [] coord;
1002
1003 return hist;
1004 }
1005
1006 //_____________________________________________________________________________
1007 TH2* AlidNdPtHelper::GenerateContCorrMatrix(TH2 *const hist1, TH2 *const hist2, char *const name) {
1008 // generate contamination correction matrix
1009 if(!hist1 || !hist2) return 0; 
1010
1011 TH2 *hist = GenerateCorrMatrix(hist1, hist2, name);
1012 if(!hist) return 0;
1013
1014 Int_t nBinsX = hist->GetNbinsX();
1015 Int_t nBinsY = hist->GetNbinsY();
1016
1017   for (Int_t i = 0; i < nBinsX+1; i++) {
1018   for (Int_t j = 0; j < nBinsY+1; j++) {
1019      Double_t cont = hist->GetBinContent(i,j);
1020      hist->SetBinContent(i,j,1.-cont);
1021      Double_t err = hist->GetBinError(i,j);
1022      hist->SetBinError(i,j,err);
1023   }
1024   }
1025
1026 return hist;
1027 }
1028
1029 //_____________________________________________________________________________
1030 TH1* AlidNdPtHelper::GenerateContCorrMatrix(TH1 *const hist1, TH1 *const hist2, char *const name) {
1031 // generate contamination correction matrix
1032 if(!hist1 || !hist2) return 0; 
1033
1034 TH1 *hist = GenerateCorrMatrix(hist1, hist2, name);
1035 if(!hist) return 0;
1036
1037 Int_t nBinsX = hist->GetNbinsX();
1038
1039   for (Int_t i = 0; i < nBinsX+1; i++) {
1040      Double_t cont = hist->GetBinContent(i);
1041      hist->SetBinContent(i,1.-cont);
1042      Double_t err = hist->GetBinError(i);
1043      hist->SetBinError(i,err);
1044   }
1045
1046 return hist;
1047 }
1048
1049 //_____________________________________________________________________________
1050 const AliESDVertex* AlidNdPtHelper::GetTPCVertexZ(AliESDEvent* const esdEvent, AlidNdPtEventCuts *const evtCuts, AlidNdPtAcceptanceCuts *const accCuts, AliESDtrackCuts *const trackCuts, Float_t fraction, Int_t ntracksMin){
1051   //
1052   // TPC Z vertexer
1053   //
1054   if(!esdEvent)
1055   { 
1056     ::Error("AlidNdPtHelper::GetTPCVertexZ()","cuts not available");
1057     return NULL;  
1058   }
1059
1060   if(!evtCuts || !accCuts || !trackCuts) 
1061   { 
1062     ::Error("AlidNdPtHelper::GetTPCVertexZ()","cuts not available");
1063     return NULL;  
1064   }
1065
1066   Double_t vtxpos[3]={evtCuts->GetMeanXv(),evtCuts->GetMeanYv(),evtCuts->GetMeanZv()};
1067   Double_t vtxsigma[3]={evtCuts->GetSigmaMeanXv(),evtCuts->GetSigmaMeanYv(),evtCuts->GetSigmaMeanZv()};
1068   AliESDVertex vtx0(vtxpos,vtxsigma);
1069
1070   Double_t maxDCAr = accCuts->GetMaxDCAr();
1071   Double_t maxDCAz = accCuts->GetMaxDCAz();
1072   Int_t minTPCClust = trackCuts->GetMinNClusterTPC();
1073
1074   //
1075   Int_t ntracks = esdEvent->GetNumberOfTracks();
1076   TVectorD ztrack(ntracks);
1077   Double_t dca[2],cov[3];
1078   Int_t counter=0;
1079   for (Int_t i=0;i <ntracks; i++){
1080     AliESDtrack *t = esdEvent->GetTrack(i);
1081     if (!t) continue;
1082     if (!t->GetTPCInnerParam()) continue;
1083     if (t->GetTPCNcls()<minTPCClust) continue;
1084     //
1085
1086     Double_t x[3]; t->GetXYZ(x);
1087     Double_t b[3]; AliTracker::GetBxByBz(x,b);
1088     const Double_t kMaxStep = 1;   //max step over the material
1089
1090     AliExternalTrackParam  *tpcTrack  = new AliExternalTrackParam(*(t->GetTPCInnerParam()));
1091     if (!tpcTrack->PropagateToDCABxByBz(&vtx0,b,kMaxStep,dca,cov)) continue;
1092
1093     //
1094     if (TMath::Abs(dca[0])>maxDCAr) continue;
1095     //if (TMath::Sqrt(cov[0])>sigmaXYcut) continue;    
1096     if (TMath::Abs(tpcTrack->GetZ())>maxDCAz) continue;
1097
1098     ztrack[counter]=tpcTrack->GetZ();
1099     counter++;    
1100
1101     if(tpcTrack) delete tpcTrack;
1102   }
1103
1104   //
1105   // Find LTM z position
1106   //
1107   Double_t mean=0, sigma=0;
1108   if (counter<ntracksMin) return 0;
1109   //
1110   Int_t nused = TMath::Nint(counter*fraction);
1111   if (nused==counter) nused=counter-1;  
1112   if (nused>1){
1113     AliMathBase::EvaluateUni(counter, ztrack.GetMatrixArray(), mean,sigma, TMath::Nint(counter*fraction));
1114     sigma/=TMath::Sqrt(nused);
1115   }else{
1116     mean  = TMath::Mean(counter, ztrack.GetMatrixArray());
1117     sigma = TMath::RMS(counter, ztrack.GetMatrixArray());
1118     sigma/=TMath::Sqrt(counter-1);
1119   }
1120   vtxpos[2]=mean;
1121   vtxsigma[2]=sigma;
1122   const AliESDVertex* vertex = new AliESDVertex(vtxpos, vtxsigma);
1123   return vertex;
1124 }
1125
1126 //_____________________________________________________________________________
1127 Int_t  AlidNdPtHelper::GetSPDMBTrackMult(AliESDEvent* const esdEvent, Float_t deltaThetaCut, Float_t deltaPhiCut) 
1128 {
1129   //
1130   // SPD track multiplicity
1131   //
1132
1133   // get tracklets
1134   const AliMultiplicity* mult = esdEvent->GetMultiplicity();
1135   if (!mult)
1136      return 0;
1137
1138   // get multiplicity from SPD tracklets
1139   Int_t inputCount = 0; 
1140   for (Int_t i=0; i<mult->GetNumberOfTracklets(); ++i)
1141   {
1142     //printf("%d %f %f %f\n", i, mult->GetTheta(i), mult->GetPhi(i), mult->GetDeltaPhi(i));
1143
1144      Float_t phi = mult->GetPhi(i);
1145      if (phi < 0)
1146        phi += TMath::Pi() * 2;
1147      Float_t deltaPhi = mult->GetDeltaPhi(i);
1148      Float_t deltaTheta = mult->GetDeltaTheta(i);
1149
1150      if (TMath::Abs(deltaPhi) > 1)
1151        printf("WARNING: Very high Delta Phi: %d %f %f %f\n", i, mult->GetTheta(i), mult->GetPhi(i), deltaPhi);
1152
1153      if (deltaThetaCut > 0. && TMath::Abs(deltaTheta) > deltaThetaCut)
1154         continue;
1155
1156      if (deltaPhiCut > 0. && TMath::Abs(deltaPhi) > deltaPhiCut)
1157         continue;
1158       
1159      ++inputCount;
1160   }
1161
1162 return inputCount;
1163 }
1164
1165 //_____________________________________________________________________________
1166 Int_t  AlidNdPtHelper::GetSPDMBPrimTrackMult(AliESDEvent* const esdEvent, AliStack* const stack, Float_t deltaThetaCut, Float_t deltaPhiCut) 
1167 {
1168   //
1169   // SPD track multiplicity
1170   //
1171
1172   // get tracklets
1173   const AliMultiplicity* mult = esdEvent->GetMultiplicity();
1174   if (!mult)
1175      return 0;
1176
1177   // get multiplicity from SPD tracklets
1178   Int_t inputCount = 0; 
1179   for (Int_t i=0; i<mult->GetNumberOfTracklets(); ++i)
1180   {
1181     //printf("%d %f %f %f\n", i, mult->GetTheta(i), mult->GetPhi(i), mult->GetDeltaPhi(i));
1182
1183      Float_t phi = mult->GetPhi(i);
1184      if (phi < 0)
1185        phi += TMath::Pi() * 2;
1186      Float_t deltaPhi = mult->GetDeltaPhi(i);
1187      Float_t deltaTheta = mult->GetDeltaTheta(i);
1188
1189      if (TMath::Abs(deltaPhi) > 1)
1190        printf("WARNING: Very high Delta Phi: %d %f %f %f\n", i, mult->GetTheta(i), mult->GetPhi(i), deltaPhi);
1191
1192      if (deltaThetaCut > 0. && TMath::Abs(deltaTheta) > deltaThetaCut)
1193         continue;
1194
1195      if (deltaPhiCut > 0. && TMath::Abs(deltaPhi) > deltaPhiCut)
1196         continue;
1197
1198
1199      if (mult->GetLabel(i, 0) < 0 || mult->GetLabel(i, 0) != mult->GetLabel(i, 1) || 
1200          !stack->IsPhysicalPrimary(mult->GetLabel(i, 0)))
1201         continue;
1202
1203       
1204      ++inputCount;
1205   }
1206
1207 return inputCount;
1208 }
1209
1210