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