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