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