AliPhysics  fffcdf3 (fffcdf3)
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
AliAnalysisTaskFilteredTree.cxx
Go to the documentation of this file.
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 /*
17  Class to process/filter reconstruction information from ESD, ESD friends, MC and provide them for later reprocessing
18  Filtering schema - low pt part is downscaled - to have flat pt spectra of selected topologies (tracks and V0s)
19  Downscaling schema is controlled by downscaling factors
20  Usage:
21  1.) Filtering on Lego train
22  2.) expert QA for tracking (resolution efficiency)
23  3.) pt resolution studies using V0s
24  4.) dEdx calibration using V0s
25  5.) pt resolution and dEdx studies using cosmic
26  +
27  6.) Info used for later raw data OFFLINE triggering (highPt, V0, laser, cosmic, high dEdx)
28 
29  Exported trees (with full objects and derived variables):
30  1.) "highPt" - filtered trees with esd tracks, derived variables(propagated tracks), optional MC info +optional space points
31  2.) "V0s" - - filtered trees with selected V0s (rough KF chi2 cut), KF particle and corresponding esd tracks + optional space points
32  3.) "Laser" - dump laser tracks with space points if exists
33  4.) "CosmicTree" - cosmic track candidate (random or triggered) + esdTracks(up/down)+ optional points
34  5.) "dEdx" - tree with high dEdx tpc tracks
35 */
36 
37 #include "iostream"
38 #include "TSystem.h"
39 #include <TPDGCode.h>
40 #include <TDatabasePDG.h>
41 
42 #include "TChain.h"
43 #include "TTreeStream.h"
44 #include "TTree.h"
45 #include "TH1F.h"
46 #include "TH3.h"
47 #include "TCanvas.h"
48 #include "TList.h"
49 #include "TObjArray.h"
50 #include "TFile.h"
51 #include "TMatrixD.h"
52 #include "TRandom3.h"
53 
54 #include "AliHeader.h"
55 #include "AliGenEventHeader.h"
56 #include "AliInputEventHandler.h"
57 #include "AliStack.h"
58 #include "AliTrackReference.h"
59 #include "AliTrackPointArray.h"
60 #include "AliSysInfo.h"
61 
62 #include "AliPhysicsSelection.h"
63 #include "AliAnalysisTask.h"
64 #include "AliAnalysisManager.h"
65 #include "AliESDEvent.h"
66 #include "AliESDfriend.h"
67 #include "AliMCEvent.h"
68 #include "AliESDInputHandler.h"
69 #include "AliESDVertex.h"
70 #include "AliTracker.h"
71 #include "AliVTrack.h"
72 #include "AliGeomManager.h"
73 
74 #include "AliCentrality.h"
75 #include "AliESDVZERO.h"
76 #include "AliMultiplicity.h"
77 
78 #include "AliESDtrackCuts.h"
79 #include "AliMCEventHandler.h"
80 #include "AliFilteredTreeEventCuts.h"
81 #include "AliFilteredTreeAcceptanceCuts.h"
82 
83 #include "AliAnalysisTaskFilteredTree.h"
84 #include "AliKFParticle.h"
85 #include "AliESDv0.h"
86 #include "AliPID.h"
87 #include "AliPIDResponse.h"
88 #include "TVectorD.h"
89 #include "TStatToolkit.h"
90 using namespace std;
91 
92 ClassImp(AliAnalysisTaskFilteredTree)
93 
94  //_____________________________________________________________________________
95  AliAnalysisTaskFilteredTree::AliAnalysisTaskFilteredTree(const char *name)
96  : AliAnalysisTaskSE(name)
97  , fESD(0)
98  , fMC(0)
99  , fESDfriend(0)
100  , fOutput(0)
101  , fPitList(0)
102  , fUseMCInfo(kFALSE)
103  , fUseESDfriends(kFALSE)
104  , fReducePileUp(kTRUE)
105  , fFillTree(kTRUE)
106  , fFilteredTreeEventCuts(0)
107  , fFilteredTreeAcceptanceCuts(0)
108  , fFilteredTreeRecAcceptanceCuts(0)
109  , fEsdTrackCuts(0)
110  , fTrigger(AliTriggerAnalysis::kMB1)
111  , fAnalysisMode(kTPCAnalysisMode)
112  , fTreeSRedirector(0)
113  , fCentralityEstimator(0)
114  , fLowPtTrackDownscaligF(0)
115  , fLowPtV0DownscaligF(0)
116  , fFriendDownscaling(-3.)
117  , fProcessAll(kFALSE)
118  , fProcessCosmics(kFALSE)
119  , fProcessITSTPCmatchOut(kFALSE) // swittch to process ITS/TPC standalone tracks
120  , fHighPtTree(0)
121  , fV0Tree(0)
122  , fdEdxTree(0)
123  , fLaserTree(0)
124  , fMCEffTree(0)
125  , fCosmicPairsTree(0)
126  , fPtResPhiPtTPC(0)
127  , fPtResPhiPtTPCc(0)
128  , fPtResPhiPtTPCITS(0)
129  , fPtResEtaPtTPC(0)
130  , fPtResEtaPtTPCc(0)
131  , fPtResEtaPtTPCITS(0)
132  , fPtResCentPtTPC(0)
133  , fPtResCentPtTPCc(0)
134  , fPtResCentPtTPCITS(0)
135  , fCurrentFileName("")
136  , fDummyTrack(0)
137 {
138  // Constructor
139 
140  // Define input and output slots here
141  DefineOutput(1, TTree::Class());
142  DefineOutput(2, TTree::Class());
143  DefineOutput(3, TTree::Class());
144  DefineOutput(4, TTree::Class());
145  DefineOutput(5, TTree::Class());
146  DefineOutput(6, TTree::Class());
147  DefineOutput(7, TList::Class());
148 }
149 
150 //_____________________________________________________________________________
151 AliAnalysisTaskFilteredTree::~AliAnalysisTaskFilteredTree()
152 {
153  //
154  // Destructor
155  //
156  delete fFilteredTreeEventCuts;
157  delete fFilteredTreeAcceptanceCuts;
158  delete fFilteredTreeRecAcceptanceCuts;
159  delete fEsdTrackCuts;
160 }
161 
162 //____________________________________________________________________________
163 Bool_t AliAnalysisTaskFilteredTree::Notify()
164 {
165  //
166  //
167 
168  static Int_t count = 0;
169  count++;
170  TTree *chain = (TChain*)GetInputData(0);
171  if(chain){
172  Printf("Processing %d. file: %s", count, chain->GetCurrentFile()->GetName());
173  }
174  fCurrentFileName=chain->GetCurrentFile()->GetName(); // This is the way to get path infor processing on the Lego train
175  //
176  if (fCurrentFileName.String().CountChar('/')<1){
177  // in case path to input is not absolute use other handle to guess chunk ID
178  // in case filtering done locally, use the Env varaible AliEn output path
179  // this will work only in case we are doing alien processing
180  TString fns = gSystem->Getenv("OutputDir");
181  if (!fns.IsNull() && !fns.EndsWith("/")) fns += "/";
182  fns += chain->GetCurrentFile()->GetName();
183  Printf("Processing %d. file: %s", count, fns.Data());
184  fCurrentFileName = fns.Data();
185  }
186 
187  return kTRUE;
188 }
189 
190 //_____________________________________________________________________________
191 void AliAnalysisTaskFilteredTree::UserCreateOutputObjects()
192 {
193  // Create histograms
194  // Called once
195 
196  //
197  //get the output file to make sure the trees will be associated to it
198  OpenFile(1);
199  fTreeSRedirector = new TTreeSRedirector();
200 
201  //
202  // Create trees
203  fV0Tree = ((*fTreeSRedirector)<<"V0s").GetTree();
204  fHighPtTree = ((*fTreeSRedirector)<<"highPt").GetTree();
205  fdEdxTree = ((*fTreeSRedirector)<<"dEdx").GetTree();
206  fLaserTree = ((*fTreeSRedirector)<<"Laser").GetTree();
207  fMCEffTree = ((*fTreeSRedirector)<<"MCEffTree").GetTree();
208  fCosmicPairsTree = ((*fTreeSRedirector)<<"CosmicPairs").GetTree();
209 
210  if (!fDummyTrack) {
211  fDummyTrack=new AliESDtrack();
212  }
213 
214  // histogram booking
215 
216  Double_t minPt = 0.1;
217  Double_t maxPt = 100.;
218  Int_t nbinsPt = 30;
219 
220  Double_t logminPt = TMath::Log10(minPt);
221  Double_t logmaxPt = TMath::Log10(maxPt);
222  Double_t binwidth = (logmaxPt-logminPt)/nbinsPt;
223  Double_t *binsPt = new Double_t[nbinsPt+1];
224  binsPt[0] = minPt;
225  for (Int_t i=1;i<=nbinsPt;i++) {
226  binsPt[i] = minPt + TMath::Power(10,logminPt+i*binwidth);
227  }
228 
229  // 1pT resol cov matrix bins
230  Double_t min1PtRes = 0.;
231  Double_t max1PtRes = 0.3;
232  Int_t nbins1PtRes = 300;
233  Double_t bins1PtRes[301];
234  for (Int_t i=0;i<=nbins1PtRes;i++) {
235  bins1PtRes[i] = min1PtRes + i*(max1PtRes-min1PtRes)/nbins1PtRes;
236  }
237 
238  // phi bins
239  Double_t minPhi = 0.;
240  Double_t maxPhi = 6.5;
241  Int_t nbinsPhi = 100;
242  Double_t binsPhi[101];
243  for (Int_t i=0;i<=nbinsPhi;i++) {
244  binsPhi[i] = minPhi + i*(maxPhi-minPhi)/nbinsPhi;
245  }
246 
247  // eta bins
248  Double_t minEta = -1.;
249  Double_t maxEta = 1.;
250  Int_t nbinsEta = 20;
251  Double_t binsEta[21];
252  for (Int_t i=0;i<=nbinsEta;i++) {
253  binsEta[i] = minEta + i*(maxEta-minEta)/nbinsEta;
254  }
255 
256  // mult bins
257  Double_t minCent = 0.;
258  Double_t maxCent = 100;
259  Int_t nbinsCent = 20;
260  Double_t binsCent[101];
261  for (Int_t i=0;i<=nbinsCent;i++) {
262  binsCent[i] = minCent + i*(maxCent-minCent)/nbinsCent;
263  }
264 
265  fPtResPhiPtTPC = new TH3D("fPtResPhiPtTPC","pt rel. resolution from cov. matrix TPC tracks",nbinsPt,binsPt,nbinsPhi,binsPhi,nbins1PtRes,bins1PtRes);
266  fPtResPhiPtTPCc = new TH3D("fPtResPhiPtTPCc","pt rel. resolution from cov. matrix TPC constrained tracks",nbinsPt,binsPt,nbinsPhi,binsPhi,nbins1PtRes,bins1PtRes);
267  fPtResPhiPtTPCITS = new TH3D("fPtResPhiPtTPCITS","pt rel. resolution from cov. matrix TPC+ITS tracks",nbinsPt,binsPt,nbinsPhi,binsPhi,nbins1PtRes,bins1PtRes);
268 
269  fPtResEtaPtTPC = new TH3D("fPtResEtaPtTPC","pt rel. resolution from cov. matrix TPC tracks",nbinsPt,binsPt,nbinsEta,binsEta,nbins1PtRes,bins1PtRes);
270  fPtResEtaPtTPCc = new TH3D("fPtResEtaPtTPCc","pt rel. resolution from cov. matrix TPC constrained tracks",nbinsPt,binsPt,nbinsEta,binsEta,nbins1PtRes,bins1PtRes);
271  fPtResEtaPtTPCITS = new TH3D("fPtResEtaPtTPCITS","pt rel. resolution from cov. matrix TPC+ITS tracks",nbinsPt,binsPt,nbinsEta,binsEta,nbins1PtRes,bins1PtRes);
272 
273  fPtResCentPtTPC = new TH3D("fPtResCentPtTPC","pt rel. resolution from cov. matrix TPC tracks",nbinsPt,binsPt,nbinsCent,binsCent,nbins1PtRes,bins1PtRes);
274  fPtResCentPtTPCc = new TH3D("fPtResCentPtTPCc","pt rel. resolution from cov. matrix TPC constrained tracks",nbinsPt,binsPt,nbinsCent,binsCent,nbins1PtRes,bins1PtRes);
275  fPtResCentPtTPCITS = new TH3D("fPtResCentPtTPCITS","pt rel. resolution from cov. matrix TPC+ITS tracks",nbinsPt,binsPt,nbinsCent,binsCent,nbins1PtRes,bins1PtRes);
276 
277 
278  fOutput = new TList;
279  if(!fOutput) return;
280  fOutput->SetOwner();
281 
282  fOutput->Add(fPtResPhiPtTPC);
283  fOutput->Add(fPtResPhiPtTPCc);
284  fOutput->Add(fPtResPhiPtTPCITS);
285  fOutput->Add(fPtResEtaPtTPC);
286  fOutput->Add(fPtResEtaPtTPCc);
287  fOutput->Add(fPtResEtaPtTPCITS);
288  fOutput->Add(fPtResCentPtTPC);
289  fOutput->Add(fPtResCentPtTPCc);
290  fOutput->Add(fPtResCentPtTPCITS);
291 
292  // post data to outputs
293 
294  PostData(1,fV0Tree);
295  PostData(2,fHighPtTree);
296  PostData(3,fdEdxTree);
297  PostData(4,fLaserTree);
298  PostData(5,fMCEffTree);
299  PostData(6,fCosmicPairsTree);
300 
301  PostData(7,fOutput);
302 }
303 
304 //_____________________________________________________________________________
305 void AliAnalysisTaskFilteredTree::UserExec(Option_t *)
306 {
307  //
308  // Called for each event
309  //
310 
311  // ESD event
312  fESD = (AliESDEvent*) (InputEvent());
313  if (!fESD) {
314  Printf("ERROR: ESD event not available");
315  return;
316  }
317  //if MC info available - use it.
318  fMC = MCEvent();
319  if (fMC){
320  // Bug fix 28.05.2016 - do not trust to presence of MC handler, check if the content is valid
321  // - proper solution (autodetection of MC information) to be implemented
322  if (fMC->Stack()==NULL) {
323  fMC=NULL;
324  }else{
325  AliInfo("ToFix: MC stack not available. Prefered MCEvent() will return 0");
326  }
327  }
328 
329  if(fUseESDfriends) {
330  //fESDfriend = dynamic_cast<AliESDfriend*>(fESD->FindListObject("AliESDfriend"));
331  fESDfriend = ESDfriend();
332  if(!fESDfriend) {
333  Printf("ERROR: ESD friends not available");
334  }
335  }
336  AliInputEventHandler* inputHandler = (AliInputEventHandler*) AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler();
337  if (!inputHandler){
338  return;
339  }
340 
341  //if set, use the environment variables to set the downscaling factors
342  //AliAnalysisTaskFilteredTree_fLowPtTrackDownscaligF
343  //AliAnalysisTaskFilteredTree_fLowPtV0DownscaligF
344  //AliAnalysisTaskFilteredTree_fFriendDownscaling
345  TString env;
346  env = gSystem->Getenv("AliAnalysisTaskFilteredTree_fLowPtTrackDownscaligF");
347  if (!env.IsNull()){
348  fLowPtTrackDownscaligF=env.Atof();
349  AliInfo(Form("fLowPtTrackDownscaligF=%f",fLowPtTrackDownscaligF));
350  }
351  env = gSystem->Getenv("AliAnalysisTaskFilteredTree_fLowPtV0DownscaligF");
352  if (!env.IsNull()){
353  fLowPtV0DownscaligF=env.Atof();
354  AliInfo(Form("fLowPtV0DownscaligF=%f",fLowPtTrackDownscaligF));
355  }
356  env = gSystem->Getenv("AliAnalysisTaskFilteredTree_fFriendDownscaling");
357  if (!env.IsNull()){
358  fFriendDownscaling=env.Atof();
359  AliInfo(Form(" fFriendDownscaling=%f",fFriendDownscaling));
360  }
361  //
362  //
363  //
364  if(fProcessAll) {
365  ProcessAll(fESD,fMC,fESDfriend); // all track stages and MC
366  }
367  else {
368  Process(fESD,fMC,fESDfriend); // only global and TPC tracks
369  }
370  //
371  ProcessV0(fESD,fMC,fESDfriend);
372  ProcessLaser(fESD,fMC,fESDfriend);
373  ProcessdEdx(fESD,fMC,fESDfriend);
374  if (fProcessCosmics) { ProcessCosmics(fESD,fESDfriend); }
375  if(fMC) {
376  ProcessMCEff(fESD,fMC,fESDfriend);
377  //ProcessMC(); //TODO - enable MC detailed view switch after holidays
378  }
379  if (fProcessITSTPCmatchOut) ProcessITSTPCmatchOut(fESD, fESDfriend);
380  printf("processed event %d\n", Int_t(Entry()));
381 }
382 
383 //_____________________________________________________________________________
384 void AliAnalysisTaskFilteredTree::ProcessCosmics(AliESDEvent *const event, AliESDfriend* esdFriend)
385 {
386  //
387  // Find cosmic pairs (triggered or random)
388  //
389  //
390  AliESDVertex *vertexSPD = (AliESDVertex *)event->GetPrimaryVertexSPD();
391  AliESDVertex *vertexTPC = (AliESDVertex *)event->GetPrimaryVertexTPC();
392  const Double_t kMinPt=0.8;
393  const Double_t kMinPtMax=0.8;
394  const Double_t kMinNcl=50;
395  const Double_t kMaxDelta[5]={2,600,0.02,0.02,0.1};
396  Int_t ntracks=event->GetNumberOfTracks();
397  UInt_t specie = event->GetEventSpecie(); // skip laser events
398  if (specie==AliRecoParam::kCalib) return;
399  Int_t ntracksFriend = esdFriend ? esdFriend->GetNumberOfTracks() : 0;
400 
401 
402  for (Int_t itrack0=0;itrack0<ntracks;itrack0++) {
403  AliESDtrack *track0 = event->GetTrack(itrack0);
404  if (!track0) continue;
405  if (!track0->IsOn(AliESDtrack::kTPCrefit)) continue;
406 
407  if (TMath::Abs(AliTracker::GetBz())>1 && track0->Pt() < kMinPt) continue;
408  if (track0->Pt() < kMinPt) continue;
409  if (track0->GetTPCncls() < kMinNcl) continue;
410  if (TMath::Abs(track0->GetY())<kMaxDelta[0]) continue;
411  if (track0->GetKinkIndex(0)>0) continue;
412  const Double_t * par0=track0->GetParameter(); //track param at rhe DCA
413  //rm primaries
414  //
415  //track0->GetImpactParametersTPC(dcaTPC,covTPC);
416  //if (TMath::Abs(dcaTPC[0])<kMaxDelta[0]) continue;
417  //if (TMath::Abs(dcaTPC[1])<kMaxDelta[0]*2) continue;
418  // const AliExternalTrackParam * trackIn0 = track0->GetInnerParam();
419  AliESDfriendTrack* friendTrack0=NULL;
420  if (esdFriend &&!esdFriend->TestSkipBit()){
421  if (itrack0<ntracksFriend){
422  friendTrack0 = esdFriend->GetTrack(itrack0);
423  } //this guy can be NULL
424  }
425 
426  for (Int_t itrack1=itrack0+1;itrack1<ntracks;itrack1++) {
427  AliESDtrack *track1 = event->GetTrack(itrack1);
428  if (!track1) continue;
429  if (!track1->IsOn(AliESDtrack::kTPCrefit)) continue;
430  if (track1->GetKinkIndex(0)>0) continue;
431  if ((TMath::Abs(AliTracker::GetBz())>1) && (track1->Pt() < kMinPt)) continue;
432  if (track1->Pt() < kMinPt) continue;
433  if (track1->GetTPCncls()<kMinNcl) continue;
434  if (TMath::Abs(AliTracker::GetBz())>1 && TMath::Max(track1->Pt(), track0->Pt())<kMinPtMax) continue;
435  if (TMath::Abs(track1->GetY())<kMaxDelta[0]) continue;
436  //track1->GetImpactParametersTPC(dcaTPC,covTPC);
437  // if (TMath::Abs(dcaTPC[0])<kMaxDelta[0]) continue;
438  //if (TMath::Abs(dcaTPC[1])<kMaxDelta[0]*2) continue;
439  //
440  const Double_t* par1=track1->GetParameter(); //track param at rhe DCA
441  //
442  Bool_t isPair=kTRUE;
443  for (Int_t ipar=0; ipar<5; ipar++){
444  if (ipar==4&&TMath::Abs(AliTracker::GetBz())<1) continue; // 1/pt not defined for B field off
445  if (TMath::Abs(TMath::Abs(par0[ipar])-TMath::Abs(par1[ipar]))>kMaxDelta[ipar]) isPair=kFALSE;
446  }
447  if (!isPair) continue;
448  if (TMath::Abs(TMath::Abs(track0->GetAlpha()-track1->GetAlpha())-TMath::Pi())>kMaxDelta[2]) isPair=kFALSE;
449  //delta with correct sign
450  /*
451  TCut cut0="abs(t1.fP[0]+t0.fP[0])<2"
452  TCut cut3="abs(t1.fP[3]+t0.fP[3])<0.02"
453  TCut cut4="abs(t1.fP[4]+t0.fP[4])<0.2"
454  */
455  if (TMath::Abs(par0[0]+par1[0])>kMaxDelta[0]) isPair=kFALSE; //delta y opposite sign
456  if (TMath::Abs(par0[3]+par1[3])>kMaxDelta[3]) isPair=kFALSE; //delta tgl opposite sign
457  if (TMath::Abs(AliTracker::GetBz())>1 && TMath::Abs(par0[4]+par1[4])>kMaxDelta[4]) isPair=kFALSE; //delta 1/pt opposite sign
458  if (!isPair) continue;
459  TString filename(AliAnalysisManager::GetAnalysisManager()->GetTree()->GetCurrentFile()->GetName());
460  Int_t eventNumber = event->GetEventNumberInFile();
461  //
462  //
463  Int_t ntracksSPD = vertexSPD->GetNContributors();
464  Int_t ntracksTPC = vertexTPC->GetNContributors();
465  Int_t runNumber = event->GetRunNumber();
466  Int_t timeStamp = event->GetTimeStamp();
467  ULong64_t triggerMask = event->GetTriggerMask();
468  Float_t magField = event->GetMagneticField();
469  TObjString triggerClass = event->GetFiredTriggerClasses().Data();
470 
471  // Global event id calculation using orbitID, bunchCrossingID and periodID
472  ULong64_t orbitID = (ULong64_t)event->GetOrbitNumber();
473  ULong64_t bunchCrossID = (ULong64_t)event->GetBunchCrossNumber();
474  ULong64_t periodID = (ULong64_t)event->GetPeriodNumber();
475  ULong64_t gid = ((periodID << 36) | (orbitID << 12) | bunchCrossID);
476 
477 
478  AliESDfriendTrack* friendTrack1=NULL;
479  if (esdFriend &&!esdFriend->TestSkipBit()){
480  if (itrack1<ntracksFriend){
481  friendTrack1 = esdFriend->GetTrack(itrack1);
482  } //this guy can be NULL
483  }
484 
485  //
486  AliESDfriendTrack *friendTrackStore0=friendTrack0; // store friend track0 for later processing
487  AliESDfriendTrack *friendTrackStore1=friendTrack1; // store friend track1 for later processing
488  if (fFriendDownscaling>=1){ // downscaling number of friend tracks
489  if (gRandom->Rndm()>1./fFriendDownscaling){
490  friendTrackStore0 = 0;
491  friendTrackStore1 = 0;
492  }
493  }
494  if (fFriendDownscaling<=0){
495  if (((*fTreeSRedirector)<<"CosmicPairs").GetTree()){
496  TTree * tree = ((*fTreeSRedirector)<<"CosmicPairs").GetTree();
497  if (tree){
498  Double_t sizeAll=tree->GetZipBytes();
499  TBranch * br= tree->GetBranch("friendTrack0.fPoints");
500  Double_t sizeFriend=(br!=NULL)?br->GetZipBytes():0;
501  br= tree->GetBranch("friendTrack0.fCalibContainer");
502  if (br) sizeFriend+=br->GetZipBytes();
503  if (sizeFriend*TMath::Abs(fFriendDownscaling)>sizeAll) {
504  friendTrackStore0=0;
505  friendTrackStore1=0;
506  }
507  }
508  }
509  }
510  if(!fFillTree) return;
511  if(!fTreeSRedirector) return;
512  (*fTreeSRedirector)<<"CosmicPairs"<<
513  "gid="<<gid<< // global id of track
514  "fileName.="<<&fCurrentFileName<< // file name
515  "runNumber="<<runNumber<< // run number
516  "evtTimeStamp="<<timeStamp<< // time stamp of event
517  "evtNumberInFile="<<eventNumber<< // event number
518  "trigger="<<triggerMask<< // trigger mask
519  "triggerClass="<<&triggerClass<< // trigger class
520  "Bz="<<magField<< // magnetic field
521  //
522  "multSPD="<<ntracksSPD<< // event ultiplicity
523  "multTPC="<<ntracksTPC<< //
524  "vertSPD.="<<vertexSPD<< // primary vertex -SPD
525  "vertTPC.="<<vertexTPC<< // primary vertex -TPC
526  "t0.="<<track0<< // first half of comsic trak
527  "t1.="<<track1<< // second half of cosmic track
528  "friendTrack0.="<<friendTrackStore0<< // friend information first track + points
529  "friendTrack1.="<<friendTrackStore1<< // frined information first track + points
530  "\n";
531  }
532  }
533 }
534 
535 
536 //_____________________________________________________________________________
537 void AliAnalysisTaskFilteredTree::Process(AliESDEvent *const esdEvent, AliMCEvent * const mcEvent, AliESDfriend *const /*esdFriend*/)
538 {
539  //
540  // Select real events with high-pT tracks
541  //
542  static Int_t downscaleCounter=0;
543  // get selection cuts
544  AliFilteredTreeEventCuts *evtCuts = GetEventCuts();
545  AliFilteredTreeAcceptanceCuts *accCuts = GetAcceptanceCuts();
546  AliESDtrackCuts *esdTrackCuts = GetTrackCuts();
547 
548  if(!evtCuts || !accCuts || !esdTrackCuts) {
549  Printf("ERROR cuts not available");
550  return;
551  }
552 
553  // trigger selection
554  Bool_t isEventTriggered = kTRUE;
555  AliPhysicsSelection *physicsSelection = NULL;
556  AliTriggerAnalysis* triggerAnalysis = NULL;
557 
558  //
559  AliInputEventHandler* inputHandler = (AliInputEventHandler*) AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler();
560  // trigger
561  if(evtCuts->IsTriggerRequired())
562  {
563  // always MB
564  isEventTriggered = inputHandler->IsEventSelected() & AliVEvent::kMB;
565 
566  physicsSelection = static_cast<AliPhysicsSelection*> (inputHandler->GetEventSelection());
567  if(!physicsSelection) return;
568  //SetPhysicsTriggerSelection(physicsSelection);
569 
570  if (isEventTriggered && (GetTrigger() == AliTriggerAnalysis::kV0AND)) {
571  // set trigger (V0AND)
572  triggerAnalysis = physicsSelection->GetTriggerAnalysis();
573  if(!triggerAnalysis) return;
574  isEventTriggered = triggerAnalysis->IsOfflineTriggerFired(esdEvent, GetTrigger());
575  }
576  }
577 
578  // centrality determination
579  Float_t centralityF = -1;
580  AliCentrality *esdCentrality = esdEvent->GetCentrality();
581  centralityF = esdCentrality->GetCentralityPercentile(fCentralityEstimator.Data());
582 
583  // use MC information
584  AliHeader* header = 0;
585  AliGenEventHeader* genHeader = 0;
586  AliStack* stack = 0;
587  TArrayF vtxMC(3);
588 
589  Int_t multMCTrueTracks = 0;
590  if(mcEvent)
591  {
592  // get MC event header
593  header = mcEvent->Header();
594  if (!header) {
595  AliDebug(AliLog::kError, "Header not available");
596  return;
597  }
598  // MC particle stack
599  stack = mcEvent->Stack();
600  if (!stack) {
601  AliDebug(AliLog::kError, "Stack not available");
602  return;
603  }
604 
605  // get MC vertex
606  genHeader = header->GenEventHeader();
607  if (!genHeader) {
608  AliDebug(AliLog::kError, "Could not retrieve genHeader from Header");
609  return;
610  }
611  genHeader->PrimaryVertex(vtxMC);
612 
613  // multipliticy of all MC primary tracks
614  // in Zv, pt and eta ranges)
615  multMCTrueTracks = GetMCTrueTrackMult(mcEvent,evtCuts,accCuts);
616  } // end bUseMC
617 
618  // get reconstructed vertex
619  //const AliESDVertex* vtxESD = 0;
620  AliESDVertex* vtxESD = 0;
621  if(GetAnalysisMode() == kTPCAnalysisMode) {
622  vtxESD = (AliESDVertex*)esdEvent->GetPrimaryVertexTPC();
623  }
624  else if(GetAnalysisMode() == kTPCITSAnalysisMode) {
625  vtxESD = (AliESDVertex*)esdEvent->GetPrimaryVertexTracks();
626  }
627  else {
628  return;
629  }
630 
631  AliESDVertex* vtxTPC = (AliESDVertex*)esdEvent->GetPrimaryVertexTPC();
632  AliESDVertex* vtxSPD = (AliESDVertex*)esdEvent->GetPrimaryVertexSPD();
633 
634  if(!vtxESD) return;
635  if(!vtxTPC) return;
636  if(!vtxSPD) return;
637 
638  Bool_t isEventOK = evtCuts->AcceptEvent(esdEvent,mcEvent,vtxESD);
639  //printf("isEventOK %d, isEventTriggered %d, status %d, vz %f \n",isEventOK, isEventTriggered, vtxESD->GetStatus(), vtxESD->GetZ());
640  //printf("GetAnalysisMode() %d \n",GetAnalysisMode());
641  Int_t ntracks = esdEvent->GetNumberOfTracks();
642 
643  // check event cuts
644  if(isEventOK && isEventTriggered)
645  {
646 
647  //
648  // get IR information
649  //
650  AliESDHeader *esdHeader = 0;
651  esdHeader = esdEvent->GetHeader();
652  if(!esdHeader) return;
653  //Int_t ir1 = esdHeader->GetTriggerIREntries(); //all ir-s
654  //Int_t ir2 = esdHeader->GetTriggerIREntries(-1,1); // int2 set, 180 ms time interval
655 
656  // Use when Peter commit the changes in the header
657  Int_t ir1 = 0;
658  Int_t ir2 = 0;
659 
660  //
661  //Double_t vert[3] = {0};
662  //vert[0] = vtxESD->GetX();
663  //vert[1] = vtxESD->GetY();
664  //vert[2] = vtxESD->GetZ();
665  Int_t mult = vtxESD->GetNContributors();
666  Int_t multSPD = vtxSPD->GetNContributors();
667  Int_t multTPC = vtxTPC->GetNContributors();
668 
669  Float_t bz = esdEvent->GetMagneticField();
670  Int_t runNumber = esdEvent->GetRunNumber();
671  Int_t evtTimeStamp = esdEvent->GetTimeStamp();
672  Int_t evtNumberInFile = esdEvent->GetEventNumberInFile();
673 
674  // Global event id calculation using orbitID, bunchCrossingID and periodID
675  ULong64_t orbitID = (ULong64_t)esdEvent->GetOrbitNumber();
676  ULong64_t bunchCrossID = (ULong64_t)esdEvent->GetBunchCrossNumber();
677  ULong64_t periodID = (ULong64_t)esdEvent->GetPeriodNumber();
678  ULong64_t gid = ((periodID << 36) | (orbitID << 12) | bunchCrossID);
679 
680 
681  // high pT tracks
682  for (Int_t iTrack = 0; iTrack < esdEvent->GetNumberOfTracks(); iTrack++)
683  {
684  AliESDtrack *track = esdEvent->GetTrack(iTrack);
685  if(!track) continue;
686  if(track->Charge()==0) continue;
687  if(!esdTrackCuts->AcceptTrack(track)) continue;
688  if(!accCuts->AcceptTrack(track)) continue;
689 
690  // downscale low-pT tracks
691  Double_t scalempt= TMath::Min(track->Pt(),10.);
692  Double_t downscaleF = gRandom->Rndm();
693  downscaleF *= fLowPtTrackDownscaligF;
694  if( downscaleCounter>0 && TMath::Exp(2*scalempt)<downscaleF) continue;
695  //printf("TMath::Exp(2*scalempt) %e, downscaleF %e \n",TMath::Exp(2*scalempt), downscaleF);
696 
697  AliExternalTrackParam * tpcInner = (AliExternalTrackParam *)(track->GetTPCInnerParam());
698  if (!tpcInner) continue;
699  // transform to the track reference frame
700  Bool_t isOK = kFALSE;
701  isOK = tpcInner->Rotate(track->GetAlpha());
702  isOK = tpcInner->PropagateTo(track->GetX(),esdEvent->GetMagneticField());
703  if(!isOK) continue;
704 
705  // Dump to the tree
706  // vertex
707  // TPC-ITS tracks
708  //
709  TObjString triggerClass = esdEvent->GetFiredTriggerClasses().Data();
710  if(!fFillTree) return;
711  if(!fTreeSRedirector) return;
712  downscaleCounter++;
713  (*fTreeSRedirector)<<"highPt"<<
714  "gid="<<gid<<
715  "fileName.="<<&fCurrentFileName<<
716  "runNumber="<<runNumber<<
717  "evtTimeStamp="<<evtTimeStamp<<
718  "evtNumberInFile="<<evtNumberInFile<<
719  "triggerClass="<<&triggerClass<< // trigger
720  "Bz="<<bz<< // magnetic field
721  "vtxESD.="<<vtxESD<<
722  "ntracksESD="<<ntracks<< // number of tracks in the ESD
723  "IRtot="<<ir1<< // interaction record history info
724  "IRint2="<<ir2<<
725  "mult="<<mult<< // multiplicity of tracks pointing to the primary vertex
726  "multSPD="<<multSPD<< // multiplicity of tracks pointing to the SPD primary vertex
727  "multTPC="<<multTPC<< // multiplicity of tracks pointing to the TPC primary vertex
728  "esdTrack.="<<track<<
729  "centralityF="<<centralityF<<
730  "\n";
731  }
732  }
733 
734 }
735 
736 
737 //_____________________________________________________________________________
738 void AliAnalysisTaskFilteredTree::ProcessLaser(AliESDEvent *const esdEvent, AliMCEvent * const /*mcEvent*/, AliESDfriend *const esdFriend)
739 {
740  //
741  // Process laser events -> dump tracks and clusters to the special tree
742  //
743  const Double_t kMinPt = 5;
744  if(!fFillTree) return;
745  if(!fTreeSRedirector) return;
746  const AliESDHeader* esdHeader = esdEvent->GetHeader();
747  if(esdHeader && esdHeader->GetEventSpecie()==AliRecoParam::kCalib) {
748  Int_t countLaserTracks = 0;
749  Int_t runNumber = esdEvent->GetRunNumber();
750  Int_t evtTimeStamp = esdEvent->GetTimeStamp();
751  Int_t evtNumberInFile = esdEvent->GetEventNumberInFile();
752  Float_t bz = esdEvent->GetMagneticField();
753  TObjString triggerClass = esdEvent->GetFiredTriggerClasses().Data();
754  // Global event id calculation using orbitID, bunchCrossingID and periodID
755  ULong64_t orbitID = (ULong64_t)esdEvent->GetOrbitNumber();
756  ULong64_t bunchCrossID = (ULong64_t)esdEvent->GetBunchCrossNumber();
757  ULong64_t periodID = (ULong64_t)esdEvent->GetPeriodNumber();
758  ULong64_t gid = ((periodID << 36) | (orbitID << 12) | bunchCrossID);
759  for (Int_t iTrack = 0; iTrack < esdEvent->GetNumberOfTracks(); iTrack++){
760  AliESDtrack *track = esdEvent->GetTrack(iTrack);
761  if(!track) continue;
762  if (track->GetTPCInnerParam()==NULL) continue;
763  if(track->GetTPCInnerParam()) countLaserTracks++;
764  AliESDfriendTrack* friendTrack=NULL;
765  // suppress beam background and CE random reacks
766  if (track->GetInnerParam()->Pt()<kMinPt) continue;
767  Bool_t skipTrack=gRandom->Rndm()>1/(1+TMath::Abs(fFriendDownscaling));
768  if (skipTrack) continue;
769  if (esdFriend) {if (!esdFriend->TestSkipBit()) friendTrack = esdFriend->GetTrack(iTrack);} //this guy can be NULL
770  (*fTreeSRedirector)<<"Laser"<<
771  "gid="<<gid<< // global identifier of event
772  "fileName.="<<&fCurrentFileName<< //
773  "runNumber="<<runNumber<<
774  "evtTimeStamp="<<evtTimeStamp<<
775  "evtNumberInFile="<<evtNumberInFile<<
776  "triggerClass="<<&triggerClass<< // trigger
777  "Bz="<<bz<< // magnetic field
778  "multTPCtracks="<<countLaserTracks<< // multiplicity of tracks
779  "track.="<<track<< // track parameters
780  "friendTrack.="<<friendTrack<< // friend track information
781  "\n";
782  }
783  }
784 }
785 
786 //_____________________________________________________________________________
787 void AliAnalysisTaskFilteredTree::ProcessAll(AliESDEvent *const esdEvent, AliMCEvent * const mcEvent, AliESDfriend *const esdFriend)
788 {
789  //
790  // Select real events with high-pT tracks
791  // Calculate and stor in the output tree:
792  // TPC constrained tracks
793  // InnerParams constrained tracks
794  // TPC-ITS tracks
795  // ITSout-InnerParams tracks
796  // chi2 distance between TPC constrained and TPC-ITS tracks
797  // chi2 distance between TPC refitted constrained and TPC-ITS tracks
798  // chi2 distance between ITSout and InnerParams tracks
799  // MC information:
800  // track references at ITSin, TPCin; InnerParam at first TPC track reference,
801  // particle ID, mother ID, production mechanism ...
802  //
803  // get selection cuts
804  static Int_t downscaleCounter=0;
805  AliFilteredTreeEventCuts *evtCuts = GetEventCuts();
806  AliFilteredTreeAcceptanceCuts *accCuts = GetAcceptanceCuts();
807  AliESDtrackCuts *esdTrackCuts = GetTrackCuts();
808 
809  if(!evtCuts || !accCuts || !esdTrackCuts) {
810  AliDebug(AliLog::kError, "cuts not available");
811  return;
812  }
813 
814  // trigger selection
815  Bool_t isEventTriggered = kTRUE;
816  AliPhysicsSelection *physicsSelection = NULL;
817  AliTriggerAnalysis* triggerAnalysis = NULL;
818 
819  //
820  AliInputEventHandler* inputHandler = (AliInputEventHandler*) AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler();
821  AliPIDResponse *pidResponse = inputHandler->GetPIDResponse();
822 
823 
824  // trigger
825  if(evtCuts->IsTriggerRequired())
826  {
827  // always MB
828  isEventTriggered = inputHandler->IsEventSelected() & AliVEvent::kMB;
829 
830  physicsSelection = static_cast<AliPhysicsSelection*> (inputHandler->GetEventSelection());
831  if(!physicsSelection) {AliInfo("no physics selection"); return;}
832  //SetPhysicsTriggerSelection(physicsSelection);
833 
834  if (isEventTriggered && (GetTrigger() == AliTriggerAnalysis::kV0AND)) {
835  // set trigger (V0AND)
836  triggerAnalysis = physicsSelection->GetTriggerAnalysis();
837  if(!triggerAnalysis) {AliInfo("no trigger analysis");return;}
838  isEventTriggered = triggerAnalysis->IsOfflineTriggerFired(esdEvent, GetTrigger());
839  }
840  }
841 
842  // centrality determination
843  Float_t centralityF = -1;
844  AliCentrality *esdCentrality = esdEvent->GetCentrality();
845  centralityF = esdCentrality->GetCentralityPercentile(fCentralityEstimator.Data());
846 
847  // use MC information
848  AliHeader* header = 0;
849  AliGenEventHeader* genHeader = 0;
850  AliStack* stack = 0;
851  TArrayF vtxMC(3);
852  Int_t mcStackSize=0;
853 
854  Int_t multMCTrueTracks = 0;
855  if(mcEvent)
856  {
857  // get MC event header
858  header = mcEvent->Header();
859  if (!header) {
860  AliDebug(AliLog::kError, "Header not available");
861  return;
862  }
863  // MC particle stack
864  stack = mcEvent->Stack();
865  if (!stack) {
866  AliDebug(AliLog::kError, "Stack not available");
867  return;
868  }
869  mcStackSize=stack->GetNtrack();
870 
871  // get MC vertex
872  genHeader = header->GenEventHeader();
873  if (!genHeader) {
874  AliDebug(AliLog::kError, "Could not retrieve genHeader from Header");
875  return;
876  }
877  genHeader->PrimaryVertex(vtxMC);
878 
879  // multipliticy of all MC primary tracks
880  // in Zv, pt and eta ranges)
881  multMCTrueTracks = GetMCTrueTrackMult(mcEvent,evtCuts,accCuts);
882 
883  } // end bUseMC
884 
885  // get reconstructed vertex
886  //const AliESDVertex* vtxESD = 0;
887  AliESDVertex* vtxESD = 0;
888  if(GetAnalysisMode() == kTPCAnalysisMode) {
889  vtxESD = (AliESDVertex*)esdEvent->GetPrimaryVertexTPC();
890  }
891  else if(GetAnalysisMode() == kTPCITSAnalysisMode) {
892  vtxESD = (AliESDVertex*)esdEvent->GetPrimaryVertexTracks();
893  }
894  else {
895  AliInfo("no ESD vertex");
896  return;
897  }
898 
899  if(!vtxESD) return;
900 
901  Bool_t isEventOK = evtCuts->AcceptEvent(esdEvent,mcEvent,vtxESD);
902 
903  //
904  // Vertex info comparison and track multiplicity
905  //
906  AliESDVertex *vertexSPD = (AliESDVertex *)esdEvent->GetPrimaryVertexSPD();
907  AliESDVertex *vertexTPC = (AliESDVertex *)esdEvent->GetPrimaryVertexTPC();
908  Int_t contSPD = vertexSPD->GetNContributors();
909  Int_t contTPC = vertexTPC->GetNContributors();
910  TVectorD vertexPosTPC(3), vertexPosSPD(3);
911  vertexSPD->GetXYZ(vertexPosSPD.GetMatrixArray());
912  vertexTPC->GetXYZ(vertexPosTPC.GetMatrixArray());
913  Int_t ntracksTPC=0;
914  Int_t ntracksITS=0;
915  for (Int_t iTrack = 0; iTrack < esdEvent->GetNumberOfTracks(); iTrack++){
916  AliESDtrack *track = esdEvent->GetTrack(iTrack);
917  if(!track) continue;
918  if (track->IsOn(AliVTrack::kTPCrefit)) ntracksTPC++;
919  if (track->IsOn(AliVTrack::kITSrefit)) ntracksITS++;
920  }
921 
922  //printf("isEventOK %d, isEventTriggered %d \n",isEventOK, isEventTriggered);
923  //printf("GetAnalysisMode() %d \n",GetAnalysisMode());
924  Int_t ntracks = esdEvent->GetNumberOfTracks();
925  // Global event id calculation using orbitID, bunchCrossingID and periodID
926  ULong64_t orbitID = (ULong64_t)esdEvent->GetOrbitNumber();
927  ULong64_t bunchCrossID = (ULong64_t)esdEvent->GetBunchCrossNumber();
928  ULong64_t periodID = (ULong64_t)esdEvent->GetPeriodNumber();
929  ULong64_t gid = ((periodID << 36) | (orbitID << 12) | bunchCrossID);
930  TObjString triggerClass = esdEvent->GetFiredTriggerClasses().Data();
931  Float_t bz = esdEvent->GetMagneticField();
932  Int_t runNumber = esdEvent->GetRunNumber();
933  Int_t evtTimeStamp = esdEvent->GetTimeStamp();
934  Int_t evtNumberInFile = esdEvent->GetEventNumberInFile();
935  Int_t mult = vtxESD->GetNContributors();
936  (*fTreeSRedirector)<<"eventInfoTracks"<<
937  "gid="<<gid<<
938  "fileName.="<<&fCurrentFileName<< // name of the chunk file (hopefully full)
939  "runNumber="<<runNumber<< // runNumber
940  "evtTimeStamp="<<evtTimeStamp<< // time stamp of event (in seconds)
941  "evtNumberInFile="<<evtNumberInFile<< // event number
942  "triggerClass="<<&triggerClass<< // trigger class as a string
943  "Bz="<<bz<< // solenoid magnetic field in the z direction (in kGaus)
944  "mult="<<mult<< // multiplicity of tracks pointing to the primary vertex
945  "ntracks="<<ntracks<< // number of the esd tracks (to take into account the pileup in the TPC)
946  "isEventOK="<<isEventOK<< // flag - AliFilteredTreeEventCuts - track dumped only for selected events
947  "isEventTriggered="<<isEventTriggered<< // flag - if tigger required - track dumped only for selected events
948  "\n";
949 
950 
951 
952  // check event cuts
953  if(isEventOK && isEventTriggered)
954  {
955  //
956  // get IR information
957  //
958  AliESDHeader *esdHeader = 0;
959  esdHeader = esdEvent->GetHeader();
960  if(!esdHeader) {AliInfo("no esdHeader");return;}
961  //Int_t ir1 = esdHeader->GetTriggerIREntries(); //all ir-s
962  //Int_t ir2 = esdHeader->GetTriggerIREntries(-1,1); // int2 set, 180 ms time interval
963  //
964  Int_t ir1 = 0;
965  Int_t ir2 = 0;
966 
967  //
968  Double_t vert[3] = {0};
969  vert[0] = vtxESD->GetX();
970  vert[1] = vtxESD->GetY();
971  vert[2] = vtxESD->GetZ();
972  Int_t mult = vtxESD->GetNContributors();
973  Int_t numberOfTracks=esdEvent->GetNumberOfTracks();
974  // high pT tracks
975  for (Int_t iTrack = 0; iTrack < numberOfTracks; iTrack++)
976  {
977  AliESDtrack *track = esdEvent->GetTrack(iTrack);
978  AliESDfriendTrack* friendTrack=NULL;
979  Int_t numberOfFriendTracks=0;
980  if (esdFriend) numberOfFriendTracks=esdFriend->GetNumberOfTracks();
981  if (esdFriend && iTrack<numberOfFriendTracks) {if (!esdFriend->TestSkipBit()) friendTrack = esdFriend->GetTrack(iTrack);} //this guy can be NULL
982  if(!track) continue;
983  if(track->Charge()==0) continue;
984  if(!esdTrackCuts->AcceptTrack(track)) continue;
985  if(!accCuts->AcceptTrack(track)) continue;
986 
987  // downscale low-pT tracks
988  Double_t scalempt= TMath::Min(track->Pt(),10.);
989  Double_t downscaleF = gRandom->Rndm();
990  downscaleF *= fLowPtTrackDownscaligF;
991  if( downscaleCounter>0 && TMath::Exp(2*scalempt)<downscaleF) continue;
992  //printf("TMath::Exp(2*scalempt) %e, downscaleF %e \n",TMath::Exp(2*scalempt), downscaleF);
993 
994  // Dump to the tree
995  // vertex
996  // TPC constrained tracks
997  // InnerParams constrained tracks
998  // TPC-ITS tracks
999  // ITSout-InnerParams tracks
1000  // chi2 distance between TPC constrained and TPC-ITS tracks
1001  // chi2 distance between TPC refitted constrained and TPC-ITS tracks
1002  // chi2 distance between ITSout and InnerParams tracks
1003  // MC information
1004 
1005  Double_t x[3]; track->GetXYZ(x);
1006  Double_t b[3]; AliTracker::GetBxByBz(x,b);
1007 
1008  //
1009  // Transform TPC inner params to track reference frame
1010  //
1011  Bool_t isOKtpcInner = kFALSE;
1012  AliExternalTrackParam * tpcInner = (AliExternalTrackParam *)(track->GetTPCInnerParam());
1013  if (tpcInner) {
1014  // transform to the track reference frame
1015  isOKtpcInner = tpcInner->Rotate(track->GetAlpha());
1016  isOKtpcInner = tpcInner->PropagateTo(track->GetX(),esdEvent->GetMagneticField());
1017  }
1018 
1019  //
1020  // Constrain TPC inner to vertex
1021  // clone TPCinner has to be deleted
1022  //
1023  Bool_t isOKtpcInnerC = kFALSE;
1024  AliExternalTrackParam * tpcInnerC = new AliExternalTrackParam(*(track->GetTPCInnerParam()));
1025  if (tpcInnerC) {
1026  isOKtpcInnerC = ConstrainTPCInner(tpcInnerC,vtxESD,b);
1027  isOKtpcInnerC = tpcInnerC->Rotate(track->GetAlpha());
1028  isOKtpcInnerC = tpcInnerC->PropagateTo(track->GetX(),esdEvent->GetMagneticField());
1029  }
1030 
1031  //
1032  // Constrain TPC refitted tracks at inner TPC wall (InnerParams) to vertex
1033  // Clone track InnerParams has to be deleted
1034  //
1035  Bool_t isOKtrackInnerC = kTRUE;
1036  AliExternalTrackParam * trackInnerC = NULL;
1037  AliExternalTrackParam * trackInnerV = new AliExternalTrackParam(*(track->GetInnerParam()));
1038  isOKtrackInnerC=AliTracker::PropagateTrackToBxByBz(trackInnerV,3,track->GetMass(),3,kFALSE);
1039  isOKtrackInnerC&= trackInnerV->Rotate(track->GetAlpha());
1040  isOKtrackInnerC&= trackInnerV->PropagateTo(track->GetX(),esdEvent->GetMagneticField());
1041 
1042  if (isOKtrackInnerC) {
1043  trackInnerC = new AliExternalTrackParam(*trackInnerV);
1044  isOKtrackInnerC&= ConstrainTrackInner(trackInnerC,vtxESD,track->GetMass(),b);
1045  isOKtrackInnerC&= trackInnerC->Rotate(track->GetAlpha());
1046  isOKtrackInnerC&= trackInnerC->PropagateTo(track->GetX(),esdEvent->GetMagneticField());
1047  }
1048  //
1049  // calculate chi2 between vi and vj vectors
1050  // with covi and covj covariance matrices
1051  // chi2ij = (vi-vj)^(T)*(covi+covj)^(-1)*(vi-vj)
1052  //
1053  TMatrixD deltaT(5,1), deltaTtrackC(5,1);
1054  TMatrixD delta(1,5), deltatrackC(1,5);
1055  TMatrixD covarM(5,5), covarMtrackC(5,5);
1056  TMatrixD chi2(1,1);
1057  TMatrixD chi2trackC(1,1);
1058 
1059  if(isOKtpcInnerC && isOKtrackInnerC)
1060  {
1061  for (Int_t ipar=0; ipar<5; ipar++) {
1062  deltaT(ipar,0)=tpcInnerC->GetParameter()[ipar]-track->GetParameter()[ipar];
1063  delta(0,ipar)=tpcInnerC->GetParameter()[ipar]-track->GetParameter()[ipar];
1064 
1065  deltaTtrackC(ipar,0)=trackInnerC->GetParameter()[ipar]-track->GetParameter()[ipar];
1066  deltatrackC(0,ipar)=trackInnerC->GetParameter()[ipar]-track->GetParameter()[ipar];
1067 
1068  for (Int_t jpar=0; jpar<5; jpar++) {
1069  Int_t index=track->GetIndex(ipar,jpar);
1070  covarM(ipar,jpar)=track->GetCovariance()[index]+tpcInnerC->GetCovariance()[index];
1071  covarMtrackC(ipar,jpar)=track->GetCovariance()[index]+trackInnerC->GetCovariance()[index];
1072  }
1073  }
1074 
1075  // chi2 distance TPC constrained and TPC+ITS
1076  TMatrixD covarMInv = covarM.Invert();
1077  TMatrixD mat2 = covarMInv*deltaT;
1078  chi2 = delta*mat2;
1079  //chi2.Print();
1080 
1081  // chi2 distance TPC refitted constrained and TPC+ITS
1082  TMatrixD covarMInvtrackC = covarMtrackC.Invert();
1083  TMatrixD mat2trackC = covarMInvtrackC*deltaTtrackC;
1084  chi2trackC = deltatrackC*mat2trackC;
1085  //chi2trackC.Print();
1086  }
1087  //
1088  // Find nearest combined and ITS standaalone tracks
1089  AliExternalTrackParam paramITS; // nearest ITS track - chi2 distance at vertex
1090  AliExternalTrackParam paramITSC; // nearest ITS track - to constrained track chi2 distance at vertex
1091  AliExternalTrackParam paramComb; // nearest comb. tack - chi2 distance at inner wall
1092  Int_t indexNearestITS = GetNearestTrack((trackInnerV!=NULL)? trackInnerV:track, iTrack, esdEvent,0,0,paramITS);
1093  if (indexNearestITS<0) indexNearestITS = GetNearestTrack((trackInnerV!=NULL)? trackInnerV:track, iTrack, esdEvent,2,0,paramITS);
1094  Int_t indexNearestITSC = GetNearestTrack((trackInnerC!=NULL)? trackInnerC:track, iTrack, esdEvent,0,0,paramITSC);
1095  if (indexNearestITSC<0) indexNearestITS = GetNearestTrack((trackInnerC!=NULL)? trackInnerC:track, iTrack, esdEvent,2,0,paramITSC);
1096  Int_t indexNearestComb = GetNearestTrack(track->GetInnerParam(), iTrack, esdEvent,1,1,paramComb);
1097 
1098  //
1099  // Propagate ITSout to TPC inner wall
1100  // and calculate chi2 distance to track (InnerParams)
1101  //
1102  const Double_t kTPCRadius=85;
1103  const Double_t kStep=3;
1104 
1105  // clone track InnerParams has to be deleted
1106  Bool_t isOKtrackInnerC2 = kFALSE;
1107  AliExternalTrackParam *trackInnerC2 = new AliExternalTrackParam(*(track->GetInnerParam()));
1108  if (trackInnerC2) {
1109  isOKtrackInnerC2 = AliTracker::PropagateTrackToBxByBz(trackInnerC2,kTPCRadius,track->GetMass(),kStep,kFALSE);
1110  }
1111 
1112  Bool_t isOKouterITSc = kFALSE;
1113  AliExternalTrackParam *outerITSc = NULL;
1114  TMatrixD chi2OuterITS(1,1);
1115 
1116  if(esdFriend && !esdFriend->TestSkipBit())
1117  {
1118  // propagate ITSout to TPC inner wall
1119  if(friendTrack)
1120  {
1121 
1122  outerITSc = NULL;
1123  if (friendTrack->GetITSOut()) outerITSc = new AliExternalTrackParam(*(friendTrack->GetITSOut()));
1124  if(outerITSc)
1125  {
1126  isOKouterITSc = AliTracker::PropagateTrackToBxByBz(outerITSc,kTPCRadius,track->GetMass(),kStep,kFALSE);
1127  isOKouterITSc = outerITSc->Rotate(trackInnerC2->GetAlpha());
1128  isOKouterITSc = outerITSc->PropagateTo(trackInnerC2->GetX(),esdEvent->GetMagneticField());
1129 
1130  //
1131  // calculate chi2 between outerITS and innerParams
1132  // cov without z-coordinate at the moment
1133  //
1134  TMatrixD deltaTouterITS(4,1);
1135  TMatrixD deltaouterITS(1,4);
1136  TMatrixD covarMouterITS(4,4);
1137 
1138  if(isOKtrackInnerC2 && isOKouterITSc) {
1139  Int_t kipar = 0;
1140  Int_t kjpar = 0;
1141  for (Int_t ipar=0; ipar<5; ipar++) {
1142  if(ipar!=1) {
1143  deltaTouterITS(kipar,0)=outerITSc->GetParameter()[ipar]-trackInnerC2->GetParameter()[ipar];
1144  deltaouterITS(0,kipar)=outerITSc->GetParameter()[ipar]-trackInnerC2->GetParameter()[ipar];
1145  }
1146 
1147  kjpar=0;
1148  for (Int_t jpar=0; jpar<5; jpar++) {
1149  Int_t index=outerITSc->GetIndex(ipar,jpar);
1150  if(ipar !=1 || jpar!=1) {
1151  covarMouterITS(kipar,kjpar)=outerITSc->GetCovariance()[index]+trackInnerC2->GetCovariance()[index];
1152  }
1153  if(jpar!=1) kjpar++;
1154  }
1155  if(ipar!=1) kipar++;
1156  }
1157 
1158  // chi2 distance ITSout and InnerParams
1159  TMatrixD covarMInvouterITS = covarMouterITS.Invert();
1160  TMatrixD mat2outerITS = covarMInvouterITS*deltaTouterITS;
1161  chi2OuterITS = deltaouterITS*mat2outerITS;
1162  //chi2OuterITS.Print();
1163  }
1164  }
1165  }
1166  }
1167 
1168  //
1169  // MC info
1170  //
1171  TParticle *particle=NULL, *particleTPC=NULL, *particleITS=NULL;
1172  TParticle *particleMother=NULL, *particleMotherTPC=NULL, *particleMotherITS=NULL;
1173  Int_t mech=-1, mechTPC=-1, mechITS=-1;
1174  Bool_t isPrim=kFALSE, isPrimTPC=kFALSE, isPrimITS=kFALSE;
1175  Bool_t isFromStrangess=kFALSE, isFromStrangessTPC=kFALSE, isFromStrangessITS=kFALSE;
1176  Bool_t isFromConversion=kFALSE, isFromConversionTPC=kFALSE, isFromConversionITS=kFALSE;
1177  Bool_t isFromMaterial=kFALSE, isFromMaterialTPC=kFALSE, isFromMaterialITS=kFALSE;
1178 
1179  AliTrackReference *refTPCIn = NULL;
1180  AliTrackReference *refTPCOut = NULL;
1181  AliTrackReference *refITS = NULL;
1182  AliTrackReference *refTRD = NULL;
1183  AliTrackReference *refTOF = NULL;
1184  AliTrackReference *refEMCAL = NULL;
1185  AliTrackReference *refPHOS = NULL;
1186  Int_t nrefTPC=0, nrefTRD=0, nrefTOF=0, nrefITS=0, nrefEMCAL=0, nrefPHOS=0;
1187 
1188  Bool_t isOKtrackInnerC3 = kFALSE;
1189  AliExternalTrackParam *trackInnerC3 = new AliExternalTrackParam(*(track->GetInnerParam()));
1190  if(mcEvent && stack)
1191  {
1192  do //artificial loop (once) to make the continue statements jump out of the MC part
1193  {
1194  multMCTrueTracks = GetMCTrueTrackMult(mcEvent,evtCuts,accCuts);
1195  //
1196  // global track
1197  //
1198  Int_t label = TMath::Abs(track->GetLabel());
1199  if (label >= mcStackSize) continue;
1200  particle = stack->Particle(label);
1201  if (!particle) continue;
1202  if(particle && particle->GetPDG() && particle->GetPDG()->Charge()!=0.)
1203  {
1204  particleMother = GetMother(particle,stack);
1205  mech = particle->GetUniqueID();
1206  isPrim = stack->IsPhysicalPrimary(label);
1207  isFromStrangess = IsFromStrangeness(label,stack);
1208  isFromConversion = IsFromConversion(label,stack);
1209  isFromMaterial = IsFromMaterial(label,stack);
1210  }
1211 
1212  //
1213  // TPC track
1214  //
1215  Int_t labelTPC = TMath::Abs(track->GetTPCLabel());
1216  if (labelTPC >= mcStackSize) continue;
1217  particleTPC = stack->Particle(labelTPC);
1218  if (!particleTPC) continue;
1219  if(particleTPC && particleTPC->GetPDG() && particleTPC->GetPDG()->Charge()!=0.)
1220  {
1221  particleMotherTPC = GetMother(particleTPC,stack);
1222  mechTPC = particleTPC->GetUniqueID();
1223  isPrimTPC = stack->IsPhysicalPrimary(labelTPC);
1224  isFromStrangessTPC = IsFromStrangeness(labelTPC,stack);
1225  isFromConversionTPC = IsFromConversion(labelTPC,stack);
1226  isFromMaterialTPC = IsFromMaterial(labelTPC,stack);
1227  }
1228 
1229  //
1230  // store first track reference
1231  // for TPC track
1232  //
1233  TParticle *part=0;
1234  TClonesArray *trefs=0;
1235  Int_t status = mcEvent->GetParticleAndTR(TMath::Abs(labelTPC), part, trefs);
1236 
1237  if(status>0 && part && trefs && part->GetPDG() && part->GetPDG()->Charge()!=0.)
1238  {
1239  Int_t nTrackRef = trefs->GetEntries();
1240  //printf("nTrackRef %d \n",nTrackRef);
1241 
1242  Int_t countITS = 0;
1243  for (Int_t iref = 0; iref < nTrackRef; iref++)
1244  {
1245  AliTrackReference *ref = (AliTrackReference *)trefs->At(iref);
1246 
1247  // Ref. in the middle ITS
1248  if(ref && ref->Label()==label && ref->DetectorId()==AliTrackReference::kITS)
1249  {
1250  nrefITS++;
1251  if(!refITS && countITS==2) {
1252  refITS = ref;
1253  //printf("refITS %p \n",refITS);
1254  }
1255  countITS++;
1256  }
1257 
1258  // TPC
1259  if(ref && ref->Label()==label && ref->DetectorId()==AliTrackReference::kTPC)
1260  {
1261  nrefTPC++;
1262  refTPCOut=ref;
1263  if(!refTPCIn) {
1264  refTPCIn = ref;
1265  //printf("refTPCIn %p \n",refTPCIn);
1266  //break;
1267  }
1268  }
1269  // TRD
1270  if(ref && ref->Label()==label && ref->DetectorId()==AliTrackReference::kTRD)
1271  {
1272  nrefTRD++;
1273  if(!refTRD) {
1274  refTRD = ref;
1275  }
1276  }
1277  // TOF
1278  if(ref && ref->Label()==label && ref->DetectorId()==AliTrackReference::kTOF)
1279  {
1280  nrefTOF++;
1281  if(!refTOF) {
1282  refTOF = ref;
1283  }
1284  }
1285  // EMCAL
1286  if(ref && ref->Label()==label && ref->DetectorId()==AliTrackReference::kEMCAL)
1287  {
1288  nrefEMCAL++;
1289  if(!refEMCAL) {
1290  refEMCAL = ref;
1291  }
1292  }
1293  // PHOS
1294  // if(ref && ref->Label()==label && ref->DetectorId()==AliTrackReference::kPHOS)
1295  // {
1296  // nrefPHOS++;
1297  // if(!refPHOS) {
1298  // refPHOS = ref;
1299  // }
1300  // }
1301  }
1302 
1303  // transform inner params to TrackRef
1304  // reference frame
1305  if(refTPCIn && trackInnerC3)
1306  {
1307  Double_t kRefPhi = TMath::ATan2(refTPCIn->Y(),refTPCIn->X());
1308  isOKtrackInnerC3 = trackInnerC3->Rotate(kRefPhi);
1309  isOKtrackInnerC3 = AliTracker::PropagateTrackToBxByBz(trackInnerC3,refTPCIn->R(),track->GetMass(),kStep,kFALSE);
1310  }
1311  }
1312 
1313  //
1314  // ITS track
1315  //
1316  Int_t labelITS = TMath::Abs(track->GetITSLabel());
1317  if (labelITS >= mcStackSize) continue;
1318  particleITS = stack->Particle(labelITS);
1319  if (!particleITS) continue;
1320  if(particleITS && particleITS->GetPDG() && particleITS->GetPDG()->Charge()!=0.)
1321  {
1322  particleMotherITS = GetMother(particleITS,stack);
1323  mechITS = particleITS->GetUniqueID();
1324  isPrimITS = stack->IsPhysicalPrimary(labelITS);
1325  isFromStrangessITS = IsFromStrangeness(labelITS,stack);
1326  isFromConversionITS = IsFromConversion(labelITS,stack);
1327  isFromMaterialITS = IsFromMaterial(labelITS,stack);
1328  }
1329  }
1330  while (0);
1331  }
1332 
1333  //
1334  Bool_t dumpToTree=kFALSE;
1335 
1336  if(isOKtpcInnerC && isOKtrackInnerC) dumpToTree = kTRUE;
1337  //if(fUseESDfriends && isOKtrackInnerC2 && isOKouterITSc) dumpToTree = kTRUE;
1338  if(isOKtrackInnerC2 && isOKouterITSc) dumpToTree = kTRUE;
1339  if(mcEvent && isOKtrackInnerC3) dumpToTree = kTRUE;
1340  TObjString triggerClass = esdEvent->GetFiredTriggerClasses().Data();
1341  if (fReducePileUp){
1342  //
1343  // 18.03 - Reduce pile-up chunks, done outside of the ESDTrackCuts for 2012/2013 data pile-up about 95 % of tracks
1344  // Done only in case no MC info
1345  //
1346  Float_t dcaTPC[2];
1347  track->GetImpactParametersTPC(dcaTPC[0],dcaTPC[1]);
1348  Bool_t isRoughPrimary = TMath::Abs(dcaTPC[1])<10;
1349  Bool_t hasOuter=(track->IsOn(AliVTrack::kITSin))||(track->IsOn(AliVTrack::kTOFout))||(track->IsOn(AliVTrack::kTRDin));
1350  Bool_t keepPileUp=gRandom->Rndm()<0.05;
1351  if ( (!hasOuter) && (!isRoughPrimary) && (!keepPileUp)){
1352  dumpToTree=kFALSE;
1353  }
1354  }
1355 
1356  //init dummy objects
1357  static AliESDVertex dummyvtxESD;
1358  //if (!dummyvtxESD)
1359  //{
1360  // dummyvtxESD=new AliESDVertex();
1361  // //dummyvtxESD->SetNContributors(1);
1362  // //UShort_t pindices[1]; pindices[0]=0;
1363  // //dummyvtxESD->SetIndices(1,pindices);
1364  // //dummyvtxESD->SetNContributors(0);
1365  //}
1366  static AliExternalTrackParam dummyexternaltrackparam;
1367  //if (!dummyexternaltrackparam) dummyexternaltrackparam=new AliExternalTrackParam();
1368  static AliTrackReference dummytrackreference;
1369  //if (!dummytrackreference) dummytrackreference=new AliTrackReference();
1370  static TParticle dummyparticle;
1371  //if (!dummyparticle) dummyparticle=new TParticle();
1372 
1373  //assign the dummy objects if needed
1374  if (!track) {track=fDummyTrack;}
1375  AliESDfriendTrack *friendTrackStore=friendTrack; // store friend track for later processing
1376  if (fFriendDownscaling>=1){ // downscaling number of friend tracks
1377  friendTrackStore = (gRandom->Rndm()<1./fFriendDownscaling)? friendTrack:0;
1378  }
1379  if (fFriendDownscaling<=0){
1380  if (((*fTreeSRedirector)<<"highPt").GetTree()){
1381  TTree * tree = ((*fTreeSRedirector)<<"highPt").GetTree();
1382  if (tree){
1383  Double_t sizeAll=tree->GetZipBytes();
1384  TBranch * br= tree->GetBranch("friendTrack.fPoints");
1385  Double_t sizeFriend=(br!=NULL)?br->GetZipBytes():0;
1386  br= tree->GetBranch("friendTrack.fCalibContainer");
1387  if (br) sizeFriend+=br->GetZipBytes();
1388  if (sizeFriend*TMath::Abs(fFriendDownscaling)>sizeAll) friendTrackStore=0;
1389  }
1390  }
1391  }
1392 
1393 
1394  // if (!friendTrackStore && fFriendDownscaling<=1) {friendTrack=fDummyFriendTrack;}
1395  if (!vtxESD) {vtxESD=&dummyvtxESD;}
1396  if (mcEvent)
1397  {
1398  if (!refTPCIn) {refTPCIn=&dummytrackreference;}
1399  if (!refITS) {refITS=&dummytrackreference;}
1400  if (!particle) {particle=&dummyparticle;}
1401  if (!particleMother) {particleMother=&dummyparticle;}
1402  if (!particleTPC) {particleTPC=&dummyparticle;}
1403  if (!particleMotherTPC) {particleMotherTPC=&dummyparticle;}
1404  if (!particleITS) {particleITS=&dummyparticle;}
1405  if (!particleMotherITS) {particleMotherITS=&dummyparticle;}
1406  }
1407 
1408 
1409  // fill histograms
1410  FillHistograms(track, tpcInnerC, centralityF, (Double_t)chi2(0,0));
1411  TVectorD tofClInfo(5); // starting at 2014 - TOF infdo not part of the AliESDtrack
1412  tofClInfo[0]=track->GetTOFsignal();
1413  tofClInfo[1]=track->GetTOFsignalToT();
1414  tofClInfo[2]=track->GetTOFsignalRaw();
1415  tofClInfo[3]=track->GetTOFsignalDz();
1416  tofClInfo[4]=track->GetTOFsignalDx();
1417 
1418  //get the nSigma information; NB particle number ID in the vectors follow the convention of AliPID
1419  const Int_t nSpecies=AliPID::kSPECIES;
1420  TVectorD tpcNsigma(nSpecies);
1421  TVectorD tofNsigma(nSpecies);
1422  TVectorD tpcPID(nSpecies); // bayes
1423  TVectorD tofPID(nSpecies);
1424  if(pidResponse){
1425  for (Int_t ispecie=0; ispecie<nSpecies; ++ispecie) {
1426  if (ispecie == Int_t(AliPID::kMuon)) continue;
1427  tpcNsigma[ispecie] = pidResponse->NumberOfSigmas(AliPIDResponse::kTPC, track, (AliPID::EParticleType)ispecie);
1428  tofNsigma[ispecie] = pidResponse->NumberOfSigmas(AliPIDResponse::kTOF, track, (AliPID::EParticleType)ispecie);
1429  //
1430  }
1431  pidResponse->ComputePIDProbability(AliPIDResponse::kTPC, track, nSpecies, tpcPID.GetMatrixArray());
1432  pidResponse->ComputePIDProbability(AliPIDResponse::kTOF, track, nSpecies, tofPID.GetMatrixArray());
1433  }
1434  if(fTreeSRedirector && dumpToTree && fFillTree) {
1435  downscaleCounter++;
1436  (*fTreeSRedirector)<<"highPt"<<
1437  "downscaleCounter="<<downscaleCounter<<
1438  "gid="<<gid<<
1439  "fileName.="<<&fCurrentFileName<< // name of the chunk file (hopefully full)
1440  "runNumber="<<runNumber<< // runNumber
1441  "evtTimeStamp="<<evtTimeStamp<< // time stamp of event (in seconds)
1442  "evtNumberInFile="<<evtNumberInFile<< // event number
1443  "triggerClass="<<&triggerClass<< // trigger class as a string
1444  "Bz="<<bz<< // solenoid magnetic field in the z direction (in kGaus)
1445  "vtxESD.="<<vtxESD<< // vertexer ESD tracks (can be biased by TPC pileup tracks)
1446  "IRtot="<<ir1<< // interaction record (trigger) counters - coutner 1
1447  "IRint2="<<ir2<< // interaction record (trigger) coutners - counter 2
1448  "mult="<<mult<< // multiplicity of tracks pointing to the primary vertex
1449  "ntracks="<<ntracks<< // number of the esd tracks (to take into account the pileup in the TPC)
1450  // important variables for the pile-up studies
1451  "contTPC="<< contTPC<< // number of contributors to the TPC primary vertex candidate
1452  "contSPD="<< contSPD<< // number of contributors to the SPD primary vertex candidate
1453  "vertexPosTPC.="<<&vertexPosTPC<< // TPC vertex position
1454  "vertexPosSPD.="<<&vertexPosSPD<< // SPD vertex position
1455  "ntracksTPC="<<ntracksTPC<< // total number of the TPC tracks which were refitted
1456  "ntracksITS="<<ntracksITS<< // total number of the ITS tracks which were refitted
1457  //
1458  "esdTrack.="<<track<< // esdTrack as used in the physical analysis
1459  "tofClInfo.="<<&tofClInfo<< // tof info
1460  // "friendTrack.="<<friendTrack<< // esdFriendTrack associated to the esdTrack
1461  "tofNsigma.="<<&tofNsigma<<
1462  "tpcNsigma.="<<&tpcNsigma<<
1463  "tofPID.="<<&tofPID<< // bayesian PID - without priors
1464  "tpcPID.="<<&tpcPID<< // bayesian PID - without priors
1465 
1466  "friendTrack.="<<friendTrackStore<< // esdFriendTrack associated to the esdTrack
1467  "extTPCInnerC.="<<tpcInnerC<< // TPC track from the first tracking iteration propagated and updated at vertex
1468  "extInnerParamV.="<<trackInnerV<< // TPC+TRD inner param after refit propagate to vertex
1469  "extInnerParamC.="<<trackInnerC<< // TPC+TRD inner param after refit propagate and updated at vertex
1470  "extInnerParam.="<<trackInnerC2<< // TPC+TRD inner param after refit propagate to refernce TPC layer
1471  "extOuterITS.="<<outerITSc<< // ITS outer track propagated to the TPC refernce radius
1472  "extInnerParamRef.="<<trackInnerC3<< // TPC+TRD inner param after refit propagated to the first TPC reference
1473  "chi2TPCInnerC="<<chi2(0,0)<< // chi2 of tracks ???
1474  "chi2InnerC="<<chi2trackC(0,0)<< // chi2s of tracks TPCinner to the combined
1475  "chi2OuterITS="<<chi2OuterITS(0,0)<< // chi2s of tracks TPC at inner wall to the ITSout
1476  "centralityF="<<centralityF;
1477  // info for 2 track resolution studies and matching efficency studies
1478  //
1479  (*fTreeSRedirector)<<"highPt"<<
1480  "paramITS.="<<&paramITS<< // nearest ITS track - chi2 distance at vertex
1481  "paramITSC.="<<&paramITSC<< // nearest ITS track - to constrained track chi2 distance at vertex
1482  "paramComb.="<<&paramComb<< // nearest comb. tack - chi2 distance at inner wall
1483  "indexNearestITS="<<indexNearestITS<< // index of nearest ITS track
1484  "indexNearestITSC="<<indexNearestITSC<< // index of nearest ITS track for constrained track
1485  "indexNearestComb="<<indexNearestComb; // index of nearest track for constrained track
1486 
1487  if (mcEvent){
1488  static AliTrackReference refDummy;
1489  if (!refITS) refITS = &refDummy;
1490  if (!refTRD) refTRD = &refDummy;
1491  if (!refTOF) refTOF = &refDummy;
1492  if (!refEMCAL) refEMCAL = &refDummy;
1493  if (!refPHOS) refPHOS = &refDummy;
1494  downscaleCounter++;
1495  (*fTreeSRedirector)<<"highPt"<<
1496  "multMCTrueTracks="<<multMCTrueTracks<< // mC track multiplicities
1497  "nrefITS="<<nrefITS<< // number of track references in the ITS
1498  "nrefTPC="<<nrefTPC<< // number of track references in the TPC
1499  "nrefTRD="<<nrefTRD<< // number of track references in the TRD
1500  "nrefTOF="<<nrefTOF<< // number of track references in the TOF
1501  "nrefEMCAL="<<nrefEMCAL<< // number of track references in the TOF
1502  "nrefPHOS="<<nrefPHOS<< // number of track references in the TOF
1503  "refTPCIn.="<<refTPCIn<<
1504  "refTPCOut.="<<refTPCOut<<
1505  "refITS.="<<refITS<<
1506  "refTRD.="<<refTRD<<
1507  "refTOF.="<<refTOF<<
1508  "refEMCAL.="<<refEMCAL<<
1509  "refPHOS.="<<refPHOS<<
1510  "particle.="<<particle<<
1511  "particleMother.="<<particleMother<<
1512  "mech="<<mech<<
1513  "isPrim="<<isPrim<<
1514  "isFromStrangess="<<isFromStrangess<<
1515  "isFromConversion="<<isFromConversion<<
1516  "isFromMaterial="<<isFromMaterial<<
1517  "particleTPC.="<<particleTPC<<
1518  "particleMotherTPC.="<<particleMotherTPC<<
1519  "mechTPC="<<mechTPC<<
1520  "isPrimTPC="<<isPrimTPC<<
1521  "isFromStrangessTPC="<<isFromStrangessTPC<<
1522  "isFromConversionTPC="<<isFromConversionTPC<<
1523  "isFromMaterialTPC="<<isFromMaterialTPC<<
1524  "particleITS.="<<particleITS<<
1525  "particleMotherITS.="<<particleMotherITS<<
1526  "mechITS="<<mechITS<<
1527  "isPrimITS="<<isPrimITS<<
1528  "isFromStrangessITS="<<isFromStrangessITS<<
1529  "isFromConversionITS="<<isFromConversionITS<<
1530  "isFromMaterialITS="<<isFromMaterialITS;
1531  }
1532  //finish writing the entry
1533  AliInfo("writing tree highPt");
1534  (*fTreeSRedirector)<<"highPt"<<"\n";
1535  }
1536  AliSysInfo::AddStamp("filteringTask",iTrack,numberOfTracks,numberOfFriendTracks,(friendTrackStore)?0:1);
1537  delete tpcInnerC;
1538  delete trackInnerC;
1539  delete trackInnerC2;
1540  delete outerITSc;
1541  delete trackInnerC3;
1542  }
1543  }
1544 }
1545 
1546 
1547 //_____________________________________________________________________________
1548 void AliAnalysisTaskFilteredTree::ProcessMCEff(AliESDEvent *const esdEvent, AliMCEvent * const mcEvent, AliESDfriend *const /*esdFriend*/)
1549 {
1550  //
1551  // Fill tree for efficiency studies MC only
1552  static Int_t downscaleCounter=0;
1553  AliInfo("we start!");
1554  if(!mcEvent) {
1555  AliDebug(AliLog::kError, "mcEvent not available");
1556  return;
1557  }
1558 
1559  // get selection cuts
1560  AliFilteredTreeEventCuts *evtCuts = GetEventCuts();
1561  AliFilteredTreeAcceptanceCuts *accCuts = GetAcceptanceCuts();
1562  AliESDtrackCuts *esdTrackCuts = GetTrackCuts();
1563 
1564  if(!evtCuts || !accCuts || !esdTrackCuts) {
1565  AliDebug(AliLog::kError, "cuts not available");
1566  return;
1567  }
1568 
1569  // trigger selection
1570  Bool_t isEventTriggered = kTRUE;
1571  AliPhysicsSelection *physicsSelection = NULL;
1572  AliTriggerAnalysis* triggerAnalysis = NULL;
1573 
1574  //
1575  AliInputEventHandler* inputHandler = (AliInputEventHandler*) AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler();
1576 
1577 
1578  // trigger
1579  if(evtCuts->IsTriggerRequired())
1580  {
1581  // always MB
1582  isEventTriggered = inputHandler->IsEventSelected() & AliVEvent::kMB;
1583  physicsSelection = static_cast<AliPhysicsSelection*> (inputHandler->GetEventSelection());
1584  if(!physicsSelection) return;
1585 
1586  if (isEventTriggered && (GetTrigger() == AliTriggerAnalysis::kV0AND)) {
1587  // set trigger (V0AND)
1588  triggerAnalysis = physicsSelection->GetTriggerAnalysis();
1589  if(!triggerAnalysis) return;
1590  isEventTriggered = triggerAnalysis->IsOfflineTriggerFired(esdEvent, GetTrigger());
1591  }
1592  }
1593 
1594  // centrality determination
1595  Float_t centralityF = -1;
1596  AliCentrality *esdCentrality = esdEvent->GetCentrality();
1597  centralityF = esdCentrality->GetCentralityPercentile(fCentralityEstimator.Data());
1598 
1599  // use MC information
1600  AliHeader* header = 0;
1601  AliGenEventHeader* genHeader = 0;
1602  AliStack* stack = 0;
1603  Int_t mcStackSize=0;
1604  TArrayF vtxMC(3);
1605 
1606  Int_t multMCTrueTracks = 0;
1607  //
1608  if(!mcEvent) {
1609  AliDebug(AliLog::kError, "mcEvent not available");
1610  return;
1611  }
1612  // get MC event header
1613  header = mcEvent->Header();
1614  if (!header) {
1615  AliDebug(AliLog::kError, "Header not available");
1616  return;
1617  }
1618  // MC particle stack
1619  stack = mcEvent->Stack();
1620  if (!stack) {
1621  AliDebug(AliLog::kError, "Stack not available");
1622  return;
1623  }
1624  mcStackSize=stack->GetNtrack();
1625 
1626  // get MC vertex
1627  genHeader = header->GenEventHeader();
1628  if (!genHeader) {
1629  AliDebug(AliLog::kError, "Could not retrieve genHeader from Header");
1630  return;
1631  }
1632  genHeader->PrimaryVertex(vtxMC);
1633 
1634  // multiplicity of all MC primary tracks
1635  // in Zv, pt and eta ranges)
1636  multMCTrueTracks = GetMCTrueTrackMult(mcEvent,evtCuts,accCuts);
1637 
1638 
1639  // get reconstructed vertex
1640  //const AliESDVertex* vtxESD = 0;
1641  AliESDVertex* vtxESD = 0;
1642  if(GetAnalysisMode() == kTPCAnalysisMode) {
1643  vtxESD = (AliESDVertex*)esdEvent->GetPrimaryVertexTPC();
1644  }
1645  else if(GetAnalysisMode() == kTPCITSAnalysisMode) {
1646  vtxESD = (AliESDVertex*)esdEvent->GetPrimaryVertexTracks();
1647  }
1648  else {
1649  return;
1650  }
1651 
1652  if(!vtxESD) return;
1653  //
1654  // Vertex info comparison and track multiplicity
1655  //
1656  AliESDVertex *vertexSPD = (AliESDVertex *)esdEvent->GetPrimaryVertexSPD();
1657  AliESDVertex *vertexTPC = (AliESDVertex *)esdEvent->GetPrimaryVertexTPC();
1658  Int_t contSPD = vertexSPD->GetNContributors();
1659  Int_t contTPC = vertexTPC->GetNContributors();
1660  TVectorD vertexPosTPC(3), vertexPosSPD(3);
1661  vertexSPD->GetXYZ(vertexPosSPD.GetMatrixArray());
1662  vertexTPC->GetXYZ(vertexPosTPC.GetMatrixArray());
1663  Int_t ntracksTPC=0;
1664  Int_t ntracksITS=0;
1665  for (Int_t iTrack = 0; iTrack < esdEvent->GetNumberOfTracks(); iTrack++){
1666  AliESDtrack *track = esdEvent->GetTrack(iTrack);
1667  if(!track) continue;
1668  if (track->IsOn(AliVTrack::kTPCrefit)) ntracksTPC++;
1669  if (track->IsOn(AliVTrack::kITSrefit)) ntracksITS++;
1670  }
1671 
1672  Bool_t isEventOK = evtCuts->AcceptEvent(esdEvent,mcEvent,vtxESD);
1673 
1674  TObjString triggerClass = esdEvent->GetFiredTriggerClasses().Data();
1675 
1676  // check event cuts
1677  if(isEventOK && isEventTriggered)
1678  {
1679  //if(!stack) return;
1680 
1681  //
1682  // MC info
1683  //
1684  TParticle *particle=NULL;
1685  TParticle *particleMother=NULL;
1686  Int_t mech=-1;
1687 
1688  // reconstruction event info
1689  Double_t vert[3] = {0};
1690  vert[0] = vtxESD->GetX();
1691  vert[1] = vtxESD->GetY();
1692  vert[2] = vtxESD->GetZ();
1693  Int_t mult = vtxESD->GetNContributors();
1694  Double_t bz = esdEvent->GetMagneticField();
1695  Double_t runNumber = esdEvent->GetRunNumber();
1696  Double_t evtTimeStamp = esdEvent->GetTimeStamp();
1697  Int_t evtNumberInFile = esdEvent->GetEventNumberInFile();
1698  // loop over MC stack
1699  for (Int_t iMc = 0; iMc < mcStackSize; ++iMc)
1700  {
1701  particle = stack->Particle(iMc);
1702  if (!particle)
1703  continue;
1704 
1705  // only charged particles
1706  if(!particle->GetPDG()) continue;
1707  Double_t charge = particle->GetPDG()->Charge()/3.;
1708  if (TMath::Abs(charge) < 0.001)
1709  continue;
1710 
1711  // only primary particles
1712  Bool_t prim = stack->IsPhysicalPrimary(iMc);
1713  if(!prim) continue;
1714 
1715  // downscale low-pT particles
1716  Double_t scalempt= TMath::Min(particle->Pt(),10.);
1717  Double_t downscaleF = gRandom->Rndm();
1718  downscaleF *= fLowPtTrackDownscaligF;
1719  if (downscaleCounter>0 && TMath::Exp(2*scalempt)<downscaleF) continue;
1720  // is particle in acceptance
1721  if(!accCuts->AcceptTrack(particle)) continue;
1722 
1723  // check if particle reconstructed
1724  Bool_t isRec = kFALSE;
1725  Int_t trackIndex = -1;
1726  Int_t trackLoopIndex = -1;
1727  Int_t isESDtrackCut= 0;
1728  Int_t isAccCuts = 0;
1729  Int_t nRec = 0; // how many times reconstructed
1730  Int_t nFakes = 0; // how many times reconstructed as a fake track
1731  AliESDtrack *recTrack = NULL;
1732 
1733  for (Int_t iTrack = 0; iTrack < esdEvent->GetNumberOfTracks(); iTrack++)
1734  {
1735  AliESDtrack *track = esdEvent->GetTrack(iTrack);
1736  if(!track) continue;
1737  if(track->Charge()==0) continue;
1738  //
1739  Int_t label = TMath::Abs(track->GetLabel());
1740  if (label >= mcStackSize) continue;
1741  if(label == iMc) {
1742  Bool_t isAcc=esdTrackCuts->AcceptTrack(track);
1743  if (isAcc) isESDtrackCut=1;
1744  if (accCuts->AcceptTrack(track)) isAccCuts=1;
1745  isRec = kTRUE;
1746  trackIndex = iTrack;
1747 
1748  if (recTrack){
1749  if (track->GetTPCncls()<recTrack->GetTPCncls()) continue; // in case looper tracks use longer track
1750  if (!isAcc) continue;
1751  trackLoopIndex = iTrack;
1752  }
1753  recTrack = esdEvent->GetTrack(trackIndex);
1754  nRec++;
1755  if(track->GetLabel()<0) nFakes++;
1756 
1757  continue;
1758  }
1759  }
1760 
1761  // Store information in the output tree
1762  if (trackLoopIndex>-1) {
1763  recTrack = esdEvent->GetTrack(trackLoopIndex);
1764  } else if (trackIndex >-1) {
1765  recTrack = esdEvent->GetTrack(trackIndex);
1766  } else {
1767  recTrack = fDummyTrack;
1768  }
1769 
1770  particleMother = GetMother(particle,stack);
1771  mech = particle->GetUniqueID();
1772 
1773  //MC particle track length
1774  Double_t tpcTrackLength = 0.;
1775  AliMCParticle *mcParticle = (AliMCParticle*) mcEvent->GetTrack(iMc);
1776  if(mcParticle) {
1777  Int_t counter=0;
1778  tpcTrackLength = mcParticle->GetTPCTrackLength(bz,0.05,counter,3.0);
1779  }
1780 
1781 
1782  //
1783  if(fTreeSRedirector && fFillTree) {
1784  downscaleCounter++;
1785  (*fTreeSRedirector)<<"MCEffTree"<<
1786  "fileName.="<<&fCurrentFileName<<
1787  "triggerClass.="<<&triggerClass<<
1788  "runNumber="<<runNumber<<
1789  "evtTimeStamp="<<evtTimeStamp<<
1790  "evtNumberInFile="<<evtNumberInFile<< //
1791  "Bz="<<bz<< // magnetic field
1792  "vtxESD.="<<vtxESD<< // vertex info
1793  //
1794  "mult="<<mult<< // primary vertex 9whatewe found) multiplicity
1795  "multMCTrueTracks="<<multMCTrueTracks<< // mC track multiplicities
1796  // important variables for the pile-up studies
1797  "contTPC="<< contTPC<< // number of contributors to the TPC primary vertex candidate
1798  "contSPD="<< contSPD<< // number of contributors to the SPD primary vertex candidate
1799  "vertexPosTPC.="<<&vertexPosTPC<< // TPC vertex position
1800  "vertexPosSPD.="<<&vertexPosSPD<< // SPD vertex position
1801  "ntracksTPC="<<ntracksTPC<< // total number of the TPC tracks which were refitted
1802  "ntracksITS="<<ntracksITS<< // total number of the ITS tracks which were refitted
1803  //
1804  //
1805  "isAcc0="<<isESDtrackCut<< // track accepted by ESD track cuts
1806  "isAcc1="<<isAccCuts<< // track accepted by acceptance cuts flag
1807  "esdTrack.="<<recTrack<< // reconstructed track (only the longest from the loopers)
1808  "isRec="<<isRec<< // track was reconstructed
1809  "tpcTrackLength="<<tpcTrackLength<< // track length in the TPC r projection
1810  "particle.="<<particle<< // particle properties
1811  "particleMother.="<<particleMother<< // particle mother
1812  "mech="<<mech<< // production mechanizm
1813  "nRec="<<nRec<< // how many times reconstruted
1814  "nFakes="<<nFakes<< // how many times reconstructed as a fake track
1815  "\n";
1816  }
1817 
1818  //if(trackIndex <0 && recTrack) delete recTrack; recTrack=0;
1819  }
1820  }
1821 
1822 }
1823 
1824 //_____________________________________________________________________________
1825 Bool_t AliAnalysisTaskFilteredTree::IsHighDeDxParticle(AliESDtrack * track) {
1826  //
1827  // check if particle is Z > 1
1828  //
1829  if (track->GetTPCNcls() < 60) return kFALSE;
1830  Double_t mom = track->GetInnerParam()->GetP();
1831  if (mom < 0.2) return kFALSE; // protection against unexpected behavior of Aleph parameterization
1832  Float_t dca[2], bCov[3];
1833  track->GetImpactParameters(dca,bCov);
1834  //
1835 
1836  Double_t triggerDeDx = 4*AliExternalTrackParam::BetheBlochAleph((mom*2)/(0.938*3),1.0288,31.9806,5.04114e-11,2.13096,2.38541);
1837 
1838  if (track->GetTPCsignal() > triggerDeDx && track->GetTPCsignal()<1000 && TMath::Abs(dca[0])<3.) return kTRUE;
1839 
1840  return kFALSE;
1841 }
1842 
1843 //_____________________________________________________________________________
1844 void AliAnalysisTaskFilteredTree::ProcessV0(AliESDEvent *const esdEvent, AliMCEvent * const mcEvent, AliESDfriend *const esdFriend)
1845 {
1846  //
1847  // Select real events with V0 (K0s and Lambda and Gamma) high-pT candidates
1848  //
1849  static Int_t downscaleCounter=0;
1850  if(!esdEvent) {
1851  AliDebug(AliLog::kError, "esdEvent not available");
1852  return;
1853  }
1854 
1855  // get selection cuts
1856  AliFilteredTreeEventCuts *evtCuts = GetEventCuts();
1857  AliFilteredTreeAcceptanceCuts *accCuts = GetAcceptanceCuts();
1858  AliESDtrackCuts *esdTrackCuts = GetTrackCuts();
1859 
1860  if(!evtCuts || !accCuts || !esdTrackCuts) {
1861  AliDebug(AliLog::kError, "cuts not available");
1862  return;
1863  }
1864 
1865  // trigger selection
1866  Bool_t isEventTriggered = kTRUE;
1867  AliPhysicsSelection *physicsSelection = NULL;
1868  AliTriggerAnalysis* triggerAnalysis = NULL;
1869 
1870  //
1871  AliInputEventHandler* inputHandler = (AliInputEventHandler*) AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler();
1872  AliPIDResponse *pidResponse = inputHandler->GetPIDResponse();
1873 
1874  // trigger
1875  if(evtCuts->IsTriggerRequired())
1876  {
1877  // always MB
1878  isEventTriggered = inputHandler->IsEventSelected() & AliVEvent::kMB;
1879 
1880  physicsSelection = static_cast<AliPhysicsSelection*> (inputHandler->GetEventSelection());
1881  if(!physicsSelection) return;
1882  //SetPhysicsTriggerSelection(physicsSelection);
1883 
1884  if (isEventTriggered && (GetTrigger() == AliTriggerAnalysis::kV0AND)) {
1885  // set trigger (V0AND)
1886  triggerAnalysis = physicsSelection->GetTriggerAnalysis();
1887  if(!triggerAnalysis) return;
1888  isEventTriggered = triggerAnalysis->IsOfflineTriggerFired(esdEvent, GetTrigger());
1889  }
1890  }
1891 
1892  // centrality determination
1893  Float_t centralityF = -1;
1894  AliCentrality *esdCentrality = esdEvent->GetCentrality();
1895  centralityF = esdCentrality->GetCentralityPercentile(fCentralityEstimator.Data());
1896 
1897 
1898  // get reconstructed vertex
1899  //const AliESDVertex* vtxESD = 0;
1900  AliESDVertex* vtxESD = 0;
1901  if(GetAnalysisMode() == kTPCAnalysisMode) {
1902  vtxESD = (AliESDVertex*)esdEvent->GetPrimaryVertexTPC();
1903  }
1904  else if(GetAnalysisMode() == kTPCITSAnalysisMode) {
1905  vtxESD = (AliESDVertex*)esdEvent->GetPrimaryVertexTracks();
1906  }
1907  else {
1908  return;
1909  }
1910 
1911  if(!vtxESD) return;
1912 
1913  Bool_t isEventOK = evtCuts->AcceptEvent(esdEvent,mcEvent,vtxESD);
1914  //printf("isEventOK %d, isEventTriggered %d \n",isEventOK, isEventTriggered);
1915  //printf("GetAnalysisMode() %d \n",GetAnalysisMode());
1916  Int_t ntracks = esdEvent->GetNumberOfTracks();
1917  // Global event id calculation using orbitID, bunchCrossingID and periodID
1918  ULong64_t orbitID = (ULong64_t)esdEvent->GetOrbitNumber();
1919  ULong64_t bunchCrossID = (ULong64_t)esdEvent->GetBunchCrossNumber();
1920  ULong64_t periodID = (ULong64_t)esdEvent->GetPeriodNumber();
1921  ULong64_t gid = ((periodID << 36) | (orbitID << 12) | bunchCrossID);
1922  TObjString triggerClass = esdEvent->GetFiredTriggerClasses().Data();
1923  Float_t bz = esdEvent->GetMagneticField();
1924  Int_t run = esdEvent->GetRunNumber();
1925  Int_t time = esdEvent->GetTimeStamp();
1926  Int_t evtNumberInFile = esdEvent->GetEventNumberInFile();
1927  Int_t nV0s = esdEvent->GetNumberOfV0s();
1928  Int_t mult = vtxESD->GetNContributors();
1929  (*fTreeSRedirector)<<"eventInfoV0"<<
1930  "gid="<<gid<<
1931  "fileName.="<<&fCurrentFileName<< // name of the chunk file (hopefully full)
1932  "run="<<run<< // runNumber
1933  "time="<<time<< // time stamp of event (in seconds)
1934  "evtNumberInFile="<<evtNumberInFile<< // event number
1935  "triggerClass="<<&triggerClass<< // trigger class as a string
1936  "Bz="<<bz<< // solenoid magnetic field in the z direction (in kGaus)
1937  "mult="<<mult<< // multiplicity of tracks pointing to the primary vertex
1938  "ntracks="<<ntracks<< // number of the esd tracks (to take into account the pileup in the TPC)
1939  "nV0s="<<nV0s<< // number of V0s
1940  "isEventOK="<<isEventOK<< // flag - AliFilteredTreeEventCuts - track dumped only for selected events
1941  "isEventTriggered="<<isEventTriggered<< // flag - if tigger required - track dumped only for selected events
1942  "\n";
1943 
1944 
1945 
1946 
1947  // check event cuts
1948  if(isEventOK && isEventTriggered) {
1949  //
1950  // Dump the pt downscaled V0 into the tree
1951  //
1952  Int_t ntracks = esdEvent->GetNumberOfTracks();
1953  Int_t evNr=esdEvent->GetEventNumberInFile();
1954 
1955 
1956  for (Int_t iv0=0; iv0<nV0s; iv0++){
1957 
1958  AliESDv0 * v0 = esdEvent->GetV0(iv0);
1959  if (!v0) continue;
1960  AliESDtrack * track0 = esdEvent->GetTrack(v0->GetIndex(0));
1961  AliESDtrack * track1 = esdEvent->GetTrack(v0->GetIndex(1));
1962  if (!track0) continue;
1963  if (!track1) continue;
1964  AliESDfriendTrack* friendTrack0=NULL;
1965  AliESDfriendTrack* friendTrack1=NULL;
1966  if (esdFriend) {
1967  if (!esdFriend->TestSkipBit()){
1968  Int_t ntracksFriend = esdFriend->GetNumberOfTracks();
1969  if (v0->GetIndex(0)<ntracksFriend){
1970  friendTrack0 = esdFriend->GetTrack(v0->GetIndex(0)); //this guy can be NULL
1971  }
1972  if (v0->GetIndex(1)<ntracksFriend){
1973  friendTrack1 = esdFriend->GetTrack(v0->GetIndex(1)); //this guy can be NULL
1974  }
1975  }
1976  }
1977  if (track0->GetSign()<0) {
1978  track1 = esdEvent->GetTrack(v0->GetIndex(0));
1979  track0 = esdEvent->GetTrack(v0->GetIndex(1));
1980  }
1981 
1982  //
1983  AliESDfriendTrack *friendTrackStore0=friendTrack0; // store friend track0 for later processing
1984  AliESDfriendTrack *friendTrackStore1=friendTrack1; // store friend track1 for later processing
1985  if (fFriendDownscaling>=1){ // downscaling number of friend tracks
1986  if (gRandom->Rndm()>1./fFriendDownscaling){
1987  friendTrackStore0 = 0;
1988  friendTrackStore1 = 0;
1989  }
1990  }
1991  if (fFriendDownscaling<=0){
1992  if (((*fTreeSRedirector)<<"V0s").GetTree()){
1993  TTree * tree = ((*fTreeSRedirector)<<"V0s").GetTree();
1994  if (tree){
1995  Double_t sizeAll=tree->GetZipBytes();
1996  TBranch * br= tree->GetBranch("friendTrack0.fPoints");
1997  Double_t sizeFriend=(br!=NULL)?br->GetZipBytes():0;
1998  br= tree->GetBranch("friendTrack0.fCalibContainer");
1999  if (br) sizeFriend+=br->GetZipBytes();
2000  if (sizeFriend*TMath::Abs(fFriendDownscaling)>sizeAll) {
2001  friendTrackStore0=0;
2002  friendTrackStore1=0;
2003  }
2004  }
2005  }
2006  }
2007 
2008  //
2009  Bool_t isDownscaled = IsV0Downscaled(v0);
2010  if (downscaleCounter>0 && isDownscaled) continue;
2011  AliKFParticle kfparticle; //
2012  Int_t type=GetKFParticle(v0,esdEvent,kfparticle);
2013  if (type==0) continue;
2014  TObjString triggerClass = esdEvent->GetFiredTriggerClasses().Data();
2015 
2016  if(!fFillTree) return;
2017  if(!fTreeSRedirector) return;
2018 
2019  TVectorD tofClInfo0(5); // starting at 2014 - TOF infdo not part of the AliESDtrack
2020  TVectorD tofClInfo1(5); // starting at 2014 - TOF infdo not part of the AliESDtrack
2021  tofClInfo0[0]=track0->GetTOFsignal();
2022  tofClInfo0[1]=track0->GetTOFsignalToT();
2023  tofClInfo0[2]=track0->GetTOFsignalRaw();
2024  tofClInfo0[3]=track0->GetTOFsignalDz();
2025  tofClInfo0[4]=track0->GetTOFsignalDx();
2026  tofClInfo1[0]=track1->GetTOFsignal();
2027  tofClInfo1[1]=track1->GetTOFsignalToT();
2028  tofClInfo1[2]=track1->GetTOFsignalRaw();
2029  tofClInfo1[3]=track1->GetTOFsignalDz();
2030  tofClInfo1[4]=track1->GetTOFsignalDx();
2031 
2032  //get the nSigma information; NB particle number ID in the vectors follow the convention in AliPID
2033  const Int_t nSpecies=AliPID::kSPECIES;
2034  TVectorD tpcNsigma0(nSpecies);
2035  TVectorD tofNsigma0(nSpecies);
2036  TVectorD tpcNsigma1(nSpecies);
2037  TVectorD tofNsigma1(nSpecies);
2038  if(pidResponse){
2039  for (Int_t ispecie=0; ispecie<nSpecies; ++ispecie) {
2040  if (ispecie == Int_t(AliPID::kMuon)) continue;
2041  tpcNsigma0[ispecie] = pidResponse->NumberOfSigmas(AliPIDResponse::kTPC, track0, (AliPID::EParticleType)ispecie);
2042  tofNsigma0[ispecie] = pidResponse->NumberOfSigmas(AliPIDResponse::kTOF, track0, (AliPID::EParticleType)ispecie);
2043  tpcNsigma1[ispecie] = pidResponse->NumberOfSigmas(AliPIDResponse::kTPC, track1, (AliPID::EParticleType)ispecie);
2044  tofNsigma1[ispecie] = pidResponse->NumberOfSigmas(AliPIDResponse::kTOF, track1, (AliPID::EParticleType)ispecie);
2045  }
2046  }
2047 
2048  downscaleCounter++;
2049  (*fTreeSRedirector)<<"V0s"<<
2050  "gid="<<gid<< // global id of event
2051  "isDownscaled="<<isDownscaled<< //
2052  "triggerClass="<<&triggerClass<< // trigger
2053  "Bz="<<bz<< //
2054  "fileName.="<<&fCurrentFileName<< // full path - file name with ESD
2055  "runNumber="<<run<< //
2056  "evtTimeStamp="<<time<< // time stamp of event in secons
2057  "evtNumberInFile="<<evNr<< //
2058  "type="<<type<< // type of V0-
2059  "ntracks="<<ntracks<<
2060  "v0.="<<v0<<
2061  "kf.="<<&kfparticle<<
2062  "track0.="<<track0<< // track
2063  "track1.="<<track1<<
2064  "tofClInfo0.="<<&tofClInfo0<<
2065  "tofClInfo1.="<<&tofClInfo1<<
2066  "tofNsigma0.="<<&tofNsigma0<<
2067  "tofNsigma1.="<<&tofNsigma1<<
2068  "tpcNsigma0.="<<&tpcNsigma0<<
2069  "tpcNsigma1.="<<&tpcNsigma1<<
2070  "friendTrack0.="<<friendTrackStore0<<
2071  "friendTrack1.="<<friendTrackStore1<<
2072  "centralityF="<<centralityF<<
2073  "\n";
2074  }
2075  }
2076 }
2077 
2078 //_____________________________________________________________________________
2079 void AliAnalysisTaskFilteredTree::ProcessdEdx(AliESDEvent *const esdEvent, AliMCEvent * const mcEvent, AliESDfriend *const esdFriend)
2080 {
2081  //
2082  // Select real events with large TPC dEdx signal
2083  //
2084  static Int_t downscaleCounter=0;
2085  if(!esdEvent) {
2086  AliDebug(AliLog::kError, "esdEvent not available");
2087  return;
2088  }
2089 
2090  // get selection cuts
2091  AliFilteredTreeEventCuts *evtCuts = GetEventCuts();
2092  AliFilteredTreeAcceptanceCuts *accCuts = GetAcceptanceCuts();
2093  AliESDtrackCuts *esdTrackCuts = GetTrackCuts();
2094 
2095  if(!evtCuts || !accCuts || !esdTrackCuts) {
2096  AliDebug(AliLog::kError, "cuts not available");
2097  return;
2098  }
2099 
2100 
2101  // trigger
2102  Bool_t isEventTriggered = kTRUE;
2103  AliPhysicsSelection *physicsSelection = NULL;
2104  AliTriggerAnalysis* triggerAnalysis = NULL;
2105 
2106  //
2107  AliInputEventHandler* inputHandler = (AliInputEventHandler*) AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler();
2108  AliPIDResponse *pidResponse = inputHandler->GetPIDResponse();
2109 
2110  if(evtCuts->IsTriggerRequired())
2111  {
2112  // always MB
2113  isEventTriggered = inputHandler->IsEventSelected() & AliVEvent::kMB;
2114 
2115  physicsSelection = static_cast<AliPhysicsSelection*> (inputHandler->GetEventSelection());
2116  if(!physicsSelection) return;
2117 
2118  if (isEventTriggered && (GetTrigger() == AliTriggerAnalysis::kV0AND)) {
2119  // set trigger (V0AND)
2120  triggerAnalysis = physicsSelection->GetTriggerAnalysis();
2121  if(!triggerAnalysis) return;
2122  isEventTriggered = triggerAnalysis->IsOfflineTriggerFired(esdEvent, GetTrigger());
2123  }
2124  }
2125 
2126  // get reconstructed vertex
2127  AliESDVertex* vtxESD = 0;
2128  if(GetAnalysisMode() == kTPCAnalysisMode) {
2129  vtxESD = (AliESDVertex*)esdEvent->GetPrimaryVertexTPC();
2130  }
2131  else if(GetAnalysisMode() == kTPCITSAnalysisMode) {
2132  vtxESD = (AliESDVertex*)esdEvent->GetPrimaryVertexTracks();
2133  }
2134  else {
2135  return;
2136  }
2137  if(!vtxESD) return;
2138 
2139  Bool_t isEventOK = evtCuts->AcceptEvent(esdEvent,mcEvent,vtxESD);
2140  //printf("isEventOK %d, isEventTriggered %d \n",isEventOK, isEventTriggered);
2141  //printf("GetAnalysisMode() %d \n",GetAnalysisMode());
2142 
2143 
2144  // check event cuts
2145  if(isEventOK && isEventTriggered)
2146  {
2147  Double_t vert[3] = {0};
2148  vert[0] = vtxESD->GetX();
2149  vert[1] = vtxESD->GetY();
2150  vert[2] = vtxESD->GetZ();
2151  Int_t mult = vtxESD->GetNContributors();
2152  Double_t bz = esdEvent->GetMagneticField();
2153  Double_t runNumber = esdEvent->GetRunNumber();
2154  Double_t evtTimeStamp = esdEvent->GetTimeStamp();
2155  Int_t evtNumberInFile = esdEvent->GetEventNumberInFile();
2156 
2157  // Global event id calculation using orbitID, bunchCrossingID and periodID
2158  ULong64_t orbitID = (ULong64_t)esdEvent->GetOrbitNumber();
2159  ULong64_t bunchCrossID = (ULong64_t)esdEvent->GetBunchCrossNumber();
2160  ULong64_t periodID = (ULong64_t)esdEvent->GetPeriodNumber();
2161  ULong64_t gid = ((periodID << 36) | (orbitID << 12) | bunchCrossID);
2162 
2163  // large dEdx
2164  for (Int_t iTrack = 0; iTrack < esdEvent->GetNumberOfTracks(); iTrack++)
2165  {
2166  AliESDtrack *track = esdEvent->GetTrack(iTrack);
2167  if(!track) continue;
2168  AliESDfriendTrack* friendTrack=NULL;
2169  if (esdFriend && !esdFriend->TestSkipBit()) {
2170  Int_t ntracksFriend = esdFriend->GetNumberOfTracks();
2171  if (iTrack<ntracksFriend){
2172  friendTrack = esdFriend->GetTrack(iTrack);
2173  } //this guy can be NULL
2174  }
2175 
2176  if(track->Charge()==0) continue;
2177  if(!esdTrackCuts->AcceptTrack(track)) continue;
2178  if(!accCuts->AcceptTrack(track)) continue;
2179 
2180  if(!IsHighDeDxParticle(track)) continue;
2181  TObjString triggerClass = esdEvent->GetFiredTriggerClasses().Data();
2182 
2183  if(!fFillTree) return;
2184  if(!fTreeSRedirector) return;
2185 
2186 
2187  //get the nSigma information; NB particle number ID in the vectors follow the convention of AliPID
2188  const Int_t nSpecies=AliPID::kSPECIES;
2189  TVectorD tpcNsigma(nSpecies);
2190  TVectorD tofNsigma(nSpecies);
2191  if(pidResponse){
2192  for (Int_t ispecie=0; ispecie<nSpecies; ++ispecie) {
2193  if (ispecie == Int_t(AliPID::kMuon)) continue;
2194  tpcNsigma[ispecie] = pidResponse->NumberOfSigmas(AliPIDResponse::kTPC, track, (AliPID::EParticleType)ispecie);
2195  tofNsigma[ispecie] = pidResponse->NumberOfSigmas(AliPIDResponse::kTOF, track, (AliPID::EParticleType)ispecie);
2196  }
2197  }
2198 
2199  downscaleCounter++;
2200  (*fTreeSRedirector)<<"dEdx"<< // high dEdx tree
2201  "gid="<<gid<< // global id
2202  "fileName.="<<&fCurrentFileName<< // file name
2203  "runNumber="<<runNumber<<
2204  "evtTimeStamp="<<evtTimeStamp<<
2205  "evtNumberInFile="<<evtNumberInFile<<
2206  "triggerClass="<<&triggerClass<< // trigger
2207  "Bz="<<bz<<
2208  "vtxESD.="<<vtxESD<< //
2209  "mult="<<mult<<
2210  "esdTrack.="<<track<<
2211  "friendTrack.="<<friendTrack<<
2212  "tofNsigma.="<<&tofNsigma<<
2213  "tpcNsigma.="<<&tpcNsigma<<
2214  "\n";
2215  }
2216  }
2217 }
2218 
2219 //_____________________________________________________________________________
2220 Int_t AliAnalysisTaskFilteredTree::GetKFParticle(AliESDv0 *const v0, AliESDEvent * const event, AliKFParticle & kfparticle)
2221 {
2222  //
2223  // Create KF particle in case the V0 fullfill selection criteria
2224  //
2225  // Selection criteria
2226  // 0. algorithm cut
2227  // 1. track cut
2228  // 3. chi2 cut
2229  // 4. rough mass cut
2230  // 5. Normalized pointing angle cut
2231  //
2232  const Double_t cutMass=0.2;
2233  const Double_t kSigmaDCACut=3;
2234  //
2235  // 0.) algo cut - accept only on the fly
2236  //
2237  if (v0->GetOnFlyStatus() ==kFALSE) return 0;
2238  //
2239  // 1.) track cut
2240  //
2241  AliESDtrack * track0 = event->GetTrack(v0->GetIndex(0));
2242  AliESDtrack * track1 = event->GetTrack(v0->GetIndex(1));
2243  /*
2244  TCut cutD="abs(track0.fD/sqrt(track0.fCdd))>2&&abs(track1.fD/sqrt(track1.fCdd))>2";
2245  TCut cutTheta="abs(track0.fP[3])<1&&abs(track1.fP[3])<1";
2246  TCut cutNcl="track0.GetTPCClusterInfo(2,1)>100&&track1.GetTPCClusterInfo(2,1)>100";
2247  */
2248  if (TMath::Abs(track0->GetTgl())>1) return 0;
2249  if (TMath::Abs(track1->GetTgl())>1) return 0;
2250  if ((track0->GetTPCClusterInfo(2,1))<100) return 0;
2251  if ((track1->GetTPCClusterInfo(2,1))<100) return 0;
2252  Float_t pos0[2]={0}, cov0[3]={0};
2253  Float_t pos1[2]={0}, cov1[3]={0};
2254  track0->GetImpactParameters(pos0,cov0);
2255  track0->GetImpactParameters(pos1,cov1);
2256  //
2257  if (TMath::Abs(pos0[0])<kSigmaDCACut*TMath::Sqrt(cov0[0])) return 0;
2258  if (TMath::Abs(pos1[0])<kSigmaDCACut*TMath::Sqrt(cov1[0])) return 0;
2259  //
2260  //
2261  // 3.) Chi2 cut
2262  //
2263  Double_t chi2KF = v0->GetKFInfo(2,2,2);
2264  if (chi2KF>25) return 0;
2265  //
2266  // 4.) Rough mass cut - 0.200 GeV
2267  //
2268  static Double_t masses[2]={-1};
2269  if (masses[0]<0){
2270  masses[0] = TDatabasePDG::Instance()->GetParticle("K_S0")->Mass();
2271  masses[1] = TDatabasePDG::Instance()->GetParticle("Lambda0")->Mass();
2272  }
2273  Double_t mass00= v0->GetEffMass(0,0);
2274  Double_t mass22= v0->GetEffMass(2,2);
2275  Double_t mass42= v0->GetEffMass(4,2);
2276  Double_t mass24= v0->GetEffMass(2,4);
2277  Bool_t massOK=kFALSE;
2278  Int_t type=0;
2279  Int_t ptype=0;
2280  Double_t dmass=1;
2281  Int_t p1=0, p2=0;
2282  if (TMath::Abs(mass00-0)<cutMass) {
2283  massOK=kTRUE; type+=1;
2284  if (TMath::Abs(mass00-0)<dmass) {
2285  ptype=1;
2286  dmass=TMath::Abs(mass00-0);
2287  p1=0; p2=0;
2288  }
2289  }
2290  if (TMath::Abs(mass24-masses[1])<cutMass) {
2291  massOK=kTRUE; type+=2;
2292  if (TMath::Abs(mass24-masses[1])<dmass){
2293  dmass = TMath::Abs(mass24-masses[1]);
2294  ptype=2;
2295  p1=2; p2=4;
2296  }
2297  }
2298  if (TMath::Abs(mass42-masses[1])<cutMass) {
2299  massOK=kTRUE; type+=4;
2300  if (TMath::Abs(mass42-masses[1])<dmass){
2301  dmass = TMath::Abs(mass42-masses[1]);
2302  ptype=4;
2303  p1=4; p2=2;
2304  }
2305  }
2306  if (TMath::Abs(mass22-masses[0])<cutMass) {
2307  massOK=kTRUE; type+=8;
2308  if (TMath::Abs(mass22-masses[0])<dmass){
2309  dmass = TMath::Abs(mass22-masses[0]);
2310  ptype=8;
2311  p1=2; p2=2;
2312  }
2313  }
2314  if (type==0) return 0;
2315  //
2316  const Int_t spdg[5]={kPositron,kMuonPlus,kPiPlus, kKPlus, kProton};
2317  const AliExternalTrackParam *paramP = v0->GetParamP();
2318  const AliExternalTrackParam *paramN = v0->GetParamN();
2319  if (paramP->GetSign()<0){
2320  paramP=v0->GetParamP();
2321  paramN=v0->GetParamN();
2322  }
2323  //Double_t *pparam1 = (Double_t*)paramP->GetParameter();
2324  //Double_t *pparam2 = (Double_t*)paramN->GetParameter();
2325  //
2326  AliKFParticle kfp1( *paramP, spdg[p1] );
2327  AliKFParticle kfp2( *paramN, -1 *spdg[p2] );
2328  AliKFParticle V0KF;
2329  (V0KF)+=kfp1;
2330  (V0KF)+=kfp2;
2331  kfparticle=V0KF;
2332  //
2333  // Pointing angle
2334  //
2335  Double_t errPhi = V0KF.GetErrPhi();
2336  Double_t pointAngle= TMath::ACos(v0->GetV0CosineOfPointingAngle());
2337  if (pointAngle/errPhi>10) return 0;
2338  //
2339  return ptype;
2340 }
2341 
2342 //_____________________________________________________________________________
2343 Bool_t AliAnalysisTaskFilteredTree::IsV0Downscaled(AliESDv0 *const v0)
2344 {
2345  //
2346  // Downscale randomly low pt V0
2347  //
2348  //return kFALSE;
2349  Double_t maxPt= TMath::Max(v0->GetParamP()->Pt(), v0->GetParamN()->Pt());
2350  Double_t scalempt= TMath::Min(maxPt,10.);
2351  Double_t downscaleF = gRandom->Rndm();
2352  downscaleF *= fLowPtV0DownscaligF;
2353  //
2354  // Special treatment of the gamma conversion pt spectra is softer -
2355  Double_t mass00= v0->GetEffMass(0,0);
2356  const Double_t cutMass=0.2;
2357  if (TMath::Abs(mass00-0)<cutMass){
2358  downscaleF/=10.; // 10 times smaller downscaling for the gamma concersion candidate
2359  }
2360  //printf("V0 TMath::Exp(2*scalempt) %e, downscaleF %e \n",TMath::Exp(2*scalempt), downscaleF);
2361  if (TMath::Exp(2*scalempt)<downscaleF) return kTRUE;
2362  return kFALSE;
2363 
2364  /*
2365  TH1F his1("his1","his1",100,0,10);
2366  TH1F his2("his2","his2",100,0,10);
2367  {for (Int_t i=0; i<10000; i++){
2368  Double_t rnd=gRandom->Exp(1);
2369  Bool_t isDownscaled =TMath::Exp(rnd)<100*gRandom->Rndm();
2370  his1->Fill(rnd);
2371  if (!isDownscaled) his2->Fill(rnd);
2372  }}
2373 
2374 */
2375 
2376 }
2377 
2378 
2379 
2380 //_____________________________________________________________________________
2381 Bool_t AliAnalysisTaskFilteredTree::ConstrainTPCInner(AliExternalTrackParam *const tpcInnerC, const AliESDVertex* vtx, Double_t b[3])
2382 {
2383  // Constrain TPC inner params constrained
2384  //
2385  if(!tpcInnerC) return kFALSE;
2386  if(!vtx) return kFALSE;
2387 
2388  Double_t dz[2],cov[3];
2389  //AliESDVertex *vtx= (AliESDVertex *)esdEvent->GetPrimaryVertex();
2390  //if(!tpcInnerC->PropagateToDCA(vtx, esdEvent->GetMagneticField(), 3, dz, cov)) return kFALSE;
2391  //if(!tpcInnerC->PropagateToDCA(vtx, Bz, 3, dz, cov)) return kFALSE;
2392  if(!tpcInnerC->PropagateToDCABxByBz(vtx, b, 3, dz, cov)) return kFALSE;
2393 
2394 
2395  Double_t covar[6]; vtx->GetCovMatrix(covar);
2396  Double_t p[2]={tpcInnerC->GetParameter()[0]-dz[0],tpcInnerC->GetParameter()[1]-dz[1]};
2397  Double_t c[3]={covar[2],0.,covar[5]};
2398  Double_t chi2C=tpcInnerC->GetPredictedChi2(p,c);
2399  if (chi2C>kVeryBig) return kFALSE;
2400 
2401  if(!tpcInnerC->Update(p,c)) return kFALSE;
2402 
2403  return kTRUE;
2404 }
2405 
2406 //_____________________________________________________________________________
2407 Bool_t AliAnalysisTaskFilteredTree::ConstrainTrackInner(AliExternalTrackParam *const trackInnerC, const AliESDVertex* vtx, Double_t mass, Double_t b[3])
2408 {
2409  // Constrain TPC inner params constrained
2410  //
2411  if(!trackInnerC) return kFALSE;
2412  if(!vtx) return kFALSE;
2413 
2414  const Double_t kRadius = 2.8;
2415  const Double_t kMaxStep = 1.0;
2416 
2417  Double_t dz[2],cov[3];
2418 
2419  //AliESDVertex *vtx= (AliESDVertex *)esdEvent->GetPrimaryVertex();
2420  //if(!trackInnerC->PropagateToDCA(vtx, esdEvent->GetMagneticField(), 3, dz, cov)) return kFALSE;
2421  //if(!trackInnerC->PropagateToDCA(vtx, Bz, 3, dz, cov)) return kFALSE;
2422 
2423  if(!AliTracker::PropagateTrackToBxByBz(trackInnerC,kRadius,mass,kMaxStep,kFALSE)) return kFALSE;
2424  if(!trackInnerC->PropagateToDCABxByBz(vtx, b, 3, dz, cov)) return kFALSE;
2425 
2426  //
2427  Double_t covar[6]; vtx->GetCovMatrix(covar);
2428  Double_t p[2]={trackInnerC->GetParameter()[0]-dz[0],trackInnerC->GetParameter()[1]-dz[1]};
2429  Double_t c[3]={covar[2],0.,covar[5]};
2430  Double_t chi2C=trackInnerC->GetPredictedChi2(p,c);
2431  if (chi2C>kVeryBig) return kFALSE;
2432 
2433  if(!trackInnerC->Update(p,c)) return kFALSE;
2434 
2435  return kTRUE;
2436 }
2437 
2438 
2439 //_____________________________________________________________________________
2440 TParticle *AliAnalysisTaskFilteredTree::GetMother(TParticle *const particle, AliStack *const stack)
2441 {
2442  if(!particle) return NULL;
2443  if(!stack) return NULL;
2444 
2445  Int_t motherLabel = TMath::Abs(particle->GetMother(0));
2446  TParticle* mother = NULL;
2447  Int_t mcStackSize=stack->GetNtrack();
2448  if (motherLabel>=mcStackSize) return NULL;
2449  mother = stack->Particle(motherLabel);
2450 
2451  return mother;
2452 }
2453 
2454 //_____________________________________________________________________________
2455 Bool_t AliAnalysisTaskFilteredTree::IsFromConversion(Int_t label, AliStack *const stack)
2456 {
2457  Bool_t isFromConversion = kFALSE;
2458 
2459  if(stack) {
2460  Int_t mcStackSize=stack->GetNtrack();
2461  if (label>=mcStackSize) return kFALSE;
2462  TParticle* particle = stack->Particle(label);
2463  if (!particle) return kFALSE;
2464 
2465  if(particle && particle->GetPDG() && particle->GetPDG()->Charge()!=0)
2466  {
2467  Int_t mech = particle->GetUniqueID(); // production mechanism
2468  Bool_t isPrim = stack->IsPhysicalPrimary(label);
2469 
2470  Int_t motherLabel = TMath::Abs(particle->GetMother(0));
2471  if (motherLabel>=mcStackSize) return kFALSE;
2472  TParticle* mother = stack->Particle(motherLabel);
2473  if(mother) {
2474  Int_t motherPdg = mother->GetPdgCode();
2475 
2476  if(!isPrim && mech==5 && motherPdg==kGamma) {
2477  isFromConversion=kTRUE;
2478  }
2479  }
2480  }
2481  }
2482 
2483  return isFromConversion;
2484 }
2485 
2486 //_____________________________________________________________________________
2487 Bool_t AliAnalysisTaskFilteredTree::IsFromMaterial(Int_t label, AliStack *const stack)
2488 {
2489  Bool_t isFromMaterial = kFALSE;
2490 
2491  if(stack) {
2492  Int_t mcStackSize=stack->GetNtrack();
2493  if (label>=mcStackSize) return kFALSE;
2494  TParticle* particle = stack->Particle(label);
2495  if (!particle) return kFALSE;
2496 
2497  if(particle && particle->GetPDG() && particle->GetPDG()->Charge()!=0)
2498  {
2499  Int_t mech = particle->GetUniqueID(); // production mechanism
2500  Bool_t isPrim = stack->IsPhysicalPrimary(label);
2501 
2502  Int_t motherLabel = TMath::Abs(particle->GetMother(0));
2503  if (motherLabel>=mcStackSize) return kFALSE;
2504  TParticle* mother = stack->Particle(motherLabel);
2505  if(mother) {
2506  if(!isPrim && mech==13) {
2507  isFromMaterial=kTRUE;
2508  }
2509  }
2510  }
2511  }
2512 
2513  return isFromMaterial;
2514 }
2515 
2516 //_____________________________________________________________________________
2517 Bool_t AliAnalysisTaskFilteredTree::IsFromStrangeness(Int_t label, AliStack *const stack)
2518 {
2519  Bool_t isFromStrangeness = kFALSE;
2520 
2521  if(stack) {
2522  Int_t mcStackSize=stack->GetNtrack();
2523  if (label>=mcStackSize) return kFALSE;
2524  TParticle* particle = stack->Particle(label);
2525  if (!particle) return kFALSE;
2526 
2527  if(particle && particle->GetPDG() && particle->GetPDG()->Charge()!=0)
2528  {
2529  Int_t mech = particle->GetUniqueID(); // production mechanism
2530  Bool_t isPrim = stack->IsPhysicalPrimary(label);
2531 
2532  Int_t motherLabel = TMath::Abs(particle->GetMother(0));
2533  if (motherLabel>=mcStackSize) return kFALSE;
2534  TParticle* mother = stack->Particle(motherLabel);
2535  if(mother) {
2536  Int_t motherPdg = mother->GetPdgCode();
2537 
2538  // K+-, lambda, antilambda, K0s decays
2539  if(!isPrim && mech==4 &&
2540  (TMath::Abs(motherPdg)==kKPlus || TMath::Abs(motherPdg)==kLambda0 || motherPdg==kK0Short))
2541  {
2542  isFromStrangeness = kTRUE;
2543  }
2544  }
2545  }
2546  }
2547 
2548  return isFromStrangeness;
2549 }
2550 
2551 
2552 //_____________________________________________________________________________
2553 void AliAnalysisTaskFilteredTree::FinishTaskOutput()
2554 {
2555  //
2556  // Called one at the end
2557  // locally on working node
2558  //
2559  Bool_t deleteTrees=kTRUE;
2560  if ((AliAnalysisManager::GetAnalysisManager()))
2561  {
2562  if (AliAnalysisManager::GetAnalysisManager()->GetAnalysisType() ==
2563  AliAnalysisManager::kProofAnalysis)
2564  deleteTrees=kFALSE;
2565  }
2566  if (deleteTrees) delete fTreeSRedirector;
2567  fTreeSRedirector=NULL;
2568 }
2569 
2570 //_____________________________________________________________________________
2571 void AliAnalysisTaskFilteredTree::Terminate(Option_t *)
2572 {
2573  //
2574  // Called one at the end
2575  //
2576 }
2577 
2578 //_____________________________________________________________________________
2579 Int_t AliAnalysisTaskFilteredTree::GetMCTrueTrackMult(AliMCEvent *const mcEvent, AliFilteredTreeEventCuts *const evtCuts, AliFilteredTreeAcceptanceCuts *const accCuts)
2580 {
2581  //
2582  // calculate mc event true track multiplicity
2583  //
2584  if(!mcEvent) return 0;
2585 
2586  AliStack* stack = 0;
2587  Int_t mult = 0;
2588 
2589  // MC particle stack
2590  stack = mcEvent->Stack();
2591  if (!stack) return 0;
2592 
2593  //
2594  //printf("minZv %f, maxZv %f \n", evtCuts->GetMinZv(), evtCuts->GetMaxZv());
2595  //
2596 
2597  Bool_t isEventOK = evtCuts->AcceptMCEvent(mcEvent);
2598  if(!isEventOK) return 0;
2599 
2600  Int_t nPart = stack->GetNtrack();
2601  for (Int_t iMc = 0; iMc < nPart; ++iMc)
2602  {
2603  TParticle* particle = stack->Particle(iMc);
2604  if (!particle)
2605  continue;
2606 
2607  // only charged particles
2608  if(!particle->GetPDG()) continue;
2609  Double_t charge = particle->GetPDG()->Charge()/3.;
2610  if (TMath::Abs(charge) < 0.001)
2611  continue;
2612 
2613  // physical primary
2614  Bool_t prim = stack->IsPhysicalPrimary(iMc);
2615  if(!prim) continue;
2616 
2617  // checked accepted without pt cut
2618  //if(accCuts->AcceptTrack(particle))
2619  if( particle->Eta() > accCuts->GetMinEta() && particle->Eta() < accCuts->GetMaxEta() )
2620  {
2621  mult++;
2622  }
2623  }
2624 
2625  return mult;
2626 }
2627 
2628 //_____________________________________________________________________________
2629 void AliAnalysisTaskFilteredTree::FillHistograms(AliESDtrack* const ptrack, AliExternalTrackParam* const ptpcInnerC, Double_t centralityF, Double_t chi2TPCInnerC)
2630 {
2631  //
2632  // Fill pT relative resolution histograms for
2633  // TPC only, TPC only constrained to vertex and TPC+ITS tracking
2634  //
2635  if(!ptrack) return;
2636  if(!ptpcInnerC) return;
2637 
2638  const AliExternalTrackParam * innerParam = (AliExternalTrackParam *) ptrack->GetInnerParam();
2639  if(!innerParam) return;
2640 
2641  Float_t dxy, dz;
2642  ptrack->GetImpactParameters(dxy,dz);
2643 
2644  // TPC+ITS primary tracks
2645  if( abs(ptrack->Eta())<0.8 &&
2646  ptrack->GetTPCClusterInfo(3,1)>120 &&
2647  ptrack->IsOn(0x40) &&
2648  ptrack->GetTPCclusters(0)>0.0 &&
2649  ptrack->GetTPCnclsS()/ptrack->GetTPCclusters(0)<0.4 &&
2650  //abs(innerParam->GetX())>0.0 &&
2651  //abs(innerParam->GetY()/innerParam->GetX())<0.14 &&
2652  //abs(innerParam->GetTgl())<0.85 &&
2653  ptrack->IsOn(0x0004) &&
2654  ptrack->GetNcls(0)>0 &&
2655  ptrack->GetITSchi2()>0 &&
2656  sqrt(ptrack->GetITSchi2()/ptrack->GetNcls(0))<6 &&
2657  sqrt(chi2TPCInnerC)<6 &&
2658  (ptrack->HasPointOnITSLayer(0) || ptrack->HasPointOnITSLayer(1)) &&
2659  abs(dz)<2.0 &&
2660  abs(dxy)<(0.018+0.035*abs(ptrack->GetSigned1Pt())) )
2661  {
2662  fPtResPhiPtTPCITS->Fill(ptrack->Pt(),ptrack->Phi(),1./abs(ptrack->GetSigned1Pt())*TMath::Sqrt(ptrack->GetSigma1Pt2()));
2663  fPtResEtaPtTPCITS->Fill(ptrack->Pt(),ptrack->Eta(),1./abs(ptrack->GetSigned1Pt())*TMath::Sqrt(ptrack->GetSigma1Pt2()));
2664  fPtResCentPtTPCITS->Fill(ptrack->Pt(),centralityF,1./abs(ptrack->GetSigned1Pt())*TMath::Sqrt(ptrack->GetSigma1Pt2()));
2665  }
2666 
2667  // TPC primary tracks
2668  // and TPC constrained primary tracks
2669 
2670  AliExternalTrackParam *ptpcInner = (AliExternalTrackParam *) ptrack->GetTPCInnerParam();
2671  if(!ptpcInner) return;
2672 
2673 
2674  Float_t dxyTPC, dzTPC;
2675  ptrack->GetImpactParametersTPC(dxyTPC,dzTPC);
2676 
2677  if( abs(ptrack->Eta())<0.8 &&
2678  ptrack->GetTPCClusterInfo(3,1)>120 &&
2679  ptrack->IsOn(0x40)&&
2680  ptrack->GetTPCclusters(0)>0.0 &&
2681  ptrack->GetTPCnclsS()/ptrack->GetTPCclusters(0)<0.4 &&
2682  //abs(innerParam->GetX())>0.0 &&
2683  //abs(innerParam->GetY()/innerParam->GetX())<0.14 &&
2684  //abs(innerParam->GetTgl())<0.85 &&
2685  abs(dzTPC)<3.2 &&
2686  abs(dxyTPC)<2.4 )
2687  {
2688  // TPC only
2689  fPtResPhiPtTPC->Fill(ptpcInner->Pt(),ptpcInner->Phi(),1./abs(ptpcInner->GetSigned1Pt())*TMath::Sqrt(ptpcInner->GetSigma1Pt2()));
2690  fPtResEtaPtTPC->Fill(ptpcInner->Pt(),ptpcInner->Eta(),1./abs(ptpcInner->GetSigned1Pt())*TMath::Sqrt(ptpcInner->GetSigma1Pt2()));
2691  fPtResCentPtTPC->Fill(ptpcInner->Pt(),centralityF,1./abs(ptpcInner->GetSigned1Pt())*TMath::Sqrt(ptpcInner->GetSigma1Pt2()));
2692 
2693  // TPC constrained to vertex
2694  fPtResPhiPtTPCc->Fill(ptpcInnerC->Pt(),ptpcInnerC->Phi(),1./abs(ptpcInnerC->GetSigned1Pt())*TMath::Sqrt(ptpcInnerC->GetSigma1Pt2()));
2695  fPtResEtaPtTPCc->Fill(ptpcInnerC->Pt(),ptpcInnerC->Eta(),1./abs(ptpcInnerC->GetSigned1Pt())*TMath::Sqrt(ptpcInnerC->GetSigma1Pt2()));
2696  fPtResCentPtTPCc->Fill(ptpcInnerC->Pt(),centralityF,1./abs(ptpcInnerC->GetSigned1Pt())*TMath::Sqrt(ptpcInnerC->GetSigma1Pt2()));
2697  }
2698 }
2699 
2700 
2701 void AliAnalysisTaskFilteredTree::ProcessITSTPCmatchOut(AliESDEvent *const esdEvent, AliESDfriend *const esdFriend){
2702  //
2703  // Process ITS standalone tracks find match with closest TPC(or combined tracks) tracks
2704  // marian.ivanov@cern.ch
2705  // 0.) Init variables
2706  // 1.) GetTrack parameters at TPC inner wall
2707  // 2.) Match closest TPC track (STANDALONE/global) - chi2 match criteria
2708  //
2709  // Logic to be used in reco:
2710  // 1.) Find matching ITSalone->TPCalone
2711  // 2.) if (!TPCalone.FindClose(TPCother)) TPCalone.Addopt(ITSalone)
2712  // 3.) ff ((ITSalone.FindClose(Global)==0) CreateGlobaltrack
2713  const Double_t radiusMatch=84.; // redius to propagate
2714  //
2715  const Double_t dFastPhiCut=0.2; // 6 sigma (200 MeV) fast angular cut
2716  const Double_t dFastThetaCut=0.12; // 6 sigma (200 MeV) fast angular cut
2717  const Double_t dFastPosPhiCut=0.06; // 6 sigma (200 MeV) fast angular cut
2718  const Double_t dFastZCut=6; // 6 sigma (200 MeV) fast z difference cut
2719  const Double_t dFastPtCut=2.; // 6 sigma (200 MeV) fast 1/pt cut
2720  const Double_t chi2Cut=100; // chi2 matching cut
2721  //
2722  if (!esdFriend) return; // not ITS standalone track
2723  if (esdFriend->TestSkipBit()) return; // friends tracks not stored
2724  Int_t ntracks=esdEvent->GetNumberOfTracks();
2725  Float_t bz = esdEvent->GetMagneticField();
2726  //
2727  // 0.) Get parameters in reference radius TPC Inner wall
2728  //
2729  //
2730  TMatrixD vecPosR0(ntracks,6); // possition and momentum estimate at reference radius
2731  TMatrixD vecMomR0(ntracks,6); //
2732  TMatrixD vecPosR1(ntracks,6); // possition and momentum estimate at reference radius TPC track
2733  TMatrixD vecMomR1(ntracks,6); //
2734  Double_t xyz[3], pxyz[3]; //
2735  for (Int_t iTrack=0; iTrack<ntracks; iTrack++){
2736  AliESDtrack *track = esdEvent->GetTrack(iTrack);
2737  if(!track) continue;
2738  if (track->GetInnerParam()){
2739  const AliExternalTrackParam *trackTPC=track->GetInnerParam();
2740  trackTPC->GetXYZAt(radiusMatch,bz,xyz);
2741  trackTPC->GetPxPyPzAt(radiusMatch,bz,pxyz);
2742  for (Int_t i=0; i<3; i++){
2743  vecPosR1(iTrack,i)=xyz[i];
2744  vecMomR1(iTrack,i)=pxyz[i];
2745  }
2746  vecPosR1(iTrack,3)= TMath::ATan2(xyz[1],xyz[0]); // phi pos angle
2747  vecMomR1(iTrack,3)= TMath::ATan2(pxyz[1],pxyz[0]); // phi mom angle
2748  vecMomR1(iTrack,4)= trackTPC->GetSigned1Pt();;
2749  vecMomR1(iTrack,5)= trackTPC->GetTgl();;
2750  }
2751  AliESDfriendTrack* friendTrack=esdFriend->GetTrack(iTrack);
2752  if(!friendTrack) continue;
2753  if (friendTrack->GetITSOut()){
2754  const AliExternalTrackParam *trackITS=friendTrack->GetITSOut();
2755  trackITS->GetXYZAt(radiusMatch,bz,xyz);
2756  trackITS->GetPxPyPzAt(radiusMatch,bz,pxyz);
2757  for (Int_t i=0; i<3; i++){
2758  vecPosR0(iTrack,i)=xyz[i];
2759  vecMomR0(iTrack,i)=pxyz[i];
2760  }
2761  vecPosR0(iTrack,3)= TMath::ATan2(xyz[1],xyz[0]);
2762  vecMomR0(iTrack,3)= TMath::ATan2(pxyz[1],pxyz[0]);
2763  vecMomR0(iTrack,4)= trackITS->GetSigned1Pt();;
2764  vecMomR0(iTrack,5)= trackITS->GetTgl();;
2765  }
2766  }
2767  //
2768  // 1.) Find closest matching tracks, between the ITS standalone track
2769  // and the all other tracks
2770  // a.) caltegory - All
2771  // b.) category - without ITS
2772  //
2773  //
2774  Int_t ntracksPropagated=0;
2775  AliExternalTrackParam extTrackDummy;
2776  AliESDtrack esdTrackDummy;
2777  AliExternalTrackParam itsAtTPC;
2778  AliExternalTrackParam itsAtITSTPC;
2779  for (Int_t iTrack0=0; iTrack0<ntracks; iTrack0++){
2780  AliESDtrack *track0 = esdEvent->GetTrack(iTrack0);
2781  if(!track0) continue;
2782  if (track0->IsOn(AliVTrack::kTPCin)) continue;
2783  AliESDfriendTrack* friendTrack0=esdFriend->GetTrack(iTrack0);
2784  if (!friendTrack0) continue;
2785  //if (!track0->IsOn(AliVTrack::kITSpureSA)) continue;
2786  //if (!friendTrack0->GetITSOut()) continue; // is there flag for ITS standalone?
2787  ntracksPropagated++;
2788  //
2789  // 2.) find clostest TPCtrack
2790  // a.) all tracks
2791  Double_t minChi2All=10000000;
2792  Double_t minChi2TPC=10000000;
2793  Double_t minChi2TPCITS=10000000;
2794  Int_t indexAll=-1;
2795  Int_t indexTPC=-1;
2796  Int_t indexTPCITS=-1;
2797  Int_t ncandidates0=0; // n candidates - rough cut
2798  Int_t ncandidates1=0; // n candidates - rough + chi2 cut
2799  itsAtTPC=*(friendTrack0->GetITSOut());
2800  itsAtITSTPC=*(friendTrack0->GetITSOut());
2801  for (Int_t iTrack1=0; iTrack1<ntracks; iTrack1++){
2802  AliESDtrack *track1 = esdEvent->GetTrack(iTrack1);
2803  if(!track1) continue;
2804  if (!track1->IsOn(AliVTrack::kTPCin)) continue;
2805  // fast checks
2806  //
2807  if (TMath::Abs(vecPosR1(iTrack1,2)-vecPosR0(iTrack0,2))>dFastZCut) continue;
2808  if (TMath::Abs(vecPosR1(iTrack1,3)-vecPosR0(iTrack0,3))>dFastPosPhiCut) continue;
2809  if (TMath::Abs(vecMomR1(iTrack1,3)-vecMomR0(iTrack0,3))>dFastPhiCut) continue;
2810  if (TMath::Abs(vecMomR1(iTrack1,5)-vecMomR0(iTrack0,5))>dFastThetaCut) continue;
2811  if (TMath::Abs(vecMomR1(iTrack1,4)-vecMomR0(iTrack0,4))>dFastPtCut) continue;
2812  ncandidates0++;
2813  //
2814  const AliExternalTrackParam * param1= track1->GetInnerParam();
2815  if (!friendTrack0->GetITSOut()) continue;
2816  AliExternalTrackParam outerITS = *(friendTrack0->GetITSOut());
2817  if (!outerITS.Rotate(param1->GetAlpha())) continue;
2818  if (!outerITS.PropagateTo(param1->GetX(),bz)) continue; // assume track close to the TPC inner wall
2819  Double_t chi2 = outerITS.GetPredictedChi2(param1);
2820  if (chi2>chi2Cut) continue;
2821  ncandidates1++;
2822  if (chi2<minChi2All){
2823  minChi2All=chi2;
2824  indexAll=iTrack1;
2825  }
2826  if (chi2<minChi2TPC && track1->IsOn(AliVTrack::kITSin)==0){
2827  minChi2TPC=chi2;
2828  indexTPC=iTrack1;
2829  itsAtTPC=outerITS;
2830  }
2831  if (chi2<minChi2TPCITS && track1->IsOn(AliVTrack::kITSin)){
2832  minChi2TPCITS=chi2;
2833  indexTPCITS=iTrack1;
2834  itsAtITSTPC=outerITS;
2835  }
2836  }
2837  //
2838  AliESDtrack * trackAll= (indexAll>=0)? esdEvent->GetTrack(indexAll):&esdTrackDummy;
2839  AliESDtrack * trackTPC= (indexTPC>=0)? esdEvent->GetTrack(indexTPC):&esdTrackDummy;
2840  AliESDtrack * trackTPCITS= (indexTPCITS>=0)? esdEvent->GetTrack(indexTPCITS):&esdTrackDummy;
2841  (*fTreeSRedirector)<<"itsTPC"<<
2842  "indexAll="<<indexAll<< // index of closest track (chi2)
2843  "indexTPC="<<indexTPC<< // index of closest TPCalone tracks
2844  "indexTPCITS="<<indexTPCITS<< // index of closest cobined tracks
2845  "ncandidates0="<<ncandidates0<< // number of candidates
2846  "ncandidates1="<<ncandidates1<<
2847  //
2848  "chi2All="<<minChi2All<< // chi2 of closest tracks
2849  "chi2TPC="<<minChi2TPC<<
2850  "chi2TPCITS="<<minChi2TPCITS<<
2851  //
2852  "track0.="<<track0<< // ITS standalone tracks
2853  "trackAll.="<<trackAll<< // Closets other track
2854  "trackTPC.="<<trackTPC<< // Closest TPC only track
2855  "trackTPCITS.="<<trackTPCITS<< // closest combined track
2856  //
2857  "itsAtTPC.="<<&itsAtTPC<< // ITS track parameters at the TPC alone track frame
2858  "itsAtITSTPC.="<<&itsAtITSTPC<< // ITS track parameters at the TPC combeined track frame
2859  "\n";
2860  }
2861 }
2862 
2863 void AliAnalysisTaskFilteredTree::ProcessTrackMatch(AliESDEvent *const /*esdEvent*/, AliESDfriend *const /*esdFriend*/){
2864 /*
2865  Use TPC standalone, ITS standalone and combined tracks to categorize/resp. recover track information.
2866 
2867  Track categories:
2868  -TPC+ITS
2869  -TPC only
2870  -ITS only
2871  Topology categories:
2872  -Nice isolated tracks with full information
2873  -Overlapped tracks - Refit and sign them
2874  -Multiple found (check overlap factor) - Merge and sign
2875  -Charge particle (kink) decays - Change of direction - Sign them)
2876  Info:
2877  -Array of indexes of closest tracks in each individual category
2878  -Chi2 distance to the closest tracks at reference radius of TPCin
2879  -Overlap factors - fraction of shared clusters or missing region
2880  -Chi2 distance at DCA
2881  Information matrix:
2882  -matrix closest tracks from each categories
2883  -characterization - chi2, index,chi2, overlap ratio
2884 
2885  Decissions:
2886  0.) Kink decay or catastophic multiple scaterring
2887  (combining all track categories)
2888  - small chi2 at the DCA
2889  - significantly large deflection angle
2890  - Combinatorial algorithm - to decrease CPU time restriction of investigation to tracks with small DCA at refernce radius
2891 
2892  1.) if (TPConly && !(TPC+ITS) && ITSonly match ) {
2893  if (!kink) TPCOnly.addoptITS
2894  if (kink) TPConly sign
2895  }
2896 
2897  2.) Overlap tracks - Refit with doUnfold
2898 */
2899 }
2900 
2901 Int_t AliAnalysisTaskFilteredTree::GetNearestTrack(const AliExternalTrackParam * trackMatch, Int_t indexSkip, AliESDEvent*event, Int_t trackType, Int_t paramType, AliExternalTrackParam & paramNearest){
2902  //
2903  // Find track with closest chi2 distance (assume all track ae propagated to the DCA)
2904  // trackType = 0 - find closets ITS standalone
2905  // 1 - find closest track with TPC
2906  // 2 - closest track with ITS and TPC
2907  // paramType = 0 - global track
2908  // 1 - track at inner wall of TPC
2909  //
2910  //
2911  if (trackMatch==NULL){
2912  ::Error("AliAnalysisTaskFilteredTree::GetNearestTrack","invalid track pointer");
2913  return -1;
2914  }
2915  Int_t ntracks=event->GetNumberOfTracks();
2916  const Double_t ktglCut=0.1;
2917  const Double_t kqptCut=0.4;
2918  const Double_t kAlphaCut=0.2;
2919  //
2920  Double_t chi2Min=100000;
2921  Int_t indexMin=-1;
2922  for (Int_t itrack=0; itrack<ntracks; itrack++){
2923  if (itrack==indexSkip) continue;
2924  AliESDtrack *ptrack=event->GetTrack(itrack);
2925  if (ptrack==NULL) continue;
2926  if (trackType==0 && (ptrack->IsOn(0x1)==kFALSE || ptrack->IsOn(0x10)==kTRUE)) continue; // looks for track without TPC information
2927  if (trackType==1 && (ptrack->IsOn(0x10)==kFALSE)) continue; // looks for tracks with TPC information
2928  if (trackType==2 && (ptrack->IsOn(0x1)==kFALSE || ptrack->IsOn(0x10)==kFALSE)) continue; // looks for tracks with TPC+ITS information
2929 
2930  if (ptrack->GetKinkIndex(0)<0) continue; // skip kink daughters
2931  const AliExternalTrackParam * track=0; //
2932  if (paramType==0) track=ptrack; // Global track
2933  if (paramType==1) track=ptrack->GetInnerParam(); // TPC only track at inner wall of TPC
2934  if (track==NULL) {
2935  continue;
2936  }
2937  // first rough cuts
2938  // fP3 cut
2939  if (TMath::Abs((track->GetTgl()-trackMatch->GetTgl()))>ktglCut) continue;
2940  // fP4 cut
2941  if (TMath::Abs((track->GetSigned1Pt()-trackMatch->GetSigned1Pt()))>kqptCut) continue;
2942  // fAlpha cut
2943  //Double_t alphaDist=TMath::Abs((track->GetAlpha()-trackMatch->GetAlpha()));
2944  Double_t alphaDist=TMath::Abs(TMath::ATan2(track->Py(),track->Px())-TMath::ATan2(trackMatch->Py(),trackMatch->Py()));
2945  if (alphaDist>TMath::Pi()) alphaDist-=TMath::TwoPi();
2946  if (alphaDist>kAlphaCut) continue;
2947  // calculate and extract track with smallest chi2 distance
2948  AliExternalTrackParam param(*track);
2949  if (param.Rotate(trackMatch->GetAlpha())==kFALSE) continue;
2950  if (param.PropagateTo(trackMatch->GetX(),trackMatch->GetBz())==kFALSE) continue;
2951  Double_t chi2=trackMatch->GetPredictedChi2(&param);
2952  if (chi2<chi2Min){
2953  indexMin=itrack;
2954  chi2Min=chi2;
2955  paramNearest=param;
2956  }
2957  }
2958  return indexMin;
2959 
2960 }
2961 
2962 
2963 void AliAnalysisTaskFilteredTree::SetDefaultAliasesV0(TTree *tree){
2964  //
2965  // SetAliases and Metadata for the V0 trees
2966  //
2967  TDatabasePDG *pdg = TDatabasePDG::Instance();
2968  Double_t massLambda = pdg->GetParticle("Lambda0")->Mass();
2969  Double_t massK0 = pdg->GetParticle("K0")->Mass();
2970  Double_t massPion = pdg->GetParticle("pi+")->Mass();
2971  Double_t massProton = pdg->GetParticle("proton")->Mass();
2972  const Double_t livetimeK0=2.684341668932; // livetime in cm (surpisely missing info in PDG - see root forum)
2973  const Double_t livetimeLambda=7.8875395; // livetime in cm (missing info in PDG - see root forum)
2974  //
2975  tree->SetAlias("massPion",Form("(%f+0)",massPion));
2976  tree->SetAlias("massProton",Form("(%f+0)",massProton));
2977  tree->SetAlias("massK0",Form("(%f+0)",massK0));
2978  tree->SetAlias("massLambda",Form("(%f+0)",massLambda));
2979  //
2980  tree->SetAlias("livetimeK0",TString::Format("%f",livetimeK0));
2981  tree->SetAlias("livetimeLambda",TString::Format("%f",livetimeLambda));
2982  tree->SetAlias("livetimeLikeK0",TString::Format("exp(-v0.fRr/(sqrt((v0.P()/%f)^2+1)*%f))",massK0, livetimeK0)); //
2983  tree->SetAlias("livetimeLikeLambda",TString::Format("exp(-v0.fRr/(sqrt((v0.P()/%f)^2+1)*%f))",massLambda,livetimeLambda)); //
2984  tree->SetAlias("livetimeLikeGamma","v0.fRr/80"); //
2985  tree->SetAlias("livetimeLikeBkg","v0.fRr/80"); //
2986  // delta of mass
2987  tree->SetAlias("K0Delta","(v0.GetEffMass(2,2)-massK0)");
2988  tree->SetAlias("LDelta","(v0.GetEffMass(4,2)-massLambda)");
2989  tree->SetAlias("ALDelta","(v0.GetEffMass(2,4)-massLambda)");
2990  tree->SetAlias("EDelta","(v0.GetEffMass(0,0))");
2991  // pull of the mass
2992  tree->SetAlias("K0Pull","(v0.GetEffMass(2,2)-massK0)/v0.GetKFInfo(2,2,1)");
2993  tree->SetAlias("LPull","(v0.GetEffMass(4,2)-massLambda)/v0.GetKFInfo(4,2,1)");
2994  tree->SetAlias("ALPull","(v0.GetEffMass(2,4)-massLambda)/v0.GetKFInfo(2,4,1)");
2995  tree->SetAlias("EPull","EDelta/v0.GetKFInfo(0,0,1)");
2996  // effective pull of the mass - (empirical values from fits)
2997  tree->SetAlias("K0PullEff","K0Delta/sqrt((3.63321e-03)**2+(5.68795e-04*v0.Pt())**2)");
2998  tree->SetAlias("LPullEff","LDelta/sqrt((1.5e-03)**2+(1.8e-04*v0.Pt())**2)");
2999  tree->SetAlias("ALPullEff","ALDelta/sqrt((1.5e-03)**2+(1.8e-04*v0.Pt())**2)");
3000  tree->SetAlias("EPullEff","v0.GetEffMass(0,0)/sqrt((5e-03)**2+(1.e-04*v0.Pt())**2)");
3001  //
3002  tree->SetAlias("dEdx0DProton","AliMathBase::BetheBlochAleph(track0.fIp.P()/massProton)");
3003  tree->SetAlias("dEdx1DProton","AliMathBase::BetheBlochAleph(track1.fIp.P()/massProton)");
3004  tree->SetAlias("dEdx0DPion","AliMathBase::BetheBlochAleph(track0.fIp.P()/massPion)");
3005  tree->SetAlias("dEdx1DPion","AliMathBase::BetheBlochAleph(track1.fIp.P()/massPion)");
3006  // V0 - cuts - for PID
3007  tree->SetAlias("cutDist","sqrt((track0.fIp.fP[0]-track1.fIp.fP[0])**2+(track0.fIp.fP[1]-track1.fIp.fP[1])**2)>3");
3008  tree->SetAlias("cutLong","track0.GetTPCClusterInfo(3,1,0)-5*abs(track0.fP[4])>130&&track1.GetTPCClusterInfo(3,1,0)>130-5*abs(track0.fP[4])");
3009  tree->SetAlias("cutPID","track0.fTPCsignal>0&&track1.fTPCsignal>0");
3010  tree->SetAlias("cutResol","sqrt(track0.fC[14]/track0.fP[4])<0.15&&sqrt(track1.fC[14]/track1.fP[4])<0.15");
3011  tree->SetAlias("cutV0","cutPID&&cutLong&&cutResol");
3012  //
3013  tree->SetAlias("K0PullBkg","min(min(abs(LPull),abs(ALPull)),abs(EPull))+0");
3014  tree->SetAlias("LambdaPullBkg","min(min(abs(K0Pull),abs(ALPull)),abs(EPull)+0)");
3015  tree->SetAlias("ALambdaPullBkg","min(min(abs(K0Pull),abs(LPull)),abs(EPull)+0)");
3016  tree->SetAlias("EPullBkg","min(min(abs(K0Pull),abs(LPull)),abs(ALPull)+0)");
3017  //
3018  tree->SetAlias("K0Selected", "abs(K0Pull)<3. &&abs(K0PullEff)<3. && abs(LPull)>3 && abs(ALPull)>3 &&v0.PtArmV0()>0.11");
3019  tree->SetAlias("LambdaSelected", "abs(LPull)<3. &&abs(LPullEff)<3. && abs(K0Pull)>3 && abs(EPull)>3 && abs(EDelta)>0.05");
3020  tree->SetAlias("ALambdaSelected", "abs(ALPull)<3. &&abs(ALPullEff)<3 && abs(K0Pull)>3 && abs(EPull)>3 &&abs(EDelta)>0.05");
3021  tree->SetAlias("GammaSelected", "abs(EPull)<3 && abs(K0Pull)>3 && abs(LPull)>3 && abs(ALPull)>3");
3022 
3023  tree->SetAlias("K0Like0","exp(-K0Pull^2)*livetimeLikeK0");
3024  tree->SetAlias("LLike0","exp(-LPull^2)*livetimeLikeLambda");
3025  tree->SetAlias("ALLike0","exp(-ALPull^2)*livetimeLikeLambda");
3026  tree->SetAlias("ELike0","exp(-abs(EPull)*0.2)*livetimeLikeGamma");
3027  tree->SetAlias("BkgLike","0.000005*ntracks"); // backround coeefecint to be fitted - depends on other cuts
3028  //
3029  tree->SetAlias("V0Like","exp(-acos(v0.fPointAngle)*v0.fRr/0.36)*exp(-sqrt(kf.GetChi2())/0.5)");
3030  tree->SetAlias("ELike","(V0Like*ELike0)/(V0Like*(K0Like0+LLike0+ALLike0+ELike0)+BkgLike)");
3031  tree->SetAlias("K0Like","K0Like0/(K0Like0+LLike0+ALLike0+ELike0+BkgLike)");
3032  tree->SetAlias("LLike","LLike0/(K0Like0+LLike0+ALLike0+ELike0+BkgLike)");
3033  tree->SetAlias("ALLike","ALLike0/(K0Like0+LLike0+ALLike0+ELike0+BkgLike)");
3034  //
3035  tree->SetAlias("K0PIDPull","(abs(track0.fTPCsignal/dEdx0DPion-50)+abs(track1.fTPCsignal/dEdx1DPion-50))/5.");
3036  tree->SetAlias("mpt","1/v0.Pt()"); //
3037  tree->SetAlias("tglV0","v0.Pz()/v0.Pt()"); //
3038  tree->SetAlias("alphaV0","atan2(v0.Py(),v0.Px()+0)");
3039  tree->SetAlias("dalphaV0","alphaV0-((int(36+9*(alphaV0/pi))-36)*pi/9.)");
3040 
3041 }
3042 
3043 void AliAnalysisTaskFilteredTree::SetDefaultAliasesHighPt(TTree *tree){
3044  //
3045  // set shortcut aliases for some variables
3046  //
3047  tree->SetAlias("phiInner","atan2(esdTrack.fIp.Py(),esdTrack.fIp.Px()+0)");
3048  tree->SetAlias("secInner","9*(atan2(esdTrack.fIp.Py(),esdTrack.fIp.Px()+0)/pi)+18*(esdTrack.fIp.Py()<0)");
3049  tree->SetAlias("tgl","esdTrack.fP[3]");
3050  tree->SetAlias("alphaV","esdTrack.fAlpha");
3051  tree->SetAlias("qPt","esdTrack.fP[4]");
3052  tree->SetAlias("dalphaQ","sign(esdTrack.fP[4])*(esdTrack.fIp.fP[0]/esdTrack.fIp.fX)");
3053  TStatToolkit::AddMetadata(tree,"phiInner.Title","#phi_{TPCin}");
3054  TStatToolkit::AddMetadata(tree,"secInner.Title","sector_{TPCin}");
3055  TStatToolkit::AddMetadata(tree,"tgl.Title","#it{p_{z}}/#it{p}_{T}");
3056  TStatToolkit::AddMetadata(tree,"alphaV.Title","#phi_{vertex}");
3057  TStatToolkit::AddMetadata(tree,"qPt.Title","q/#it{p}_{T}");
3058  TStatToolkit::AddMetadata(tree,"phiInner.AxisTitle","#it{#phi}_{TPCin}");
3059  TStatToolkit::AddMetadata(tree,"secInner.AxisTitle","sector_{TPCin}");
3060  TStatToolkit::AddMetadata(tree,"tgl.AxisTitle","#it{p}_{z}/#it{p}_{t}");
3061  TStatToolkit::AddMetadata(tree,"alphaV.AxisTitle","#it{#phi}_{vertex}");
3062  TStatToolkit::AddMetadata(tree,"qPt.AxisTitle","q/#it{p}_{T} (1/GeV)");
3063 
3064  //
3065  tree->SetAlias("normChi2ITS","sqrt(esdTrack.fITSchi2/esdTrack.fITSncls)");
3066  tree->SetAlias("normChi2TPC","esdTrack.fTPCchi2/esdTrack.fTPCncls");
3067  tree->SetAlias("normChi2TRD","esdTrack.fTRDchi2/esdTrack.fTRDncls");
3068  tree->SetAlias("normDCAR","esdTrack.fdTPC/sqrt(1+esdTrack.fP[4]**2)");
3069  tree->SetAlias("normDCAZ","esdTrack.fzTPC/sqrt(1+esdTrack.fP[4]**2)");
3070  TStatToolkit::AddMetadata(tree,"normChi2ITS.Title","#sqrt{#chi2_{ITS}/N_{clITS}}");
3071  TStatToolkit::AddMetadata(tree,"normChi2TPC.Title","#chi2_{TPC}/N_{clTPC}");
3072  TStatToolkit::AddMetadata(tree,"normChi2ITS.AxisTitle","#sqrt{#chi2_{ITS}/N_{clITS}}");
3073  TStatToolkit::AddMetadata(tree,"normChi2TPC.AxisTitle","#chi2_{TPC}/N_{clTPC}");
3074  //
3075  tree->SetAlias("TPCASide","esdTrack.fIp.fP[1]>0");
3076  tree->SetAlias("TPCCSide","esdTrack.fIp.fP[1]<0");
3077  tree->SetAlias("TPCCross","esdTrack.fIp.fP[1]*esdTrack.fIp.fP[3]<0");
3078  tree->SetAlias("ITSOn","((esdTrack.fFlags&0x1)>0)");
3079  tree->SetAlias("TPCOn","((esdTrack.fFlags&0x10)>0)");
3080  tree->SetAlias("ITSRefit","((esdTrack.fFlags&0x4)>0)");
3081  tree->SetAlias("TPCRefit","((esdTrack.fFlags&0x40)>0)");
3082  tree->SetAlias("TOFOn","((esdTrack.fFlags&0x2000)>0)");
3083  tree->SetAlias("TRDOn","((esdTrack.fFlags&0x400)>0)");
3084  tree->SetAlias("ITSOn0","esdTrack.fITSncls>4&&esdTrack.HasPointOnITSLayer(0)&&esdTrack.HasPointOnITSLayer(1)");
3085  tree->SetAlias("ITSOn01","esdTrack.fITSncls>3&&(esdTrack.HasPointOnITSLayer(0)||esdTrack.HasPointOnITSLayer(1))");
3086  tree->SetAlias("nclCut","(esdTrack.GetTPCClusterInfo(3,1)+esdTrack.fTRDncls)>140-5*(abs(esdTrack.fP[4]))");
3087  tree->SetAlias("IsPrim4","sqrt((esdTrack.fD**2)/esdTrack.fCdd+(esdTrack.fZ**2)/esdTrack.fCzz)<4");
3088  tree->SetAlias("IsPrim4TPC","sqrt((esdTrack.fdTPC**2)/esdTrack.fCddTPC+(esdTrack.fzTPC**2)/esdTrack.fCzzTPC)<4");
3089 
3090 
3091  const char * chName[5]={"#it{r#phi}","#it{z}","sin(#phi)","tan(#it{#theta})", "q/#it{p}_{t}"};
3092  const char * chUnit[5]={"cm","cm","","", "(1/GeV)"};
3093  const char * refBranch=(tree->GetBranch("extInnerParamV."))? "extInnerParamV":"esdTrack.fTPCInner";
3094 
3095  for (Int_t iPar=0; iPar<5; iPar++){
3096  tree->SetAlias(TString::Format("covarP%dITS",iPar).Data(),TString::Format("sqrt(esdTrack.fC[%d]+0)",AliExternalTrackParam::GetIndex(iPar,iPar)).Data());
3097  tree->SetAlias(TString::Format("covarP%d",iPar).Data(),TString::Format("sqrt(%s.fC[%d]+0)",refBranch,AliExternalTrackParam::GetIndex(iPar,iPar)).Data());
3098  tree->SetAlias(TString::Format("covarPC%d",iPar).Data(),TString::Format("sqrt(extInnerParamC.fC[%d]+0)",AliExternalTrackParam::GetIndex(iPar,iPar)).Data());
3099  tree->SetAlias(TString::Format("covarP%dNorm",iPar).Data(),TString::Format("sqrt(%s.fC[%d]+0)/sqrt(1+(1*esdTrack.fP[4])**2)/sqrt(1+(1*esdTrack.fP[4])**2)",refBranch,AliExternalTrackParam::GetIndex(iPar,iPar)).Data());
3100  tree->SetAlias(TString::Format("covarPC%dNorm",iPar).Data(),TString::Format("sqrt(extInnerParamC.fC[%d]+0)/sqrt(1+(5*esdTrack.fP[4])**2)",AliExternalTrackParam::GetIndex(iPar,iPar)).Data());
3101 
3102  tree->SetAlias(TString::Format("deltaP%d",iPar).Data(),TString::Format("(%s.fP[%d]-esdTrack.fCp.fP[%d])",refBranch,iPar,iPar).Data());
3103  tree->SetAlias(TString::Format("deltaPC%d",iPar).Data(),TString::Format("(extInnerParamC.fP[%d]-esdTrack.fCp.fP[%d])",iPar,iPar).Data());
3104  // rough normalization of the residual sigma ~ sqrt(1+*k/pt)^2)
3105  // histogram pt normalized deltas enable us to use wider bins in q/pt and also less bins in deltas
3106  tree->SetAlias(TString::Format("deltaP%dNorm",iPar).Data(),TString::Format("(%s.fP[%d]-esdTrack.fCp.fP[%d])/sqrt(1+(1*esdTrack.fP[4])**2)",refBranch,iPar,iPar).Data());
3107  tree->SetAlias(TString::Format("deltaPC%dNorm",iPar).Data(),TString::Format("(extInnerParamC.fP[%d]-esdTrack.fCp.fP[%d])/sqrt(1.+(5.*esdTrack.fP[4])**2)",iPar,iPar).Data());
3108  //
3109  tree->SetAlias(TString::Format("pullP%d",iPar).Data(),
3110  TString::Format("(%s.fP[%d]-esdTrack.fCp.fP[%d])/sqrt(%s.fC[%d]+esdTrack.fCp.fC[%d])",refBranch,
3111  iPar,iPar,refBranch,AliExternalTrackParam::GetIndex(iPar,iPar),AliExternalTrackParam::GetIndex(iPar,iPar)).Data());
3112  tree->SetAlias(TString::Format("pullPC%d",iPar).Data(),
3113  TString::Format("(extInnerParamC.fP[%d]-esdTrack.fCp.fP[%d])/sqrt(extInnerParamC.fC[%d]+esdTrack.fCp.fC[%d])",
3114  iPar,iPar,AliExternalTrackParam::GetIndex(iPar,iPar),AliExternalTrackParam::GetIndex(iPar,iPar)).Data());
3115  TStatToolkit::AddMetadata(tree,TString::Format("deltaP%d.AxisTitle",iPar).Data(),TString::Format("%s (%s)",chName[iPar], chUnit[iPar]).Data());
3116  TStatToolkit::AddMetadata(tree,TString::Format("deltaPC%d.AxisTitle",iPar).Data(),TString::Format("%s (%s)",chName[iPar], chUnit[iPar]).Data());
3117  TStatToolkit::AddMetadata(tree,TString::Format("pullP%d.AxisTitle",iPar).Data(),TString::Format("pull %s (unit)",chName[iPar]).Data());
3118  TStatToolkit::AddMetadata(tree,TString::Format("pullPC%d.AxisTitle",iPar).Data(),TString::Format("pull %s (unit)",chName[iPar]).Data());
3119  //
3120  TStatToolkit::AddMetadata(tree,TString::Format("deltaP%d.Title",iPar).Data(),TString::Format("%s",chName[iPar]).Data());
3121  TStatToolkit::AddMetadata(tree,TString::Format("deltaPC%d.Title",iPar).Data(),TString::Format("%s",chName[iPar]).Data());
3122  TStatToolkit::AddMetadata(tree,TString::Format("pullP%d.Title",iPar).Data(),TString::Format("pull %s",chName[iPar]).Data());
3123  TStatToolkit::AddMetadata(tree,TString::Format("pullPC%d.Title",iPar).Data(),TString::Format("pull %s",chName[iPar]).Data());
3124  }
3125  TStatToolkit::AddMetadata(tree, "mult.Title","N_{prim}");
3126  TStatToolkit::AddMetadata(tree, "mult.AxisTitle","N_{prim}");
3127  TStatToolkit::AddMetadata(tree, "ntracks.Title","N_{tr}");
3128  TStatToolkit::AddMetadata(tree, "ntracks.AxisTitle","N_{tr} (prim+sec+pile-up)");
3129 }
3130 
3139 Int_t AliAnalysisTaskFilteredTree::GetMCTrackDiff(const TParticle &particle, const AliExternalTrackParam &param, TClonesArray &trackRefArray, TVectorF &mcDiff){
3140  Double_t xyz[3]={0};
3141  param.GetXYZ(xyz);
3142  // 1.) Find nearest reference
3143  Int_t nTrackRefs = trackRefArray.GetEntriesFast();
3144  Int_t nTPCRef=0;
3145  AliTrackReference *refNearest=0;
3146  TVectorF fdist(nTrackRefs);
3147  for (Int_t itrR = 0; itrR < nTrackRefs; ++itrR) {
3148  AliTrackReference* ref = static_cast<AliTrackReference*>(trackRefArray.UncheckedAt(itrR));
3149  Double_t lDist=(ref->X()-xyz[0])*(ref->X()-xyz[0])+(ref->Y()-xyz[1])*(ref->Y()-xyz[1])+(ref->Z()-xyz[2])*(ref->Z()-xyz[2]);
3150  fdist[itrR]=lDist;
3151  }
3152  Double_t dist=250*250;
3153  Int_t index0=0,index1=0;
3154  for (Int_t itrR = 1; itrR < nTrackRefs-1; ++itrR){
3155  if (fdist[itrR]<dist){
3156  dist=fdist[itrR];
3157  if (fdist[itrR-1]<fdist[itrR+1]){
3158  index0=itrR-1;
3159  index1=itrR;
3160  }else{
3161  index0=itrR;
3162  index1=itrR+1;
3163  }
3164  }
3165  refNearest=static_cast<AliTrackReference*>(trackRefArray.UncheckedAt(index0));
3166  }
3167  // 2.) Get MC snapshot interpolation
3168  if (index0<0) return -1;
3169  AliTrackReference *ref0=static_cast<AliTrackReference*>(trackRefArray.UncheckedAt(index0));
3170  AliTrackReference *ref1=static_cast<AliTrackReference*>(trackRefArray.UncheckedAt(index0));
3171  Double_t xyz0[3]={ref0->X(), ref0->Y(), ref0->Z()};
3172  Double_t pxyz0[3]={ref0->Px(), ref0->Py(), ref0->Pz()};
3173  Double_t xyz1[3]={ref1->X(), ref1->Y(), ref1->Z()};
3174  Double_t pxyz1[3]={ref1->Px(), ref1->Py(), ref1->Pz()};
3175  AliExternalTrackParam param0(xyz0,pxyz0,(Double_t*)param.GetCovariance(),param.GetSign()); // TODO - fix constructor in AliExternalTrackParam using constant pointers
3176  AliExternalTrackParam param1(xyz1,pxyz1,(Double_t*)param.GetCovariance(),param.GetSign());
3177  Int_t errorCode=0;
3178  errorCode+=(param0.Rotate(param.GetAlpha())==kFALSE);
3179  errorCode+=(param1.Rotate(param.GetAlpha())==kFALSE);
3180  if (errorCode>0){
3181  return 100+errorCode;
3182  }
3183  errorCode+=(AliTrackerBase::PropagateTrackToBxByBz(&param0,param.GetX(),particle.GetMass(),5.,kFALSE,0.9999)==kFALSE);
3184  errorCode+=(AliTrackerBase::PropagateTrackToBxByBz(&param1,param.GetX(),particle.GetMass(),5.,kFALSE,0.9999)==kFALSE);
3185  if (errorCode>0){
3186  return 1000 +errorCode;
3187  }
3188  AliTrackerBase::UpdateTrack(param0,param1);
3189  // 3.) fill diff
3190  for (Int_t iPar=0; iPar<5; iPar++){
3191  mcDiff[iPar]=param.GetParameter()[iPar]-param0.GetParameter()[iPar];
3192  }
3193  return 0;
3194 }
3195 
3208 Int_t AliAnalysisTaskFilteredTree::GetMCInfoTrack(Int_t label, std::map<std::string,float> &trackInfoF, std::map<std::string,TObject*> &trackInfoO){
3209  // 0.) define some constants
3210  const Double_t kTPCOutR=245; // used in loop counters
3211  AliStack * stack = fMC->Stack();
3212  Int_t mcStackSize=stack->GetNtrack();
3213  if (label>mcStackSize){
3214  return 1;
3215  }
3216  // 1.) particle information
3217  TParticle *particle=NULL;
3218  TClonesArray *trackRefs=0;
3219  Int_t status = fMC->GetParticleAndTR(label, particle, trackRefs);
3220  trackInfoO["p"]=particle; // particle information
3221  if (particle==NULL || particle->GetPDG() ==NULL || particle->GetPDG()->Charge()!=0.) {
3222  return 2;
3223  }
3224  // 2.) particle trajectory information
3225  Int_t nTrackRef = trackRefs->GetEntries();
3226  trackInfoF["nRef"]=nTrackRef; // number of references
3227  if (nTrackRef==0){
3228  return 4;
3229  }
3230  TVectorF refCounter(21); // counter of references
3231  TVectorF detLength(21); // track length defined as distance between first and the last track reference in detector
3232  Double_t maxRadius=0;
3233  AliTrackReference*detRef[21]={NULL};
3234  Int_t loopCounter=0; // turning point counter to check loopers
3235  AliTrackReference *lastLoopRef = NULL;
3236  for (Int_t iRef = 0; iRef < nTrackRef; iRef++) {
3237  AliTrackReference *ref = (AliTrackReference *) trackRefs->At(iRef);
3238  if (ref->Label() != label) continue;
3239  Int_t detID=ref->DetectorId();
3240  if (detID < 0) {
3241  trackInfoO["refDecay"] = ref;
3242  break;
3243  }
3244  if (ref->R()>maxRadius) maxRadius=ref->R();
3245  refCounter[detID]++;
3246  if (lastLoopRef!=NULL && ref->R()<kTPCOutR){ //loop counter
3247  Double_t dir0=ref->Px()*ref->X()+ref->Py()*ref->Y();
3248  Double_t dir1=lastLoopRef->Px()*lastLoopRef->X()+lastLoopRef->Py()*lastLoopRef->Y();
3249  if (dir0*dir1<0) {
3250  loopCounter++;
3251  lastLoopRef=ref;
3252  }
3253  }
3254  if (lastLoopRef==NULL) lastLoopRef=ref;
3255  if (refCounter[detID] >0){
3256  detLength[detID]=ref->GetLength()-detRef[detID]->GetLength();
3257  }else{
3258  detRef[detID]=ref;
3259  }
3260  }
3261  trackInfoO["refCounter"]=new TVectorF(refCounter);
3262  trackInfoO["detLength"]=new TVectorF(detLength);
3263  trackInfoF["loopCounter"]=loopCounter;
3264  trackInfoF["maxRadius"]=maxRadius;
3265  // 3.) Assign - reconstruction information
3266  // In case particle reconstructed more than once - use the best
3267  // Distance definition ???
3268  Int_t ntracks=fESD->GetNumberOfTracks();
3269  Int_t nRecITS=0, nRecTPC=0, nRecTRD=0, nRecTOF=0;
3270  AliESDtrack * esdTrack=NULL;
3271  AliESDtrack * itsTrack=0;
3272  Int_t detRecLength=0, trackIndex=-1;
3273  for (Int_t iTrack=0;iTrack<ntracks;iTrack++) {
3274  AliESDtrack *track = fESD->GetTrack(iTrack);
3275  if (track == NULL) continue;
3276  Int_t recoLabel = TMath::Abs(track->GetLabel());
3277  if (recoLabel != label) continue;
3278  // find longest combined track - only "non fake legs" counted
3279  Int_t detRecLength0 = 0;
3280  if (TMath::Abs(track->GetITSLabel()) == label) detRecLength += track->IsOn(AliESDtrack::kITSin) +
3281  track->IsOn(AliESDtrack::kITSrefit);
3282  if (TMath::Abs(track->GetTPCLabel()) == label) detRecLength += track->IsOn(AliESDtrack::kTPCin) +
3283  track->IsOn(AliESDtrack::kTPCrefit);
3284  if (TMath::Abs(track->GetTRDLabel()) == label) detRecLength += track->IsOn(AliESDtrack::kTRDout) +
3285  track->IsOn(AliESDtrack::kTRDrefit);
3286  if (track->IsOn(AliESDtrack::kITSin)) nRecITS++;
3287  if (track->IsOn(AliESDtrack::kTPCin)) nRecTPC++;
3288  if (track->IsOn(AliESDtrack::kTRDout)) nRecTRD++;
3289  // in case the same "detector length" use most precise angular pz/pt determination
3290  // TODO - pz/pt chosen because valid also for secondary particles, better to use combined chi2, More complex - closest reference to be used in that case
3291  if (detRecLength == detRecLength0) {
3292  Double_t tglParticle = particle->Pz() / particle->Pt();
3293  Double_t deltaTgl0 = esdTrack->Pz() / esdTrack->Pt() - tglParticle;
3294  Double_t deltaTgl1 = track->Pz() / track->Pt() - tglParticle;
3295  if (TMath::Abs(deltaTgl1) < TMath::Abs(deltaTgl0)) {
3296  esdTrack = track;
3297  trackIndex = iTrack;
3298  }
3299  }
3300  if (detRecLength < detRecLength0) {
3301  detRecLength = detRecLength0;
3302  esdTrack = track;
3303  trackIndex = iTrack;
3304  }
3305  if (!track->IsOn(AliESDtrack::kTPCin) && track->IsOn(AliESDtrack::kITSin)) {
3306  itsTrack = track;
3307  }
3308  }
3309  trackInfoO["esdTrack"]=esdTrack;
3310  trackInfoO["itsTrack"]=itsTrack;
3311  //4.) diff between MC and real data at reference planes
3312  TVectorF mcDiff(5);
3313  if (esdTrack!=NULL){
3314  if (GetMCTrackDiff(*particle,*(esdTrack), *trackRefs, mcDiff)==0){
3315  trackInfoO["diffesdTrack"]=new TVectorF(mcDiff);
3316  };
3317  if (esdTrack->GetTPCInnerParam()){
3318  if (GetMCTrackDiff(*particle,*(esdTrack->GetTPCInnerParam()), *trackRefs, mcDiff)==0){
3319  trackInfoO["diffTPCInnerParam"]=new TVectorF(mcDiff);
3320  };
3321  }
3322  if (esdTrack->GetInnerParam()){
3323  if (GetMCTrackDiff(*particle,*(esdTrack->GetInnerParam()), *trackRefs, mcDiff)==0){
3324  trackInfoO["diffInnerParam"]=new TVectorF(mcDiff);
3325  };
3326  }
3327  if (esdTrack->GetOuterParam()){
3328  if (GetMCTrackDiff(*particle,*(esdTrack->GetOuterParam()), *trackRefs, mcDiff)==0){
3329  trackInfoO["diffOuterParam"]=new TVectorF(mcDiff);
3330  };
3331  }
3332  if (esdTrack->GetOuterHmpParam()){
3333  if (GetMCTrackDiff(*particle,*(esdTrack->GetOuterHmpParam()), *trackRefs, mcDiff)==0){
3334  trackInfoO["diffOuterHmpParam"]=new TVectorF(mcDiff);
3335  };
3336  }
3337  }
3338  return 0;
3339 }
3340 
3341 Int_t AliAnalysisTaskFilteredTree::GetMCInfoKink(Int_t label, std::map<std::string,float> &kinkInfoF, std::map<std::string,TObject*> &kinkInfoO){
3342 
3343 }
3344 
3348 void AliAnalysisTaskFilteredTree::ProcessMC(){
3349  //
3350  AliStack * stack = fMC->Stack();
3351  if (!stack) return;
3352  Int_t mcStackSize=stack->GetNtrack();
3353  std::map<std::string,float> trackInfoF;
3354  std::map<std::string,TObject*> trackInfoO;
3355  static Int_t downscaleCounter=0;
3356  for (Int_t iMc = 0; iMc < mcStackSize; ++iMc) {
3357  TParticle *particle = stack->Particle(iMc);
3358  if (!particle) continue;
3359  // apply downscaling function
3360  Double_t scalempt= TMath::Min(particle->Pt(),10.);
3361  Double_t downscaleF = gRandom->Rndm();
3362  downscaleF *= fLowPtTrackDownscaligF;
3363  if (downscaleCounter>0 && TMath::Exp(2*scalempt)<downscaleF) continue;
3364  Int_t result = GetMCInfoTrack(iMc, trackInfoF,trackInfoO);
3365 
3366  }
3367 }
charge
Int_t charge
Definition: MuonTrackingEfficiency.C:54
AliAnalysisTaskFilteredTree::ProcessV0
void ProcessV0(AliESDEvent *const esdEvent=0, AliMCEvent *const mcEvent=0, AliESDfriend *const esdFriend=0)
Definition: AliAnalysisTaskFilteredTree.cxx:1844
pdg
Int_t pdg
Definition: charmCutsOptimization.C:34
filename
const char * filename
Definition: TestFCM.C:1
Double_t
double Double_t
Definition: External.C:58
AliAnalysisTaskFilteredTree::SetDefaultAliasesHighPt
static void SetDefaultAliasesHighPt(TTree *treeV0)
Definition: AliAnalysisTaskFilteredTree.cxx:3043
AliAnalysisTaskFilteredTree::ProcessITSTPCmatchOut
void ProcessITSTPCmatchOut(AliESDEvent *const esdEvent=0, AliESDfriend *const esdFriend=0)
Definition: AliAnalysisTaskFilteredTree.cxx:2701
AliAnalysisTaskFilteredTree::GetMCInfoTrack
Int_t GetMCInfoTrack(Int_t label, std::map< std::string, float > &trackInfoF, std::map< std::string, TObject * > &trackInfoO)
Definition: AliAnalysisTaskFilteredTree.cxx:3208
AliAnalysisTaskFilteredTree::IsFromConversion
Bool_t IsFromConversion(Int_t label, AliStack *const stack)
Definition: AliAnalysisTaskFilteredTree.cxx:2455
AliFilteredTreeEventCuts.h
AliAnalysisTaskFilteredTree::ConstrainTPCInner
Bool_t ConstrainTPCInner(AliExternalTrackParam *const tpcInnerC, const AliESDVertex *vtx, Double_t b[3])
Definition: AliAnalysisTaskFilteredTree.cxx:2381
AliAnalysisTaskFilteredTree::IsFromMaterial
Bool_t IsFromMaterial(Int_t label, AliStack *const stack)
Definition: AliAnalysisTaskFilteredTree.cxx:2487
AliAnalysisTaskFilteredTree::ConstrainTrackInner
Bool_t ConstrainTrackInner(AliExternalTrackParam *const trackInnerC, const AliESDVertex *vtx, Double_t mass, Double_t b[3])
Definition: AliAnalysisTaskFilteredTree.cxx:2407
mass
Double_t mass
Definition: charmCutsOptimization.C:35
gSystem
TSystem * gSystem
bz
Double_t bz
Definition: AliMCTreeTools.cxx:77
kRadius
const Float_t kRadius[3]
Definition: AliAnalysisTaskFullppJet.cxx:56
c
TCanvas * c
Definition: TestFitELoss.C:172
AliAnalysisTaskFilteredTree::ProcessAll
void ProcessAll(AliESDEvent *const esdEvent=0, AliMCEvent *const mcEvent=0, AliESDfriend *const esdFriend=0)
Definition: AliAnalysisTaskFilteredTree.cxx:787
AliAnalysisTaskFilteredTree::GetMCTrackDiff
static Int_t GetMCTrackDiff(const TParticle &particle, const AliExternalTrackParam &param, TClonesArray &trackRefArray, TVectorF &mcDiff)
Definition: AliAnalysisTaskFilteredTree.cxx:3139
stack
AliStack * stack
Definition: AliMCTreeTools.cxx:72
gRandom
TRandom * gRandom
AliAnalysisTaskFilteredTree::IsV0Downscaled
Bool_t IsV0Downscaled(AliESDv0 *const v0)
Definition: AliAnalysisTaskFilteredTree.cxx:2343
kProton
Definition: FastCentEstimators.C:25
AliAnalysisTaskFilteredTree::UserExec
virtual void UserExec(Option_t *option)
Definition: AliAnalysisTaskFilteredTree.cxx:305
Int_t
int Int_t
Definition: External.C:63
UInt_t
unsigned int UInt_t
Definition: External.C:33
Float_t
float Float_t
Definition: External.C:68
TH3D
Definition: External.C:252
AliAnalysisTaskFilteredTree::ProcessTrackMatch
void ProcessTrackMatch(AliESDEvent *const esdEvent=0, AliESDfriend *const esdFriend=0)
Definition: AliAnalysisTaskFilteredTree.cxx:2863
AliAnalysisTaskFilteredTree.h
GetIndex
Int_t GetIndex(TObjArray *triggersB, Int_t trigNr, Int_t centNr)
Definition: PlotMuonQA.C:2046
AliAnalysisTaskFilteredTree::ProcessCosmics
void ProcessCosmics(AliESDEvent *const esdEvent=0, AliESDfriend *esdFriend=0)
Definition: AliAnalysisTaskFilteredTree.cxx:384
AliAnalysisTaskFilteredTree::ProcessdEdx
void ProcessdEdx(AliESDEvent *const esdEvent=0, AliMCEvent *const mcEvent=0, AliESDfriend *const esdFriend=0)
Definition: AliAnalysisTaskFilteredTree.cxx:2079
kEMCAL
Definition: ComparePi0CalibResults.C:44
Data
Bool_t Data(TH1F *h, Double_t *rangefit, Bool_t writefit, Double_t &sgn, Double_t &errsgn, Double_t &bkg, Double_t &errbkg, Double_t &sgnf, Double_t &errsgnf, Double_t &sigmafit, Int_t &status)
Definition: charmCutsOptimization.C:371
AliAnalysisTaskFilteredTree::GetNearestTrack
Int_t GetNearestTrack(const AliExternalTrackParam *trackMatch, Int_t indexSkip, AliESDEvent *event, Int_t trackType, Int_t paramType, AliExternalTrackParam &paramNearest)
Definition: AliAnalysisTaskFilteredTree.cxx:2901
AliAnalysisTaskFilteredTree::GetMCInfoKink
Int_t GetMCInfoKink(Int_t label, std::map< std::string, float > &kinkInfoF, std::map< std::string, TObject * > &kinkInfoO)
Definition: AliAnalysisTaskFilteredTree.cxx:3341
AliFilteredTreeAcceptanceCuts.h
AliAnalysisTaskFilteredTree::FillHistograms
void FillHistograms(AliESDtrack *const ptrack, AliExternalTrackParam *const ptpcInnerC, Double_t centralityF, Double_t chi2TPCInnerC)
Definition: AliAnalysisTaskFilteredTree.cxx:2629
AliFilteredTreeEventCuts::AcceptMCEvent
Bool_t AcceptMCEvent(AliMCEvent *mcEvent=0)
Definition: AliFilteredTreeEventCuts.cxx:127
AliAnalysisTaskFilteredTree::GetMother
TParticle * GetMother(TParticle *const particle, AliStack *const stack)
Definition: AliAnalysisTaskFilteredTree.cxx:2440
AliAnalysisTaskFilteredTree::GetKFParticle
Int_t GetKFParticle(AliESDv0 *const v0, AliESDEvent *const event, AliKFParticle &kfparticle)
Definition: AliAnalysisTaskFilteredTree.cxx:2220
Option_t
const char Option_t
Definition: External.C:48
AliAnalysisTaskFilteredTree::ProcessLaser
void ProcessLaser(AliESDEvent *const esdEvent=0, AliMCEvent *const mcEvent=0, AliESDfriend *const esdFriend=0)
Definition: AliAnalysisTaskFilteredTree.cxx:738
Process
void Process(Int_t *pflag[23040][7], TH1 *inhisto, Double_t Nsigma=4., Int_t dnbins=200, Double_t dmaxval=-1., Int_t compteur=1)
Definition: BadChannelAnalysis.C:344
AliAnalysisTaskFilteredTree::GetMCTrueTrackMult
static Int_t GetMCTrueTrackMult(AliMCEvent *const mcEvent, AliFilteredTreeEventCuts *const evtCuts, AliFilteredTreeAcceptanceCuts *const accCuts)
Definition: AliAnalysisTaskFilteredTree.cxx:2579
AliAnalysisTaskFilteredTree::Process
void Process(AliESDEvent *const esdEvent=0, AliMCEvent *const mcEvent=0, AliESDfriend *const esdFriend=0)
Definition: AliAnalysisTaskFilteredTree.cxx:537
Bool_t
bool Bool_t
Definition: External.C:53
AliAnalysisTaskFilteredTree::IsFromStrangeness
Bool_t IsFromStrangeness(Int_t label, AliStack *const stack)
Definition: AliAnalysisTaskFilteredTree.cxx:2517
AliAnalysisTaskFilteredTree::SetDefaultAliasesV0
static void SetDefaultAliasesV0(TTree *treeV0)
Definition: AliAnalysisTaskFilteredTree.cxx:2963
AliFilteredTreeEventCuts::AcceptEvent
Bool_t AcceptEvent(AliESDEvent *event=0, AliMCEvent *mcEvent=0, const AliESDVertex *vtx=0)
Definition: AliFilteredTreeEventCuts.cxx:89
trackRefs
TClonesArray * trackRefs
Definition: AliMCTreeTools.cxx:65
OpenFile
TList * OpenFile(const char *fname)
Definition: DrawAnaELoss.C:65
AliAnalysisTaskFilteredTree::IsHighDeDxParticle
Bool_t IsHighDeDxParticle(AliESDtrack *const track)
Definition: AliAnalysisTaskFilteredTree.cxx:1825
AliAnalysisTaskFilteredTree::ProcessMCEff
void ProcessMCEff(AliESDEvent *const esdEvent=0, AliMCEvent *const mcEvent=0, AliESDfriend *const esdFriend=0)
Definition: AliAnalysisTaskFilteredTree.cxx:1548