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