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