]> 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
Coverity fix
[u/mrichter/AliRoot.git] / PWGPP / AliAnalysisTaskFilteredTree.cxx
CommitLineData
1059f583 1/**************************************************************************\r
2* Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *\r
3* *\r
4* Author: The ALICE Off-line Project. *\r
5* Contributors are mentioned in the code where appropriate. *\r
6* *\r
7* Permission to use, copy, modify and distribute this software and its *\r
8* documentation strictly for non-commercial purposes is hereby granted *\r
9* without fee, provided that the above copyright notice appears in all *\r
10* copies and that both the copyright notice and this permission notice *\r
11* appear in the supporting documentation. The authors make no claims *\r
12* about the suitability of this software for any purpose. It is *\r
13* provided "as is" without express or implied warranty. *\r
14**************************************************************************/\r
15\r
16#include "iostream"\r
17\r
18#include <TPDGCode.h>\r
19#include <TDatabasePDG.h>\r
20\r
21#include "TChain.h"\r
22#include "TTreeStream.h"\r
23#include "TTree.h"\r
24#include "TH1F.h"\r
25#include "TCanvas.h"\r
26#include "TList.h"\r
27#include "TObjArray.h"\r
28#include "TFile.h"\r
29#include "TMatrixD.h"\r
30#include "TRandom3.h"\r
31\r
32#include "AliHeader.h" \r
33#include "AliGenEventHeader.h" \r
34#include "AliInputEventHandler.h" \r
35#include "AliStack.h" \r
36#include "AliTrackReference.h" \r
37\r
38#include "AliPhysicsSelection.h"\r
39#include "AliAnalysisTask.h"\r
40#include "AliAnalysisManager.h"\r
41#include "AliESDEvent.h"\r
42#include "AliESDfriend.h"\r
43#include "AliMCEvent.h"\r
44#include "AliESDInputHandler.h"\r
45#include "AliESDVertex.h"\r
46#include "AliTracker.h"\r
47#include "AliGeomManager.h"\r
48\r
49#include "AliCentrality.h"\r
50#include "AliESDVZERO.h"\r
51#include "AliMultiplicity.h"\r
52\r
53#include "AliESDtrackCuts.h"\r
54#include "AliMCEventHandler.h"\r
55#include "AliFilteredTreeEventCuts.h"\r
56#include "AliFilteredTreeAcceptanceCuts.h"\r
57\r
58#include "AliAnalysisTaskFilteredTree.h"\r
59#include "AliKFParticle.h"\r
60#include "AliESDv0.h"\r
61\r
62using namespace std;\r
63\r
64ClassImp(AliAnalysisTaskFilteredTree)\r
65\r
66//_____________________________________________________________________________\r
67AliAnalysisTaskFilteredTree::AliAnalysisTaskFilteredTree(const char *name) \r
68 : AliAnalysisTaskSE(name)\r
69 , fESD(0)\r
70 , fMC(0)\r
71 , fESDfriend(0)\r
72 , fOutput(0)\r
73 , fPitList(0)\r
74 , fUseMCInfo(kFALSE)\r
75 , fUseESDfriends(kFALSE)\r
76 , fFilteredTreeEventCuts(0)\r
77 , fFilteredTreeAcceptanceCuts(0)\r
78 , fFilteredTreeRecAcceptanceCuts(0)\r
79 , fEsdTrackCuts(0)\r
80 , fTrigger(AliTriggerAnalysis::kMB1) \r
81 , fAnalysisMode(kTPCAnalysisMode) \r
82 , fTreeSRedirector(0)\r
83 , fCentralityEstimator(0)\r
84 , fLowPtTrackDownscaligF(0)\r
85 , fLowPtV0DownscaligF(0)\r
86 , fProcessAll(kFALSE)\r
87 , fProcessCosmics(kFALSE)\r
88 , fHighPtTree(0)\r
89 , fV0Tree(0)\r
90 , fdEdxTree(0)\r
91 , fLaserTree(0)\r
92 , fMCEffTree(0)\r
93 , fCosmicPairsTree(0)\r
94{\r
95 // Constructor\r
96\r
97 // Define input and output slots here\r
98 DefineOutput(1, TTree::Class());\r
99 DefineOutput(2, TTree::Class());\r
100 DefineOutput(3, TTree::Class());\r
101 DefineOutput(4, TTree::Class());\r
102 DefineOutput(5, TTree::Class());\r
103 DefineOutput(6, TTree::Class());\r
104}\r
105\r
106//_____________________________________________________________________________\r
107AliAnalysisTaskFilteredTree::~AliAnalysisTaskFilteredTree()\r
108{\r
a1419baa 109 //the output trees not to be deleted in case of proof analysis\r
110 Bool_t deleteTrees=kTRUE;\r
e4aa7a8e 111 if ((AliAnalysisManager::GetAnalysisManager()))\r
112 {\r
113 if (AliAnalysisManager::GetAnalysisManager()->GetAnalysisType() == \r
114 AliAnalysisManager::kProofAnalysis)\r
a1419baa 115 deleteTrees=kFALSE;\r
e4aa7a8e 116 }\r
a1419baa 117 if (deleteTrees) delete fTreeSRedirector;\r
e4aa7a8e 118\r
a1419baa 119 delete fFilteredTreeEventCuts;\r
120 delete fFilteredTreeAcceptanceCuts;\r
121 delete fFilteredTreeRecAcceptanceCuts;\r
122 delete fEsdTrackCuts;\r
1059f583 123}\r
124\r
125//____________________________________________________________________________\r
126Bool_t AliAnalysisTaskFilteredTree::Notify()\r
127{\r
128 static Int_t count = 0;\r
129 count++;\r
130 TTree *chain = (TChain*)GetInputData(0);\r
131 if(chain)\r
132 {\r
133 Printf("Processing %d. file: %s", count, chain->GetCurrentFile()->GetName());\r
134 }\r
135\r
136return kTRUE;\r
137}\r
138\r
139//_____________________________________________________________________________\r
140void AliAnalysisTaskFilteredTree::UserCreateOutputObjects()\r
141{\r
142 // Create histograms\r
143 // Called once\r
1059f583 144\r
145 //\r
a1419baa 146 //get the output file to make sure the trees will be associated to it\r
147 OpenFile(1);\r
148 fTreeSRedirector = new TTreeSRedirector();\r
1059f583 149 \r
1059f583 150 //\r
151 // Create trees\r
a1419baa 152 fV0Tree = ((*fTreeSRedirector)<<"V0s").GetTree();\r
153 fHighPtTree = ((*fTreeSRedirector)<<"highPt").GetTree();\r
154 fdEdxTree = ((*fTreeSRedirector)<<"dEdx").GetTree();\r
155 fLaserTree = ((*fTreeSRedirector)<<"Laser").GetTree();\r
156 fMCEffTree = ((*fTreeSRedirector)<<"MCEffTree").GetTree();\r
157 fCosmicPairsTree = ((*fTreeSRedirector)<<"CosmicPairs").GetTree();\r
158\r
159 PostData(1,fV0Tree);\r
160 PostData(2,fHighPtTree);\r
1059f583 161 PostData(3,fdEdxTree);\r
162 PostData(4,fLaserTree);\r
163 PostData(5,fMCEffTree);\r
164 PostData(6,fCosmicPairsTree);\r
1059f583 165}\r
166\r
167//_____________________________________________________________________________\r
168void AliAnalysisTaskFilteredTree::UserExec(Option_t *) \r
169{\r
170 //\r
171 // Called for each event\r
172 //\r
173\r
174 // ESD event\r
175 fESD = (AliESDEvent*) (InputEvent());\r
176 if (!fESD) {\r
177 Printf("ERROR: ESD event not available");\r
178 return;\r
179 }\r
180\r
181 // MC event\r
182 if(fUseMCInfo) {\r
183 fMC = MCEvent();\r
184 if (!fMC) {\r
185 Printf("ERROR: MC event not available");\r
186 return;\r
187 }\r
188 }\r
189\r
190 if(fUseESDfriends) {\r
191 fESDfriend = static_cast<AliESDfriend*>(fESD->FindListObject("AliESDfriend"));\r
192 if(!fESDfriend) {\r
193 Printf("ERROR: ESD friends not available");\r
194 }\r
195 }\r
196\r
197 //\r
198 if(fProcessAll) { \r
199 ProcessAll(fESD,fMC,fESDfriend); // all track stages and MC\r
200 }\r
201 else {\r
202 Process(fESD,fMC,fESDfriend); // only global and TPC tracks\r
203 }\r
204\r
205 //\r
206 ProcessV0(fESD,fMC,fESDfriend);\r
207 ProcessLaser(fESD,fMC,fESDfriend);\r
208 ProcessdEdx(fESD,fMC,fESDfriend);\r
209\r
210 if (fProcessCosmics) { ProcessCosmics(fESD); }\r
211 if(IsUseMCInfo()) { ProcessMCEff(fESD,fMC,fESDfriend); }\r
212}\r
213\r
214//_____________________________________________________________________________\r
215void AliAnalysisTaskFilteredTree::ProcessCosmics(AliESDEvent *const event)\r
216{\r
217 //\r
218 // Select real events with high-pT tracks \r
219 //\r
220 if(!event) {\r
221 AliDebug(AliLog::kError, "event not available");\r
222 return;\r
223 }\r
224\r
225 // \r
226 AliInputEventHandler* inputHandler = (AliInputEventHandler*) AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler();\r
227 if (!inputHandler)\r
228 {\r
229 Printf("ERROR: Could not receive input handler");\r
230 return;\r
231 }\r
232 \r
233 // get file name\r
234 TTree *chain = (TChain*)GetInputData(0);\r
235 if(!chain) { \r
236 Printf("ERROR: Could not receive input chain");\r
237 return;\r
238 }\r
239 TObjString fileName(chain->GetCurrentFile()->GetName());\r
240\r
241\r
242 // check for cosmic pairs\r
243 //\r
244 // find cosmic pairs trigger by random trigger\r
245 //\r
246 //\r
247 AliESDVertex *vertexSPD = (AliESDVertex *)event->GetPrimaryVertexSPD();\r
248 AliESDVertex *vertexTPC = (AliESDVertex *)event->GetPrimaryVertexTPC(); \r
249 const Double_t kMinPt=0.8;\r
250 const Double_t kMinPtMax=0.8;\r
251 const Double_t kMinNcl=50;\r
252 const Double_t kMaxDelta[5]={2,600,0.02,0.02,0.1};\r
253 Int_t ntracks=event->GetNumberOfTracks(); \r
254 // Float_t dcaTPC[2]={0,0};\r
255 // Float_t covTPC[3]={0,0,0};\r
256\r
257 UInt_t specie = event->GetEventSpecie(); // skip laser events\r
258 if (specie==AliRecoParam::kCalib) return;\r
259 \r
260\r
261\r
262 for (Int_t itrack0=0;itrack0<ntracks;itrack0++) {\r
263 AliESDtrack *track0 = event->GetTrack(itrack0);\r
264 if (!track0) continue;\r
265 if (!track0->IsOn(AliESDtrack::kTPCrefit)) continue;\r
266\r
267 if (TMath::Abs(AliTracker::GetBz())>1 && track0->Pt() < kMinPt) continue;\r
268 if (track0->Pt() < kMinPt) continue;\r
269 if (track0->GetTPCncls() < kMinNcl) continue;\r
270 if (TMath::Abs(track0->GetY())<kMaxDelta[0]) continue; \r
271 if (track0->GetKinkIndex(0)>0) continue;\r
272 const Double_t * par0=track0->GetParameter(); //track param at rhe DCA\r
273 //rm primaries\r
274 //\r
275 //track0->GetImpactParametersTPC(dcaTPC,covTPC);\r
276 //if (TMath::Abs(dcaTPC[0])<kMaxDelta[0]) continue;\r
277 //if (TMath::Abs(dcaTPC[1])<kMaxDelta[0]*2) continue;\r
278 // const AliExternalTrackParam * trackIn0 = track0->GetInnerParam();\r
279 for (Int_t itrack1=itrack0+1;itrack1<ntracks;itrack1++) {\r
280 AliESDtrack *track1 = event->GetTrack(itrack1);\r
281 if (!track1) continue; \r
282 if (!track1->IsOn(AliESDtrack::kTPCrefit)) continue;\r
283 if (track1->GetKinkIndex(0)>0) continue;\r
284 if ((TMath::Abs(AliTracker::GetBz())>1) && (track1->Pt() < kMinPt)) continue;\r
285 if (track1->Pt() < kMinPt) continue;\r
286 if (track1->GetTPCncls()<kMinNcl) continue;\r
287 if (TMath::Abs(AliTracker::GetBz())>1 && TMath::Max(track1->Pt(), track0->Pt())<kMinPtMax) continue;\r
288 if (TMath::Abs(track1->GetY())<kMaxDelta[0]) continue;\r
289 //track1->GetImpactParametersTPC(dcaTPC,covTPC);\r
290 // if (TMath::Abs(dcaTPC[0])<kMaxDelta[0]) continue;\r
291 //if (TMath::Abs(dcaTPC[1])<kMaxDelta[0]*2) continue;\r
292 //\r
293 const Double_t* par1=track1->GetParameter(); //track param at rhe DCA\r
294 //\r
295 Bool_t isPair=kTRUE;\r
296 for (Int_t ipar=0; ipar<5; ipar++){\r
297 if (ipar==4&&TMath::Abs(AliTracker::GetBz())<1) continue; // 1/pt not defined for B field off\r
298 if (TMath::Abs(TMath::Abs(par0[ipar])-TMath::Abs(par1[ipar]))>kMaxDelta[ipar]) isPair=kFALSE;\r
299 }\r
300 if (!isPair) continue;\r
301 if (TMath::Abs(TMath::Abs(track0->GetAlpha()-track1->GetAlpha())-TMath::Pi())>kMaxDelta[2]) isPair=kFALSE;\r
302 //delta with correct sign\r
303 /*\r
304 TCut cut0="abs(t1.fP[0]+t0.fP[0])<2"\r
305 TCut cut3="abs(t1.fP[3]+t0.fP[3])<0.02"\r
306 TCut cut4="abs(t1.fP[4]+t0.fP[4])<0.2"\r
307 */\r
308 if (TMath::Abs(par0[0]+par1[0])>kMaxDelta[0]) isPair=kFALSE; //delta y opposite sign\r
309 if (TMath::Abs(par0[3]+par1[3])>kMaxDelta[3]) isPair=kFALSE; //delta tgl opposite sign\r
310 if (TMath::Abs(AliTracker::GetBz())>1 && TMath::Abs(par0[4]+par1[4])>kMaxDelta[4]) isPair=kFALSE; //delta 1/pt opposite sign\r
311 if (!isPair) continue;\r
312 TString filename(AliAnalysisManager::GetAnalysisManager()->GetTree()->GetCurrentFile()->GetName());\r
313 Int_t eventNumber = event->GetEventNumberInFile(); \r
314 //Bool_t hasFriend = kFALSE;\r
315 //Bool_t hasITS=(track0->GetNcls(0)+track1->GetNcls(0)>4);\r
316 //printf("DUMPHPTCosmic:%s|%f|%d|%d|%d\n",filename.Data(),(TMath::Min(track0->Pt(),track1->Pt())), eventNumber,hasFriend,hasITS);\r
317 // const AliExternalTrackParam * trackIn1 = track1->GetInnerParam(); \r
318 //\r
319 // \r
320 Int_t ntracksSPD = vertexSPD->GetNContributors();\r
321 Int_t ntracksTPC = vertexTPC->GetNContributors(); \r
322 Int_t runNumber = event->GetRunNumber(); \r
323 Int_t timeStamp = event->GetTimeStamp();\r
324 ULong64_t triggerMask = event->GetTriggerMask();\r
325 Float_t magField = event->GetMagneticField();\r
326 TObjString triggerClass = event->GetFiredTriggerClasses().Data();\r
327 \r
328 //\r
329 // Dump to the tree \r
330 // vertex\r
331 // TPC-ITS tracks\r
332 //\r
a1419baa 333\r
334 //fCosmicPairsTree->Branch("fileName",&fileName,32000,0);\r
335 //fCosmicPairsTree->Branch("runNumber",&runNumber,"runNumber/I");\r
336 //fCosmicPairsTree->Branch("timeStamp",&timeStamp,"timeStamp/I");\r
337 //fCosmicPairsTree->Branch("eventNumber",&eventNumber,"eventNumber/I");\r
338 //fCosmicPairsTree->Branch("triggerMask",&triggerMask,32000,0);\r
339 //fCosmicPairsTree->Branch("triggerClass",&triggerClass,32000,0);\r
340 //fCosmicPairsTree->Branch("magField",&magField,"magField/F");\r
341 //fCosmicPairsTree->Branch("ntracksSPD",&ntracksSPD,"ntracksSPD/I");\r
342 //fCosmicPairsTree->Branch("ntracksTPC",&ntracksTPC,"ntracksTPC/I");\r
343 //fCosmicPairsTree->Branch("vertexSPD",vertexSPD,32000,0);\r
344 //fCosmicPairsTree->Branch("vertexTPC",vertexTPC,32000,0);\r
345 //fCosmicPairsTree->Branch("track0",track0,32000,0);\r
346 //fCosmicPairsTree->Branch("track1",track1,32000,0);\r
347 //\r
348 //fCosmicPairsTree->Fill();\r
349\r
1059f583 350 if(!fTreeSRedirector) return;\r
351 (*fTreeSRedirector)<<"CosmicPairs"<<\r
352 "fileName.="<<&fileName<< // file name\r
353 "runNumber="<<runNumber<< // run number \r
354 "evtTimeStamp="<<timeStamp<< // time stamp of event\r
355 "evtNumberInFile="<<eventNumber<< // event number \r
356 "trigger="<<triggerMask<< // trigger\r
04fe1e30 357 "triggerClass="<<&triggerClass<< // trigger\r
1059f583 358 "Bz="<<magField<< // magnetic field\r
359 //\r
360 "multSPD="<<ntracksSPD<<\r
361 "multTPC="<<ntracksTPC<<\r
362 "vertSPD.="<<vertexSPD<< //primary vertex -SPD\r
363 "vertTPC.="<<vertexTPC<< //primary vertex -TPC\r
364 "t0.="<<track0<< //track0\r
365 "t1.="<<track1<< //track1\r
366 "\n"; \r
367 }\r
368 }\r
369}\r
370\r
371\r
372//_____________________________________________________________________________\r
373void AliAnalysisTaskFilteredTree::Process(AliESDEvent *const esdEvent, AliMCEvent * const mcEvent, AliESDfriend *const /*esdFriend*/)\r
374{\r
375 //\r
376 // Select real events with high-pT tracks \r
377 //\r
378 if(!esdEvent) {\r
379 AliDebug(AliLog::kError, "esdEvent not available");\r
380 return;\r
381 }\r
382\r
383 // get selection cuts\r
384 AliFilteredTreeEventCuts *evtCuts = GetEventCuts(); \r
385 AliFilteredTreeAcceptanceCuts *accCuts = GetAcceptanceCuts(); \r
386 AliESDtrackCuts *esdTrackCuts = GetTrackCuts(); \r
387\r
388 if(!evtCuts || !accCuts || !esdTrackCuts) {\r
389 Printf("ERROR cuts not available");\r
390 return;\r
391 }\r
392\r
393 // trigger selection\r
394 Bool_t isEventTriggered = kTRUE;\r
395 AliPhysicsSelection *physicsSelection = NULL;\r
396 AliTriggerAnalysis* triggerAnalysis = NULL;\r
397\r
398 // \r
399 AliInputEventHandler* inputHandler = (AliInputEventHandler*) AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler();\r
400 if (!inputHandler)\r
401 {\r
402 Printf("ERROR: Could not receive input handler");\r
403 return;\r
404 }\r
405 \r
406 // get file name\r
407 TTree *chain = (TChain*)GetInputData(0);\r
408 if(!chain) { \r
409 Printf("ERROR: Could not receive input chain");\r
410 return;\r
411 }\r
412 TObjString fileName(chain->GetCurrentFile()->GetName());\r
413\r
414 // trigger\r
415 if(evtCuts->IsTriggerRequired()) \r
416 {\r
417 // always MB\r
418 isEventTriggered = inputHandler->IsEventSelected() & AliVEvent::kMB;\r
419\r
420 physicsSelection = static_cast<AliPhysicsSelection*> (inputHandler->GetEventSelection());\r
421 if(!physicsSelection) return;\r
422 //SetPhysicsTriggerSelection(physicsSelection);\r
423\r
424 if (isEventTriggered && (GetTrigger() == AliTriggerAnalysis::kV0AND)) {\r
425 // set trigger (V0AND)\r
426 triggerAnalysis = physicsSelection->GetTriggerAnalysis();\r
427 if(!triggerAnalysis) return;\r
428 isEventTriggered = triggerAnalysis->IsOfflineTriggerFired(esdEvent, GetTrigger());\r
429 }\r
430 }\r
431\r
432 // centrality determination\r
433 Float_t centralityF = -1;\r
434 AliCentrality *esdCentrality = esdEvent->GetCentrality();\r
435 centralityF = esdCentrality->GetCentralityPercentile(fCentralityEstimator.Data());\r
436\r
437 // use MC information\r
438 AliHeader* header = 0;\r
439 AliGenEventHeader* genHeader = 0;\r
440 AliStack* stack = 0;\r
441 TArrayF vtxMC(3);\r
442\r
443 Int_t multMCTrueTracks = 0;\r
444 if(IsUseMCInfo())\r
445 {\r
446 //\r
447 if(!mcEvent) {\r
448 AliDebug(AliLog::kError, "mcEvent not available");\r
449 return;\r
450 }\r
451 // get MC event header\r
452 header = mcEvent->Header();\r
453 if (!header) {\r
454 AliDebug(AliLog::kError, "Header not available");\r
455 return;\r
456 }\r
457 // MC particle stack\r
458 stack = mcEvent->Stack();\r
459 if (!stack) {\r
460 AliDebug(AliLog::kError, "Stack not available");\r
461 return;\r
462 }\r
463\r
464 // get MC vertex\r
465 genHeader = header->GenEventHeader();\r
466 if (!genHeader) {\r
467 AliDebug(AliLog::kError, "Could not retrieve genHeader from Header");\r
468 return;\r
469 }\r
470 genHeader->PrimaryVertex(vtxMC);\r
471\r
472 // multipliticy of all MC primary tracks\r
473 // in Zv, pt and eta ranges)\r
474 multMCTrueTracks = GetMCTrueTrackMult(mcEvent,evtCuts,accCuts);\r
475\r
476 } // end bUseMC\r
477\r
478 // get reconstructed vertex \r
479 //const AliESDVertex* vtxESD = 0; \r
480 AliESDVertex* vtxESD = 0; \r
481 if(GetAnalysisMode() == kTPCAnalysisMode) {\r
482 vtxESD = (AliESDVertex*)esdEvent->GetPrimaryVertexTPC();\r
483 }\r
484 else if(GetAnalysisMode() == kTPCITSAnalysisMode) {\r
485 vtxESD = (AliESDVertex*)esdEvent->GetPrimaryVertexTracks();\r
486 }\r
487 else {\r
488 return;\r
489 }\r
490\r
491 if(!vtxESD) return;\r
492\r
493 Bool_t isEventOK = evtCuts->AcceptEvent(esdEvent,mcEvent,vtxESD); \r
494 //printf("isEventOK %d, isEventTriggered %d, status %d, vz %f \n",isEventOK, isEventTriggered, vtxESD->GetStatus(), vtxESD->GetZv());\r
495 //printf("GetAnalysisMode() %d \n",GetAnalysisMode());\r
496\r
497\r
498 // check event cuts\r
499 if(isEventOK && isEventTriggered)\r
500 {\r
501\r
502 //\r
503 // get IR information\r
504 //\r
505 AliESDHeader *esdHeader = 0;\r
506 esdHeader = esdEvent->GetHeader();\r
507 if(!esdHeader) return;\r
508 //Int_t ir1 = esdHeader->GetTriggerIREntries(); //all ir-s\r
509 //Int_t ir2 = esdHeader->GetTriggerIREntries(-1,1); // int2 set, 180 ms time interval\r
510\r
511 // Use when Peter commit the changes in the header\r
512 Int_t ir1 = 0;\r
513 Int_t ir2 = 0;\r
514\r
515 //\r
516 Double_t vert[3] = {0}; \r
517 vert[0] = vtxESD->GetXv();\r
518 vert[1] = vtxESD->GetYv();\r
519 vert[2] = vtxESD->GetZv();\r
520 Int_t mult = vtxESD->GetNContributors();\r
a1419baa 521 Float_t bz = esdEvent->GetMagneticField();\r
522 Int_t runNumber = esdEvent->GetRunNumber();\r
523 Int_t evtTimeStamp = esdEvent->GetTimeStamp();\r
1059f583 524 Int_t evtNumberInFile = esdEvent->GetEventNumberInFile();\r
525\r
526 // high pT tracks\r
527 for (Int_t iTrack = 0; iTrack < esdEvent->GetNumberOfTracks(); iTrack++)\r
528 {\r
529 AliESDtrack *track = esdEvent->GetTrack(iTrack);\r
530 if(!track) continue;\r
531 if(track->Charge()==0) continue;\r
532 if(!esdTrackCuts->AcceptTrack(track)) continue;\r
533 if(!accCuts->AcceptTrack(track)) continue;\r
534 \r
535 // downscale low-pT tracks\r
536 Double_t scalempt= TMath::Min(track->Pt(),10.);\r
537 Double_t downscaleF = gRandom->Rndm();\r
538 downscaleF *= fLowPtTrackDownscaligF;\r
539 if(TMath::Exp(2*scalempt)<downscaleF) continue;\r
540 //printf("TMath::Exp(2*scalempt) %e, downscaleF %e \n",TMath::Exp(2*scalempt), downscaleF);\r
541\r
542 AliExternalTrackParam * tpcInner = (AliExternalTrackParam *)(track->GetTPCInnerParam());\r
543 if (!tpcInner) continue;\r
544 // transform to the track reference frame \r
545 Bool_t isOK = kFALSE;\r
546 isOK = tpcInner->Rotate(track->GetAlpha());\r
547 isOK = tpcInner->PropagateTo(track->GetX(),esdEvent->GetMagneticField());\r
548 if(!isOK) continue;\r
549\r
550 // Dump to the tree \r
551 // vertex\r
552 // TPC-ITS tracks\r
553 //\r
04fe1e30 554 TObjString triggerClass = esdEvent->GetFiredTriggerClasses().Data();\r
a1419baa 555\r
556 //fHighPtTree->Branch("fileName",&fileName,32000,0);\r
557 //fHighPtTree->Branch("runNumber",&runNumber,"runNumber/I");\r
558 //fHighPtTree->Branch("evtTimeStamp",&evtTimeStamp,"evtTimeStamp/I");\r
559 //fHighPtTree->Branch("evtNumberInFile",&evtNumberInFile,"evtNumberInFile/I");\r
560 //fHighPtTree->Branch("triggerClass",&triggerClass,32000,0);\r
561 //fHighPtTree->Branch("Bz",&bz,"Bz/F");\r
562 //fHighPtTree->Branch("vtxESD",vtxESD,32000,0);\r
563 //fHighPtTree->Branch("IRtot",&ir1,"IRtot/I");\r
564 //fHighPtTree->Branch("IRint2",&ir2,"IRint2/I");\r
565 //fHighPtTree->Branch("mult",&mult,"mult/I");\r
566 //fHighPtTree->Branch("esdTrack",track,32000,0);\r
567 //fHighPtTree->Branch("centralityF",&centralityF,"centralityF/F");\r
568\r
569 //fHighPtTree->Fill();\r
570\r
571 //Double_t vtxX=vtxESD->GetX();\r
572 //Double_t vtxY=vtxESD->GetY();\r
573 //Double_t vtxZ=vtxESD->GetZ();\r
1059f583 574 if(!fTreeSRedirector) return;\r
575 (*fTreeSRedirector)<<"highPt"<<\r
576 "fileName.="<<&fileName<<\r
577 "runNumber="<<runNumber<<\r
578 "evtTimeStamp="<<evtTimeStamp<<\r
579 "evtNumberInFile="<<evtNumberInFile<<\r
04fe1e30 580 "triggerClass="<<&triggerClass<< // trigger\r
1059f583 581 "Bz="<<bz<<\r
582 "vtxESD.="<<vtxESD<<\r
a1419baa 583 // "vtxESDx="<<vtxX<<\r
584 // "vtxESDy="<<vtxY<<\r
585 // "vtxESDz="<<vtxZ<<\r
1059f583 586 "IRtot="<<ir1<<\r
587 "IRint2="<<ir2<<\r
588 "mult="<<mult<<\r
589 "esdTrack.="<<track<<\r
590 "centralityF="<<centralityF<<\r
591 "\n";\r
592 }\r
593 }\r
594 \r
595}\r
596\r
597\r
598//_____________________________________________________________________________\r
599void AliAnalysisTaskFilteredTree::ProcessLaser(AliESDEvent *const esdEvent, AliMCEvent * const /*mcEvent*/, AliESDfriend *const /*esdFriend*/)\r
600{\r
601 //\r
602 // Process laser events\r
603 //\r
604 if(!esdEvent) {\r
605 AliDebug(AliLog::kError, "esdEvent not available");\r
606 return;\r
607 }\r
608\r
609 // get file name\r
610 TTree *chain = (TChain*)GetInputData(0);\r
611 if(!chain) { \r
612 Printf("ERROR: Could not receive input chain");\r
613 return;\r
614 }\r
615 TObjString fileName(chain->GetCurrentFile()->GetName());\r
616\r
617 // laser events \r
618 const AliESDHeader* esdHeader = esdEvent->GetHeader();\r
619 if(esdHeader && esdHeader->GetEventSpecie()==AliRecoParam::kCalib) \r
620 {\r
621 Int_t countLaserTracks = 0;\r
622 for (Int_t iTrack = 0; iTrack < esdEvent->GetNumberOfTracks(); iTrack++)\r
623 {\r
624 AliESDtrack *track = esdEvent->GetTrack(iTrack);\r
625 if(!track) continue;\r
626\r
627 if(track->GetTPCInnerParam()) countLaserTracks++;\r
628 }\r
629 \r
630 if(countLaserTracks > 100) \r
631 { \r
a1419baa 632 Int_t runNumber = esdEvent->GetRunNumber();\r
633 Int_t evtTimeStamp = esdEvent->GetTimeStamp();\r
1059f583 634 Int_t evtNumberInFile = esdEvent->GetEventNumberInFile();\r
a1419baa 635 Float_t bz = esdEvent->GetMagneticField();\r
04fe1e30 636 TObjString triggerClass = esdEvent->GetFiredTriggerClasses().Data();\r
1059f583 637\r
a1419baa 638 //fLaserTree->Branch("fileName",&fileName,32000,0);\r
639 //fLaserTree->Branch("runNumber",&runNumber,"runNumber/I");\r
640 //fLaserTree->Branch("evtTimeStamp",&evtTimeStamp,"evtTimeStamp/I");\r
641 //fLaserTree->Branch("evtNumberInFile",&evtNumberInFile,"evtNumberInFile/I");\r
642 //fLaserTree->Branch("triggerClass",&triggerClass,32000,0);\r
643 //fLaserTree->Branch("Bz",&bz,"Bz/F");\r
644 //fLaserTree->Branch("multTPCtracks",&countLaserTracks,"multTPCtracks/I");\r
645\r
646 //fLaserTree->Fill();\r
647\r
1059f583 648 if(!fTreeSRedirector) return;\r
649 (*fTreeSRedirector)<<"Laser"<<\r
650 "fileName.="<<&fileName<<\r
651 "runNumber="<<runNumber<<\r
652 "evtTimeStamp="<<evtTimeStamp<<\r
653 "evtNumberInFile="<<evtNumberInFile<<\r
04fe1e30 654 "triggerClass="<<&triggerClass<< // trigger\r
1059f583 655 "Bz="<<bz<<\r
656 "multTPCtracks="<<countLaserTracks<<\r
657 "\n";\r
658 }\r
659 }\r
660}\r
661\r
662//_____________________________________________________________________________\r
663void AliAnalysisTaskFilteredTree::ProcessAll(AliESDEvent *const esdEvent, AliMCEvent * const mcEvent, AliESDfriend *const esdFriend)\r
664{\r
665 //\r
666 // Select real events with high-pT tracks\r
667 // Calculate and stor in the output tree:\r
668 // TPC constrained tracks\r
669 // InnerParams constrained tracks\r
670 // TPC-ITS tracks\r
671 // ITSout-InnerParams tracks\r
672 // chi2 distance between TPC constrained and TPC-ITS tracks\r
673 // chi2 distance between TPC refitted constrained and TPC-ITS tracks\r
674 // chi2 distance between ITSout and InnerParams tracks\r
675 // MC information: \r
676 // track references at ITSin, TPCin; InnerParam at first TPC track reference, \r
677 // particle ID, mother ID, production mechanism ...\r
678 // \r
679\r
680 if(!esdEvent) {\r
681 AliDebug(AliLog::kError, "esdEvent not available");\r
682 return;\r
683 }\r
684\r
685 // get selection cuts\r
686 AliFilteredTreeEventCuts *evtCuts = GetEventCuts(); \r
687 AliFilteredTreeAcceptanceCuts *accCuts = GetAcceptanceCuts(); \r
688 AliESDtrackCuts *esdTrackCuts = GetTrackCuts(); \r
689\r
690 if(!evtCuts || !accCuts || !esdTrackCuts) {\r
691 AliDebug(AliLog::kError, "cuts not available");\r
692 return;\r
693 }\r
694\r
695 // trigger selection\r
696 Bool_t isEventTriggered = kTRUE;\r
697 AliPhysicsSelection *physicsSelection = NULL;\r
698 AliTriggerAnalysis* triggerAnalysis = NULL;\r
699\r
700 // \r
701 AliInputEventHandler* inputHandler = (AliInputEventHandler*) AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler();\r
702 if (!inputHandler)\r
703 {\r
704 Printf("ERROR: Could not receive input handler");\r
705 return;\r
706 }\r
707 \r
708 // get file name\r
709 TTree *chain = (TChain*)GetInputData(0);\r
710 if(!chain) { \r
711 Printf("ERROR: Could not receive input chain");\r
712 return;\r
713 }\r
714 TObjString fileName(chain->GetCurrentFile()->GetName());\r
715\r
716 // trigger\r
717 if(evtCuts->IsTriggerRequired()) \r
718 {\r
719 // always MB\r
720 isEventTriggered = inputHandler->IsEventSelected() & AliVEvent::kMB;\r
721\r
722 physicsSelection = static_cast<AliPhysicsSelection*> (inputHandler->GetEventSelection());\r
723 if(!physicsSelection) return;\r
724 //SetPhysicsTriggerSelection(physicsSelection);\r
725\r
726 if (isEventTriggered && (GetTrigger() == AliTriggerAnalysis::kV0AND)) {\r
727 // set trigger (V0AND)\r
728 triggerAnalysis = physicsSelection->GetTriggerAnalysis();\r
729 if(!triggerAnalysis) return;\r
730 isEventTriggered = triggerAnalysis->IsOfflineTriggerFired(esdEvent, GetTrigger());\r
731 }\r
732 }\r
733\r
734 // centrality determination\r
735 Float_t centralityF = -1;\r
736 AliCentrality *esdCentrality = esdEvent->GetCentrality();\r
737 centralityF = esdCentrality->GetCentralityPercentile(fCentralityEstimator.Data());\r
738\r
739 // use MC information\r
740 AliHeader* header = 0;\r
741 AliGenEventHeader* genHeader = 0;\r
742 AliStack* stack = 0;\r
743 TArrayF vtxMC(3);\r
744\r
745 Int_t multMCTrueTracks = 0;\r
746 if(IsUseMCInfo())\r
747 {\r
748 //\r
749 if(!mcEvent) {\r
750 AliDebug(AliLog::kError, "mcEvent not available");\r
751 return;\r
752 }\r
753 // get MC event header\r
754 header = mcEvent->Header();\r
755 if (!header) {\r
756 AliDebug(AliLog::kError, "Header not available");\r
757 return;\r
758 }\r
759 // MC particle stack\r
760 stack = mcEvent->Stack();\r
761 if (!stack) {\r
762 AliDebug(AliLog::kError, "Stack not available");\r
763 return;\r
764 }\r
765\r
766 // get MC vertex\r
767 genHeader = header->GenEventHeader();\r
768 if (!genHeader) {\r
769 AliDebug(AliLog::kError, "Could not retrieve genHeader from Header");\r
770 return;\r
771 }\r
772 genHeader->PrimaryVertex(vtxMC);\r
773\r
774 // multipliticy of all MC primary tracks\r
775 // in Zv, pt and eta ranges)\r
776 multMCTrueTracks = GetMCTrueTrackMult(mcEvent,evtCuts,accCuts);\r
777\r
778 } // end bUseMC\r
779\r
780 // get reconstructed vertex \r
781 //const AliESDVertex* vtxESD = 0; \r
782 AliESDVertex* vtxESD = 0; \r
783 if(GetAnalysisMode() == kTPCAnalysisMode) {\r
784 vtxESD = (AliESDVertex*)esdEvent->GetPrimaryVertexTPC();\r
785 }\r
786 else if(GetAnalysisMode() == kTPCITSAnalysisMode) {\r
787 vtxESD = (AliESDVertex*)esdEvent->GetPrimaryVertexTracks();\r
788 }\r
789 else {\r
790 return;\r
791 }\r
792\r
793 if(!vtxESD) return;\r
794\r
795 Bool_t isEventOK = evtCuts->AcceptEvent(esdEvent,mcEvent,vtxESD); \r
796 //printf("isEventOK %d, isEventTriggered %d \n",isEventOK, isEventTriggered);\r
797 //printf("GetAnalysisMode() %d \n",GetAnalysisMode());\r
798\r
799\r
800 // check event cuts\r
801 if(isEventOK && isEventTriggered)\r
802 {\r
803 //\r
804 // get IR information\r
805 //\r
806 AliESDHeader *esdHeader = 0;\r
807 esdHeader = esdEvent->GetHeader();\r
808 if(!esdHeader) return;\r
809 //Int_t ir1 = esdHeader->GetTriggerIREntries(); //all ir-s\r
810 //Int_t ir2 = esdHeader->GetTriggerIREntries(-1,1); // int2 set, 180 ms time interval\r
811 //\r
812 Int_t ir1 = 0;\r
813 Int_t ir2 = 0;\r
814\r
815 //\r
816 Double_t vert[3] = {0}; \r
817 vert[0] = vtxESD->GetXv();\r
818 vert[1] = vtxESD->GetYv();\r
819 vert[2] = vtxESD->GetZv();\r
820 Int_t mult = vtxESD->GetNContributors();\r
a1419baa 821 Float_t bz = esdEvent->GetMagneticField();\r
822 Int_t runNumber = esdEvent->GetRunNumber();\r
823 Int_t evtTimeStamp = esdEvent->GetTimeStamp();\r
1059f583 824 Int_t evtNumberInFile = esdEvent->GetEventNumberInFile();\r
825\r
826 // high pT tracks\r
827 for (Int_t iTrack = 0; iTrack < esdEvent->GetNumberOfTracks(); iTrack++)\r
828 {\r
829 AliESDtrack *track = esdEvent->GetTrack(iTrack);\r
830 if(!track) continue;\r
831 if(track->Charge()==0) continue;\r
832 if(!esdTrackCuts->AcceptTrack(track)) continue;\r
833 if(!accCuts->AcceptTrack(track)) continue;\r
834 \r
835 // downscale low-pT tracks\r
836 Double_t scalempt= TMath::Min(track->Pt(),10.);\r
837 Double_t downscaleF = gRandom->Rndm();\r
838 downscaleF *= fLowPtTrackDownscaligF;\r
839 if(TMath::Exp(2*scalempt)<downscaleF) continue;\r
840 //printf("TMath::Exp(2*scalempt) %e, downscaleF %e \n",TMath::Exp(2*scalempt), downscaleF);\r
841\r
842 // Dump to the tree \r
843 // vertex\r
844 // TPC constrained tracks\r
845 // InnerParams constrained tracks\r
846 // TPC-ITS tracks\r
847 // ITSout-InnerParams tracks\r
848 // chi2 distance between TPC constrained and TPC-ITS tracks\r
849 // chi2 distance between TPC refitted constrained and TPC-ITS tracks\r
850 // chi2 distance between ITSout and InnerParams tracks\r
851 // MC information\r
852 \r
853 Double_t x[3]; track->GetXYZ(x);\r
854 Double_t b[3]; AliTracker::GetBxByBz(x,b);\r
855\r
856 //\r
857 // Transform TPC inner params to track reference frame\r
858 //\r
859 Bool_t isOKtpcInner = kFALSE;\r
860 AliExternalTrackParam * tpcInner = (AliExternalTrackParam *)(track->GetTPCInnerParam());\r
861 if (tpcInner) {\r
862 // transform to the track reference frame \r
863 isOKtpcInner = tpcInner->Rotate(track->GetAlpha());\r
864 isOKtpcInner = tpcInner->PropagateTo(track->GetX(),esdEvent->GetMagneticField());\r
865 }\r
866\r
867 //\r
868 // Constrain TPC inner to vertex\r
869 // clone TPCinner has to be deleted\r
870 //\r
871 Bool_t isOKtpcInnerC = kFALSE;\r
872 AliExternalTrackParam * tpcInnerC = new AliExternalTrackParam(*(track->GetTPCInnerParam()));\r
873 if (tpcInnerC) {\r
874 isOKtpcInnerC = ConstrainTPCInner(tpcInnerC,vtxESD,b);\r
875 isOKtpcInnerC = tpcInnerC->Rotate(track->GetAlpha());\r
876 isOKtpcInnerC = tpcInnerC->PropagateTo(track->GetX(),esdEvent->GetMagneticField());\r
877 }\r
878\r
879 //\r
880 // Constrain TPC refitted tracks at inner TPC wall (InnerParams) to vertex \r
881 // Clone track InnerParams has to be deleted\r
882 //\r
883 Bool_t isOKtrackInnerC = kFALSE;\r
884 AliExternalTrackParam * trackInnerC = new AliExternalTrackParam(*(track->GetInnerParam()));\r
885 if (trackInnerC) {\r
886 isOKtrackInnerC = ConstrainTrackInner(trackInnerC,vtxESD,track->GetMass(),b);\r
887 isOKtrackInnerC = trackInnerC->Rotate(track->GetAlpha());\r
888 isOKtrackInnerC = trackInnerC->PropagateTo(track->GetX(),esdEvent->GetMagneticField());\r
889 } \r
890 \r
891 //\r
892 // calculate chi2 between vi and vj vectors\r
893 // with covi and covj covariance matrices\r
894 // chi2ij = (vi-vj)^(T)*(covi+covj)^(-1)*(vi-vj)\r
895 //\r
896 TMatrixD deltaT(5,1), deltaTtrackC(5,1);\r
897 TMatrixD delta(1,5), deltatrackC(1,5);\r
898 TMatrixD covarM(5,5), covarMtrackC(5,5);\r
899 TMatrixD chi2(1,1);\r
900 TMatrixD chi2trackC(1,1);\r
901\r
902 if(isOKtpcInnerC && isOKtrackInnerC) \r
903 {\r
904 for (Int_t ipar=0; ipar<5; ipar++) {\r
905 deltaT(ipar,0)=tpcInnerC->GetParameter()[ipar]-track->GetParameter()[ipar];\r
906 delta(0,ipar)=tpcInnerC->GetParameter()[ipar]-track->GetParameter()[ipar];\r
907\r
908 deltaTtrackC(ipar,0)=trackInnerC->GetParameter()[ipar]-track->GetParameter()[ipar];\r
909 deltatrackC(0,ipar)=trackInnerC->GetParameter()[ipar]-track->GetParameter()[ipar];\r
910\r
911 for (Int_t jpar=0; jpar<5; jpar++) {\r
912 Int_t index=track->GetIndex(ipar,jpar);\r
913 covarM(ipar,jpar)=track->GetCovariance()[index]+tpcInnerC->GetCovariance()[index];\r
914 covarMtrackC(ipar,jpar)=track->GetCovariance()[index]+trackInnerC->GetCovariance()[index];\r
915 }\r
916 }\r
917\r
918 // chi2 distance TPC constrained and TPC+ITS\r
919 TMatrixD covarMInv = covarM.Invert();\r
920 TMatrixD mat2 = covarMInv*deltaT;\r
921 chi2 = delta*mat2; \r
922 //chi2.Print();\r
923\r
924 // chi2 distance TPC refitted constrained and TPC+ITS\r
925 TMatrixD covarMInvtrackC = covarMtrackC.Invert();\r
926 TMatrixD mat2trackC = covarMInvtrackC*deltaTtrackC;\r
927 chi2trackC = deltatrackC*mat2trackC; \r
928 //chi2trackC.Print();\r
929 }\r
930\r
931\r
932 //\r
933 // Propagate ITSout to TPC inner wall \r
934 // and calculate chi2 distance to track (InnerParams)\r
935 //\r
936 const Double_t kTPCRadius=85; \r
937 const Double_t kStep=3; \r
938\r
939 // clone track InnerParams has to be deleted\r
940 Bool_t isOKtrackInnerC2 = kFALSE;\r
941 AliExternalTrackParam *trackInnerC2 = new AliExternalTrackParam(*(track->GetInnerParam()));\r
942 if (trackInnerC2) {\r
943 isOKtrackInnerC2 = AliTracker::PropagateTrackToBxByBz(trackInnerC2,kTPCRadius,track->GetMass(),kStep,kFALSE);\r
944 }\r
945\r
946 Bool_t isOKouterITSc = kFALSE;\r
947 AliExternalTrackParam *outerITSc = NULL;\r
948 TMatrixD chi2OuterITS(1,1);\r
949\r
950 if(esdFriend && esdFriend->TestSkipBit()==kFALSE) \r
951 {\r
952 // propagate ITSout to TPC inner wall\r
953 AliESDfriendTrack *friendTrack = esdFriend->GetTrack(iTrack);\r
954\r
955 if(friendTrack) \r
956 {\r
957 outerITSc = new AliExternalTrackParam(*friendTrack->GetITSOut());\r
958 if(outerITSc) \r
959 {\r
960 isOKouterITSc = AliTracker::PropagateTrackToBxByBz(outerITSc,kTPCRadius,track->GetMass(),kStep,kFALSE);\r
961 isOKouterITSc = outerITSc->Rotate(trackInnerC2->GetAlpha());\r
962 isOKouterITSc = outerITSc->PropagateTo(trackInnerC2->GetX(),esdEvent->GetMagneticField());\r
963\r
964 //\r
965 // calculate chi2 between outerITS and innerParams\r
966 // cov without z-coordinate at the moment\r
967 //\r
968 TMatrixD deltaTouterITS(4,1);\r
969 TMatrixD deltaouterITS(1,4);\r
970 TMatrixD covarMouterITS(4,4);\r
971\r
972 if(isOKtrackInnerC2 && isOKouterITSc) {\r
973 Int_t kipar = 0;\r
974 Int_t kjpar = 0;\r
975 for (Int_t ipar=0; ipar<5; ipar++) {\r
976 if(ipar!=1) {\r
977 deltaTouterITS(kipar,0)=outerITSc->GetParameter()[ipar]-trackInnerC2->GetParameter()[ipar];\r
978 deltaouterITS(0,kipar)=outerITSc->GetParameter()[ipar]-trackInnerC2->GetParameter()[ipar];\r
979 }\r
980\r
981 kjpar=0;\r
982 for (Int_t jpar=0; jpar<5; jpar++) {\r
983 Int_t index=outerITSc->GetIndex(ipar,jpar);\r
984 if(ipar !=1 || jpar!=1) {\r
985 covarMouterITS(kipar,kjpar)=outerITSc->GetCovariance()[index]+trackInnerC2->GetCovariance()[index];\r
986 }\r
987 if(jpar!=1) kjpar++;\r
988 }\r
989 if(ipar!=1) kipar++;\r
990 }\r
991\r
992 // chi2 distance ITSout and InnerParams\r
993 TMatrixD covarMInvouterITS = covarMouterITS.Invert();\r
994 TMatrixD mat2outerITS = covarMInvouterITS*deltaTouterITS;\r
995 chi2OuterITS = deltaouterITS*mat2outerITS; \r
996 //chi2OuterITS.Print();\r
997 } \r
998 }\r
999 }\r
1000 }\r
1001\r
1002 //\r
1003 // MC info\r
1004 //\r
1005 TParticle *particle=NULL, *particleTPC=NULL, *particleITS=NULL;\r
1006 TParticle *particleMother=NULL, *particleMotherTPC=NULL, *particleMotherITS=NULL;\r
1007 Int_t mech=-1, mechTPC=-1, mechITS=-1;\r
1008 Bool_t isPrim=kFALSE, isPrimTPC=kFALSE, isPrimITS=kFALSE;\r
1009 Bool_t isFromStrangess=kFALSE, isFromStrangessTPC=kFALSE, isFromStrangessITS=kFALSE;\r
1010 Bool_t isFromConversion=kFALSE, isFromConversionTPC=kFALSE, isFromConversionITS=kFALSE;\r
1011 Bool_t isFromMaterial=kFALSE, isFromMaterialTPC=kFALSE, isFromMaterialITS=kFALSE;\r
1012\r
1013 AliTrackReference *refTPCIn = NULL;\r
1014 AliTrackReference *refITS = NULL;\r
1015\r
1016 Bool_t isOKtrackInnerC3 = kFALSE;\r
1017 AliExternalTrackParam *trackInnerC3 = new AliExternalTrackParam(*(track->GetInnerParam()));\r
1018\r
1019 if(IsUseMCInfo()) \r
1020 {\r
1021 if(!stack) return;\r
1022\r
1023 //\r
1024 // global track\r
1025 //\r
1026 Int_t label = TMath::Abs(track->GetLabel()); \r
1027 particle = stack->Particle(label);\r
1028 if(particle && particle->GetPDG() && particle->GetPDG()->Charge()!=0.)\r
1029 {\r
1030 particleMother = GetMother(particle,stack);\r
1031 mech = particle->GetUniqueID();\r
1032 isPrim = stack->IsPhysicalPrimary(label);\r
1033 isFromStrangess = IsFromStrangeness(label,stack);\r
1034 isFromConversion = IsFromConversion(label,stack);\r
1035 isFromMaterial = IsFromMaterial(label,stack);\r
1036 }\r
1037\r
1038 //\r
1039 // TPC track\r
1040 //\r
1041 Int_t labelTPC = TMath::Abs(track->GetTPCLabel()); \r
1042 particleTPC = stack->Particle(labelTPC);\r
1043 if(particleTPC && particleTPC->GetPDG() && particleTPC->GetPDG()->Charge()!=0.)\r
1044 {\r
1045 particleMotherTPC = GetMother(particleTPC,stack);\r
1046 mechTPC = particleTPC->GetUniqueID();\r
1047 isPrimTPC = stack->IsPhysicalPrimary(labelTPC);\r
1048 isFromStrangessTPC = IsFromStrangeness(labelTPC,stack);\r
1049 isFromConversionTPC = IsFromConversion(labelTPC,stack);\r
1050 isFromMaterialTPC = IsFromMaterial(labelTPC,stack);\r
1051 }\r
1052\r
1053 //\r
1054 // store first track reference\r
1055 // for TPC track\r
1056 //\r
1057 TParticle *part=0;\r
1058 TClonesArray *trefs=0;\r
1059 Int_t status = mcEvent->GetParticleAndTR(track->GetTPCLabel(), part, trefs);\r
1060\r
1061 if(status>0 && part && trefs && part->GetPDG() && part->GetPDG()->Charge()!=0.) \r
1062 {\r
1063 Int_t nTrackRef = trefs->GetEntries();\r
1064 //printf("nTrackRef %d \n",nTrackRef);\r
1065\r
1066 Int_t countITS = 0;\r
1067 for (Int_t iref = 0; iref < nTrackRef; iref++) \r
1068 {\r
1069 AliTrackReference *ref = (AliTrackReference *)trefs->At(iref);\r
1070\r
1071 // Ref. in the middle ITS \r
1072 if(ref && ref->DetectorId()==AliTrackReference::kITS)\r
1073 {\r
1074 if(!refITS && countITS==2) {\r
1075 refITS = ref;\r
1076 //printf("refITS %p \n",refITS);\r
1077 }\r
1078 countITS++;\r
1079 }\r
1080\r
1081 // TPC\r
1082 if(ref && ref->DetectorId()==AliTrackReference::kTPC)\r
1083 {\r
1084 if(!refTPCIn) {\r
1085 refTPCIn = ref;\r
1086 //printf("refTPCIn %p \n",refTPCIn);\r
1087 //break;\r
1088 }\r
1089 }\r
1090 }\r
1091\r
1092 // transform inner params to TrackRef\r
1093 // reference frame\r
1094 if(refTPCIn && trackInnerC3) \r
1095 {\r
1096 Double_t kRefPhi = TMath::ATan2(refTPCIn->Y(),refTPCIn->X());\r
1097 isOKtrackInnerC3 = trackInnerC3->Rotate(kRefPhi);\r
1098 isOKtrackInnerC3 = AliTracker::PropagateTrackToBxByBz(trackInnerC3,refTPCIn->R(),track->GetMass(),kStep,kFALSE);\r
1099 }\r
1100 }\r
1101\r
1102 //\r
1103 // ITS track\r
1104 //\r
1105 Int_t labelITS = TMath::Abs(track->GetITSLabel()); \r
1106 particleITS = stack->Particle(labelITS);\r
1107 if(particleITS && particleITS->GetPDG() && particleITS->GetPDG()->Charge()!=0.)\r
1108 {\r
1109 particleMotherITS = GetMother(particleITS,stack);\r
1110 mechITS = particleITS->GetUniqueID();\r
1111 isPrimITS = stack->IsPhysicalPrimary(labelITS);\r
1112 isFromStrangessITS = IsFromStrangeness(labelITS,stack);\r
1113 isFromConversionITS = IsFromConversion(labelITS,stack);\r
1114 isFromMaterialITS = IsFromMaterial(labelITS,stack);\r
1115 }\r
1116 }\r
1117\r
1118 //\r
1119 Bool_t dumpToTree=kFALSE;\r
1120 \r
1121 if(isOKtpcInnerC && isOKtrackInnerC) dumpToTree = kTRUE;\r
1122 if(fUseESDfriends && isOKtrackInnerC2 && isOKouterITSc) dumpToTree = kTRUE;\r
1123 if(fUseMCInfo && isOKtrackInnerC3) dumpToTree = kTRUE;\r
04fe1e30 1124 TObjString triggerClass = esdEvent->GetFiredTriggerClasses().Data();\r
1059f583 1125\r
a1419baa 1126 /////////////////\r
1127 //book keeping of created dummy objects (to avoid NULL in trees)\r
1128 Bool_t newvtxESD=kFALSE;\r
1129 Bool_t newtrack=kFALSE;\r
1130 Bool_t newtpcInnerC=kFALSE;\r
1131 Bool_t newtrackInnerC=kFALSE;\r
1132 Bool_t newtrackInnerC2=kFALSE;\r
1133 Bool_t newouterITSc=kFALSE;\r
1134 Bool_t newtrackInnerC3=kFALSE;\r
1135 Bool_t newrefTPCIn=kFALSE;\r
1136 Bool_t newrefITS=kFALSE;\r
1137 Bool_t newparticle=kFALSE;\r
1138 Bool_t newparticleMother=kFALSE;\r
1139 Bool_t newparticleTPC=kFALSE;\r
1140 Bool_t newparticleMotherTPC=kFALSE;\r
1141 Bool_t newparticleITS=kFALSE;\r
1142 Bool_t newparticleMotherITS=kFALSE;\r
1143 \r
1144 //check that the vertex is there and that it is OK, \r
1145 //i.e. no null member arrays, otherwise a problem with merging\r
1146 //later on.\r
1147 //this is a very ugly hack!\r
1148 if (!vtxESD)\r
1149 {\r
1150 AliInfo("fixing the ESD vertex for streaming");\r
1151 vtxESD=new AliESDVertex(); \r
1152 //vtxESD->SetNContributors(1);\r
1153 //UShort_t pindices[1]; pindices[0]=0;\r
1154 //vtxESD->SetIndices(1,pindices);\r
1155 //vtxESD->SetNContributors(0);\r
1156 newvtxESD=kTRUE;\r
1157 }\r
1059f583 1158 //\r
a1419baa 1159 if (!track) {track=new AliESDtrack();newtrack=kTRUE;}\r
1160 if (!tpcInnerC) {tpcInnerC=new AliExternalTrackParam();newtpcInnerC=kTRUE;}\r
1161 if (!trackInnerC) {trackInnerC=new AliExternalTrackParam();newtrackInnerC=kTRUE;}\r
1162 if (!trackInnerC2) {trackInnerC2=new AliExternalTrackParam();newtrackInnerC2=kTRUE;}\r
1163 if (!outerITSc) {outerITSc=new AliExternalTrackParam();newouterITSc=kTRUE;}\r
1164 if (!trackInnerC3) {trackInnerC3=new AliExternalTrackParam();newtrackInnerC3=kTRUE;}\r
1165 if (fUseMCInfo)\r
1166 {\r
1167 if (!refTPCIn) {refTPCIn=new AliTrackReference(); newrefTPCIn=kTRUE;}\r
1168 if (!refITS) {refITS=new AliTrackReference();newrefITS=kTRUE;}\r
1169 if (!particle) {particle=new TParticle(); newparticle=kTRUE;}\r
1170 if (!particleMother) {particleMother=new TParticle();newparticleMother=kTRUE;}\r
1171 if (!particleTPC) {particleTPC=new TParticle();newparticleTPC=kTRUE;}\r
1172 if (!particleMotherTPC) {particleMotherTPC=new TParticle();newparticleMotherTPC=kTRUE;}\r
1173 if (!particleITS) {particleITS=new TParticle();newparticleITS=kTRUE;}\r
1174 if (!particleMotherITS) {particleMotherITS=new TParticle();newparticleMotherITS=kTRUE;}\r
1175 }\r
1176 /////////////////////////\r
1177\r
1178 //Double_t vtxX=vtxESD->GetX();\r
1179 //Double_t vtxY=vtxESD->GetY();\r
1180 //Double_t vtxZ=vtxESD->GetZ();\r
1181\r
1182 AliESDVertex* pvtxESD = (AliESDVertex*)vtxESD->Clone();\r
1183 AliESDtrack* ptrack=(AliESDtrack*)track->Clone();\r
1184 AliExternalTrackParam* ptpcInnerC = (AliExternalTrackParam*)tpcInnerC->Clone();\r
1185 AliExternalTrackParam* ptrackInnerC = (AliExternalTrackParam*)trackInnerC->Clone();\r
1186 AliExternalTrackParam* ptrackInnerC2 = (AliExternalTrackParam*)trackInnerC2->Clone();\r
1187 AliExternalTrackParam* pouterITSc = (AliExternalTrackParam*)outerITSc->Clone();\r
1188 AliExternalTrackParam* ptrackInnerC3 = (AliExternalTrackParam*)trackInnerC3->Clone();\r
1189\r
1059f583 1190 if(fTreeSRedirector && dumpToTree) \r
1191 {\r
a1419baa 1192\r
1059f583 1193 (*fTreeSRedirector)<<"highPt"<<\r
a1419baa 1194 "runNumber="<<runNumber<<\r
1195 "evtTimeStamp="<<evtTimeStamp<<\r
1196 "evtNumberInFile="<<evtNumberInFile<<\r
1197 "triggerClass="<<&triggerClass<< // trigger\r
1198 "Bz="<<bz<<\r
1199 "vtxESD.="<<pvtxESD<<\r
1200 //"vtxESDx="<<vtxX<<\r
1201 //"vtxESDy="<<vtxY<<\r
1202 //"vtxESDz="<<vtxZ<<\r
1203 "IRtot="<<ir1<<\r
1204 "IRint2="<<ir2<<\r
1205 "mult="<<mult<<\r
1206 "esdTrack.="<<ptrack<<\r
1207 "extTPCInnerC.="<<ptpcInnerC<<\r
1208 "extInnerParamC.="<<ptrackInnerC<<\r
1209 "extInnerParam.="<<ptrackInnerC2<<\r
1210 "extOuterITS.="<<pouterITSc<<\r
1211 "extInnerParamRef.="<<ptrackInnerC3<<\r
1212 "chi2TPCInnerC="<<chi2(0,0)<<\r
1213 "chi2InnerC="<<chi2trackC(0,0)<<\r
1214 "chi2OuterITS="<<chi2OuterITS(0,0)<<\r
1215 "centralityF="<<centralityF;\r
1216 if (fUseMCInfo)\r
1217 {\r
1218 (*fTreeSRedirector)<<"highPt"<<\r
1059f583 1219 "refTPCIn.="<<refTPCIn<<\r
1220 "refITS.="<<refITS<<\r
1059f583 1221 "particle.="<<particle<<\r
a1419baa 1222 "particleMother.="<<particleMother<<\r
1059f583 1223 "mech="<<mech<<\r
1224 "isPrim="<<isPrim<<\r
a1419baa 1225 "isFromStrangess="<<isFromStrangess<<\r
1226 "isFromConversion="<<isFromConversion<<\r
1059f583 1227 "isFromMaterial="<<isFromMaterial<<\r
1228 "particleTPC.="<<particleTPC<<\r
a1419baa 1229 "particleMotherTPC.="<<particleMotherTPC<<\r
1059f583 1230 "mechTPC="<<mechTPC<<\r
1231 "isPrimTPC="<<isPrimTPC<<\r
a1419baa 1232 "isFromStrangessTPC="<<isFromStrangessTPC<<\r
1233 "isFromConversionTPC="<<isFromConversionTPC<<\r
1059f583 1234 "isFromMaterialTPC="<<isFromMaterialTPC<<\r
1235 "particleITS.="<<particleITS<<\r
a1419baa 1236 "particleMotherITS.="<<particleMotherITS<<\r
1059f583 1237 "mechITS="<<mechITS<<\r
1238 "isPrimITS="<<isPrimITS<<\r
a1419baa 1239 "isFromStrangessITS="<<isFromStrangessITS<<\r
1240 "isFromConversionITS="<<isFromConversionITS<<\r
1241 "isFromMaterialITS="<<isFromMaterialITS;\r
1059f583 1242 }\r
a1419baa 1243 //finish writing the entry\r
1244 (*fTreeSRedirector)<<"highPt"<<"\n";\r
1245 }\r
1246\r
1247 delete pvtxESD;\r
1248 delete ptrack;\r
1249 delete ptpcInnerC;\r
1250 delete ptrackInnerC;\r
1251 delete ptrackInnerC2;\r
1252 delete pouterITSc;\r
1253 delete ptrackInnerC3;\r
1254\r
1255 ////////////////////\r
1256 //delete the dummy objects we might have created.\r
1257 if (newvtxESD) {delete vtxESD; vtxESD=NULL;}\r
1258 if (newtrack) {delete track; track=NULL;}\r
1259 if (newtpcInnerC) {delete tpcInnerC; tpcInnerC=NULL;}\r
1260 if (newtrackInnerC) {delete trackInnerC; trackInnerC=NULL;}\r
1261 if (newtrackInnerC2) {delete trackInnerC2; trackInnerC2=NULL;}\r
1262 if (newouterITSc) {delete outerITSc; outerITSc=NULL;}\r
1263 if (newtrackInnerC3) {delete trackInnerC3; trackInnerC3=NULL;}\r
1264 if (newrefTPCIn) {delete refTPCIn; refTPCIn=NULL;}\r
1265 if (newrefITS) {delete refITS; refITS=NULL;}\r
1266 if (newparticle) {delete particle; particle=NULL;}\r
1267 if (newparticleMother) {delete particleMother; particleMother=NULL;}\r
1268 if (newparticleTPC) {delete particleTPC; particleTPC=NULL;}\r
1269 if (newparticleMotherTPC) {delete particleMotherTPC; particleMotherTPC=NULL;}\r
1270 if (newparticleITS) {delete particleITS; particleITS=NULL;}\r
1271 if (newparticleMotherITS) {delete particleMotherITS; particleMotherITS=NULL;}\r
1272 ///////////////\r
1273\r
1274 delete tpcInnerC;\r
1275 delete trackInnerC;\r
1276 delete trackInnerC2;\r
1277 delete outerITSc;\r
1278 delete trackInnerC3;\r
1059f583 1279 }\r
1280 }\r
a1419baa 1281\r
1059f583 1282}\r
1283\r
1284\r
1285//_____________________________________________________________________________\r
1286void AliAnalysisTaskFilteredTree::ProcessMCEff(AliESDEvent *const esdEvent, AliMCEvent * const mcEvent, AliESDfriend *const /*esdFriend*/)\r
1287{\r
1288 //\r
1289 // Fill tree for efficiency studies MC only\r
1290\r
1291 if(!esdEvent) {\r
1292 AliDebug(AliLog::kError, "esdEvent not available");\r
1293 return;\r
1294 }\r
1295\r
1296 if(!mcEvent) {\r
1297 AliDebug(AliLog::kError, "mcEvent not available");\r
1298 return;\r
1299 }\r
1300\r
1301 // get selection cuts\r
1302 AliFilteredTreeEventCuts *evtCuts = GetEventCuts(); \r
1303 AliFilteredTreeAcceptanceCuts *accCuts = GetAcceptanceCuts(); \r
1304 AliESDtrackCuts *esdTrackCuts = GetTrackCuts(); \r
1305\r
1306 if(!evtCuts || !accCuts || !esdTrackCuts) {\r
1307 AliDebug(AliLog::kError, "cuts not available");\r
1308 return;\r
1309 }\r
1310\r
1311 // trigger selection\r
1312 Bool_t isEventTriggered = kTRUE;\r
1313 AliPhysicsSelection *physicsSelection = NULL;\r
1314 AliTriggerAnalysis* triggerAnalysis = NULL;\r
1315\r
1316 // \r
1317 AliInputEventHandler* inputHandler = (AliInputEventHandler*) AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler();\r
1318 if (!inputHandler)\r
1319 {\r
1320 Printf("ERROR: Could not receive input handler");\r
1321 return;\r
1322 }\r
1323 \r
1324 // get file name\r
1325 TTree *chain = (TChain*)GetInputData(0);\r
1326 if(!chain) { \r
1327 Printf("ERROR: Could not receive input chain");\r
1328 return;\r
1329 }\r
1330 TObjString fileName(chain->GetCurrentFile()->GetName());\r
1331\r
1332 // trigger\r
1333 if(evtCuts->IsTriggerRequired()) \r
1334 {\r
1335 // always MB\r
1336 isEventTriggered = inputHandler->IsEventSelected() & AliVEvent::kMB;\r
1337\r
1338 physicsSelection = static_cast<AliPhysicsSelection*> (inputHandler->GetEventSelection());\r
1339 if(!physicsSelection) return;\r
1340\r
1341 if (isEventTriggered && (GetTrigger() == AliTriggerAnalysis::kV0AND)) {\r
1342 // set trigger (V0AND)\r
1343 triggerAnalysis = physicsSelection->GetTriggerAnalysis();\r
1344 if(!triggerAnalysis) return;\r
1345 isEventTriggered = triggerAnalysis->IsOfflineTriggerFired(esdEvent, GetTrigger());\r
1346 }\r
1347 }\r
1348\r
1349 // centrality determination\r
1350 Float_t centralityF = -1;\r
1351 AliCentrality *esdCentrality = esdEvent->GetCentrality();\r
1352 centralityF = esdCentrality->GetCentralityPercentile(fCentralityEstimator.Data());\r
1353\r
1354 // use MC information\r
1355 AliHeader* header = 0;\r
1356 AliGenEventHeader* genHeader = 0;\r
1357 AliStack* stack = 0;\r
1358 TArrayF vtxMC(3);\r
1359\r
1360 Int_t multMCTrueTracks = 0;\r
1361 if(IsUseMCInfo())\r
1362 {\r
1363 //\r
1364 if(!mcEvent) {\r
1365 AliDebug(AliLog::kError, "mcEvent not available");\r
1366 return;\r
1367 }\r
1368 // get MC event header\r
1369 header = mcEvent->Header();\r
1370 if (!header) {\r
1371 AliDebug(AliLog::kError, "Header not available");\r
1372 return;\r
1373 }\r
1374 // MC particle stack\r
1375 stack = mcEvent->Stack();\r
1376 if (!stack) {\r
1377 AliDebug(AliLog::kError, "Stack not available");\r
1378 return;\r
1379 }\r
1380\r
1381 // get MC vertex\r
1382 genHeader = header->GenEventHeader();\r
1383 if (!genHeader) {\r
1384 AliDebug(AliLog::kError, "Could not retrieve genHeader from Header");\r
1385 return;\r
1386 }\r
1387 genHeader->PrimaryVertex(vtxMC);\r
1388\r
1389 // multipliticy of all MC primary tracks\r
1390 // in Zv, pt and eta ranges)\r
1391 multMCTrueTracks = GetMCTrueTrackMult(mcEvent,evtCuts,accCuts);\r
1392\r
1393 } // end bUseMC\r
1394\r
1395 // get reconstructed vertex \r
1396 //const AliESDVertex* vtxESD = 0; \r
1397 AliESDVertex* vtxESD = 0; \r
1398 if(GetAnalysisMode() == kTPCAnalysisMode) {\r
1399 vtxESD = (AliESDVertex*)esdEvent->GetPrimaryVertexTPC();\r
1400 }\r
1401 else if(GetAnalysisMode() == kTPCITSAnalysisMode) {\r
1402 vtxESD = (AliESDVertex*)esdEvent->GetPrimaryVertexTracks();\r
1403 }\r
1404 else {\r
1405 return;\r
1406 }\r
1407\r
1408 if(!vtxESD) return;\r
1409\r
1410 Bool_t isEventOK = evtCuts->AcceptEvent(esdEvent,mcEvent,vtxESD); \r
1411 //printf("isEventOK %d, isEventTriggered %d \n",isEventOK, isEventTriggered);\r
1412 //printf("GetAnalysisMode() %d \n",GetAnalysisMode());\r
1413\r
04fe1e30 1414 TObjString triggerClass = esdEvent->GetFiredTriggerClasses().Data();\r
1415\r
1059f583 1416 // check event cuts\r
1417 if(isEventOK && isEventTriggered)\r
1418 {\r
1419 if(IsUseMCInfo()) \r
1420 {\r
1421 if(!stack) return;\r
1422\r
1423 //\r
1424 // MC info\r
1425 //\r
1426 TParticle *particle=NULL;\r
1427 TParticle *particleMother=NULL;\r
1428 Int_t mech=-1;\r
1429\r
1430 // reco event info\r
1431 Double_t vert[3] = {0}; \r
1432 vert[0] = vtxESD->GetXv();\r
1433 vert[1] = vtxESD->GetYv();\r
1434 vert[2] = vtxESD->GetZv();\r
1435 Int_t mult = vtxESD->GetNContributors();\r
1436 Double_t bz = esdEvent->GetMagneticField();\r
1437 Double_t runNumber = esdEvent->GetRunNumber();\r
1438 Double_t evtTimeStamp = esdEvent->GetTimeStamp();\r
1439 Int_t evtNumberInFile = esdEvent->GetEventNumberInFile();\r
1440\r
1441 // loop over MC stack\r
1442 for (Int_t iMc = 0; iMc < stack->GetNtrack(); ++iMc) \r
1443 {\r
1444 particle = stack->Particle(iMc);\r
1445 if (!particle)\r
1446 continue;\r
1447\r
1448 // only charged particles\r
1449 if(!particle->GetPDG()) continue;\r
1450 Double_t charge = particle->GetPDG()->Charge()/3.;\r
1451 if (TMath::Abs(charge) < 0.001)\r
1452 continue;\r
1453\r
1454 // only primary particles\r
1455 Bool_t prim = stack->IsPhysicalPrimary(iMc);\r
1456 if(!prim) continue;\r
1457\r
1458 // downscale low-pT particles\r
1459 Double_t scalempt= TMath::Min(particle->Pt(),10.);\r
1460 Double_t downscaleF = gRandom->Rndm();\r
1461 downscaleF *= fLowPtTrackDownscaligF;\r
1462 if(TMath::Exp(2*scalempt)<downscaleF) continue;\r
1463\r
1464 // is particle in acceptance\r
1465 if(!accCuts->AcceptTrack(particle)) continue;\r
1466 \r
1467 // check if particle reconstructed\r
1468 Bool_t isRec = kFALSE;\r
1469 Int_t trackIndex = -1;\r
1470 for (Int_t iTrack = 0; iTrack < esdEvent->GetNumberOfTracks(); iTrack++)\r
1471 {\r
1472 \r
1473 AliESDtrack *track = esdEvent->GetTrack(iTrack);\r
1474 if(!track) continue;\r
1475 if(track->Charge()==0) continue;\r
1476 if(esdTrackCuts->AcceptTrack(track) && accCuts->AcceptTrack(track)) \r
1477 {\r
1478 Int_t label = TMath::Abs(track->GetLabel());\r
1479 if(label == iMc) {\r
1480 isRec = kTRUE;\r
1481 trackIndex = iTrack;\r
1482 break;\r
1483 }\r
1484 } \r
1485 }\r
1486\r
1487 // Store information in the output tree\r
1488 AliESDtrack *recTrack = NULL; \r
1489 if(trackIndex>-1) { \r
1490 recTrack = esdEvent->GetTrack(trackIndex); \r
1491 } else {\r
1492 recTrack = new AliESDtrack(); \r
1493 } \r
1494\r
1495 particleMother = GetMother(particle,stack);\r
1496 mech = particle->GetUniqueID();\r
1497\r
1498 //MC particle track length\r
1499 Double_t tpcTrackLength = 0.;\r
1500 AliMCParticle *mcParticle = (AliMCParticle*) mcEvent->GetTrack(iMc);\r
1501 if(mcParticle) {\r
1502 Int_t counter;\r
1503 tpcTrackLength = mcParticle->GetTPCTrackLength(bz,0.05,counter,3.0);\r
1504 } \r
1505\r
1506\r
1507 //\r
1508 if(fTreeSRedirector) {\r
1509 (*fTreeSRedirector)<<"MCEffTree"<<\r
1510 "fileName.="<<&fileName<<\r
04fe1e30 1511 "triggerClass.="<<&triggerClass<<\r
1059f583 1512 "runNumber="<<runNumber<<\r
1513 "evtTimeStamp="<<evtTimeStamp<<\r
1514 "evtNumberInFile="<<evtNumberInFile<<\r
1515 "Bz="<<bz<<\r
1516 "vtxESD.="<<vtxESD<<\r
1517 "mult="<<mult<<\r
1518 "esdTrack.="<<recTrack<<\r
1519 "isRec="<<isRec<<\r
1520 "tpcTrackLength="<<tpcTrackLength<<\r
1521 "particle.="<<particle<<\r
1522 "particleMother.="<<particleMother<<\r
1523 "mech="<<mech<<\r
1524 "\n";\r
1525 }\r
1526\r
1527 if(trackIndex <0 && recTrack) delete recTrack; recTrack=0;\r
1528 }\r
1529 }\r
1530 }\r
1531 \r
1532}\r
1533\r
1534//_____________________________________________________________________________\r
1535Bool_t AliAnalysisTaskFilteredTree::IsHighDeDxParticle(AliESDtrack * track) {\r
1536 //\r
1537 // check if particle is Z > 1 \r
1538 //\r
1539 if (track->GetTPCNcls() < 60) return kFALSE;\r
1540 Double_t mom = track->GetInnerParam()->GetP();\r
1541 if (mom < 0.2) return kFALSE; // protection against unexpected behavior of Aleph parameterization\r
1542 Float_t dca[2], bCov[3];\r
1543 track->GetImpactParameters(dca,bCov);\r
1544 //\r
1545\r
1546 Double_t triggerDeDx = 4*AliExternalTrackParam::BetheBlochAleph((mom*2)/(0.938*3),1.0288,31.9806,5.04114e-11,2.13096,2.38541);\r
1547\r
1548 if (track->GetTPCsignal() > triggerDeDx && track->GetTPCsignal()<1000 && TMath::Abs(dca[0])<3.) return kTRUE;\r
1549\r
1550 return kFALSE;\r
1551}\r
1552\r
1553//_____________________________________________________________________________\r
1554void AliAnalysisTaskFilteredTree::ProcessV0(AliESDEvent *const esdEvent, AliMCEvent * const mcEvent, AliESDfriend *const /*esdFriend*/)\r
1555{\r
1556 //\r
1557 // Select real events with V0 (K0s and Lambda) high-pT candidates\r
1558 //\r
1559 if(!esdEvent) {\r
1560 AliDebug(AliLog::kError, "esdEvent not available");\r
1561 return;\r
1562 }\r
1563\r
1564 // get selection cuts\r
1565 AliFilteredTreeEventCuts *evtCuts = GetEventCuts(); \r
1566 AliFilteredTreeAcceptanceCuts *accCuts = GetAcceptanceCuts(); \r
1567 AliESDtrackCuts *esdTrackCuts = GetTrackCuts(); \r
1568\r
1569 if(!evtCuts || !accCuts || !esdTrackCuts) {\r
1570 AliDebug(AliLog::kError, "cuts not available");\r
1571 return;\r
1572 }\r
1573\r
1574 // trigger selection\r
1575 Bool_t isEventTriggered = kTRUE;\r
1576 AliPhysicsSelection *physicsSelection = NULL;\r
1577 AliTriggerAnalysis* triggerAnalysis = NULL;\r
1578\r
1579 // \r
1580 AliInputEventHandler* inputHandler = (AliInputEventHandler*) AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler();\r
1581 if (!inputHandler)\r
1582 {\r
1583 Printf("ERROR: Could not receive input handler");\r
1584 return;\r
1585 }\r
1586 \r
1587 // get file name\r
1588 TTree *chain = (TChain*)GetInputData(0);\r
1589 if(!chain) { \r
1590 Printf("ERROR: Could not receive input chain");\r
1591 return;\r
1592 }\r
1593 TObjString fileName(chain->GetCurrentFile()->GetName());\r
1594\r
1595 // trigger\r
1596 if(evtCuts->IsTriggerRequired()) \r
1597 {\r
1598 // always MB\r
1599 isEventTriggered = inputHandler->IsEventSelected() & AliVEvent::kMB;\r
1600\r
1601 physicsSelection = static_cast<AliPhysicsSelection*> (inputHandler->GetEventSelection());\r
1602 if(!physicsSelection) return;\r
1603 //SetPhysicsTriggerSelection(physicsSelection);\r
1604\r
1605 if (isEventTriggered && (GetTrigger() == AliTriggerAnalysis::kV0AND)) {\r
1606 // set trigger (V0AND)\r
1607 triggerAnalysis = physicsSelection->GetTriggerAnalysis();\r
1608 if(!triggerAnalysis) return;\r
1609 isEventTriggered = triggerAnalysis->IsOfflineTriggerFired(esdEvent, GetTrigger());\r
1610 }\r
1611 }\r
1612\r
1613 // centrality determination\r
1614 Float_t centralityF = -1;\r
1615 AliCentrality *esdCentrality = esdEvent->GetCentrality();\r
1616 centralityF = esdCentrality->GetCentralityPercentile(fCentralityEstimator.Data());\r
1617\r
1618\r
1619 // get reconstructed vertex \r
1620 //const AliESDVertex* vtxESD = 0; \r
1621 AliESDVertex* vtxESD = 0; \r
1622 if(GetAnalysisMode() == kTPCAnalysisMode) {\r
1623 vtxESD = (AliESDVertex*)esdEvent->GetPrimaryVertexTPC();\r
1624 }\r
1625 else if(GetAnalysisMode() == kTPCITSAnalysisMode) {\r
1626 vtxESD = (AliESDVertex*)esdEvent->GetPrimaryVertexTracks();\r
1627 }\r
1628 else {\r
1629 return;\r
1630 }\r
1631\r
1632 if(!vtxESD) return;\r
1633\r
1634 Bool_t isEventOK = evtCuts->AcceptEvent(esdEvent,mcEvent,vtxESD); \r
1635 //printf("isEventOK %d, isEventTriggered %d \n",isEventOK, isEventTriggered);\r
1636 //printf("GetAnalysisMode() %d \n",GetAnalysisMode());\r
1637\r
1638 // check event cuts\r
1639 if(isEventOK && isEventTriggered) {\r
1640 //\r
1641 // Dump the pt downscaled V0 into the tree\r
1642 // \r
1643 Int_t ntracks = esdEvent->GetNumberOfTracks();\r
1644 Int_t nV0s = esdEvent->GetNumberOfV0s();\r
1645 Int_t run = esdEvent->GetRunNumber();\r
1646 Int_t time= esdEvent->GetTimeStamp();\r
1647 Int_t evNr=esdEvent->GetEventNumberInFile();\r
1648 Double_t bz = esdEvent->GetMagneticField();\r
1649\r
1650\r
1651 for (Int_t iv0=0; iv0<nV0s; iv0++){\r
1652 AliESDv0 * v0 = esdEvent->GetV0(iv0);\r
1653 if (!v0) continue;\r
1654 AliESDtrack * track0 = esdEvent->GetTrack(v0->GetIndex(0));\r
1655 AliESDtrack * track1 = esdEvent->GetTrack(v0->GetIndex(1));\r
1656 if (!track0) continue;\r
1657 if (!track1) continue;\r
1658 if (track0->GetSign()<0) {\r
1659 track1 = esdEvent->GetTrack(v0->GetIndex(0));\r
1660 track0 = esdEvent->GetTrack(v0->GetIndex(1));\r
1661 }\r
1662 //\r
1663 Bool_t isDownscaled = IsV0Downscaled(v0);\r
1664 if (isDownscaled) continue;\r
1665 AliKFParticle kfparticle; //\r
1666 Int_t type=GetKFParticle(v0,esdEvent,kfparticle);\r
1667 if (type==0) continue; \r
04fe1e30 1668 TObjString triggerClass = esdEvent->GetFiredTriggerClasses().Data();\r
1059f583 1669\r
1670 if(!fTreeSRedirector) return;\r
1671 (*fTreeSRedirector)<<"V0s"<<\r
1672 "isDownscaled="<<isDownscaled<<\r
04fe1e30 1673 "triggerClass="<<&triggerClass<< // trigger\r
1059f583 1674 "Bz="<<bz<<\r
1675 "fileName.="<<&fileName<<\r
1676 "runNumber="<<run<<\r
1677 "evtTimeStamp="<<time<<\r
1678 "evtNumberInFile="<<evNr<<\r
1679 "type="<<type<<\r
1680 "ntracks="<<ntracks<<\r
1681 "v0.="<<v0<<\r
1682 "kf.="<<&kfparticle<<\r
1683 "track0.="<<track0<<\r
1684 "track1.="<<track1<<\r
1685 "centralityF="<<centralityF<<\r
1686 "\n";\r
1687 }\r
1688 }\r
1689}\r
1690\r
1691//_____________________________________________________________________________\r
1692void AliAnalysisTaskFilteredTree::ProcessdEdx(AliESDEvent *const esdEvent, AliMCEvent * const mcEvent, AliESDfriend *const /*esdFriend*/)\r
1693{\r
1694 //\r
1695 // Select real events with large TPC dEdx signal\r
1696 //\r
1697 if(!esdEvent) {\r
1698 AliDebug(AliLog::kError, "esdEvent not available");\r
1699 return;\r
1700 }\r
1701\r
1702 // get selection cuts\r
1703 AliFilteredTreeEventCuts *evtCuts = GetEventCuts(); \r
1704 AliFilteredTreeAcceptanceCuts *accCuts = GetAcceptanceCuts(); \r
1705 AliESDtrackCuts *esdTrackCuts = GetTrackCuts(); \r
1706\r
1707 if(!evtCuts || !accCuts || !esdTrackCuts) {\r
1708 AliDebug(AliLog::kError, "cuts not available");\r
1709 return;\r
1710 }\r
1711\r
1712 // get file name\r
1713 TTree *chain = (TChain*)GetInputData(0);\r
1714 if(!chain) { \r
1715 Printf("ERROR: Could not receive input chain");\r
1716 return;\r
1717 }\r
1718 TObjString fileName(chain->GetCurrentFile()->GetName());\r
1719\r
1720 // trigger\r
1721 Bool_t isEventTriggered = kTRUE;\r
1722 AliPhysicsSelection *physicsSelection = NULL;\r
1723 AliTriggerAnalysis* triggerAnalysis = NULL;\r
1724\r
1725 // \r
1726 AliInputEventHandler* inputHandler = (AliInputEventHandler*) AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler();\r
1727 if (!inputHandler)\r
1728 {\r
1729 Printf("ERROR: Could not receive input handler");\r
1730 return;\r
1731 }\r
1732\r
1733 if(evtCuts->IsTriggerRequired()) \r
1734 {\r
1735 // always MB\r
1736 isEventTriggered = inputHandler->IsEventSelected() & AliVEvent::kMB;\r
1737\r
1738 physicsSelection = static_cast<AliPhysicsSelection*> (inputHandler->GetEventSelection());\r
1739 if(!physicsSelection) return;\r
1740\r
1741 if (isEventTriggered && (GetTrigger() == AliTriggerAnalysis::kV0AND)) {\r
1742 // set trigger (V0AND)\r
1743 triggerAnalysis = physicsSelection->GetTriggerAnalysis();\r
1744 if(!triggerAnalysis) return;\r
1745 isEventTriggered = triggerAnalysis->IsOfflineTriggerFired(esdEvent, GetTrigger());\r
1746 }\r
1747 }\r
1748\r
1749 // get reconstructed vertex \r
1750 AliESDVertex* vtxESD = 0; \r
1751 if(GetAnalysisMode() == kTPCAnalysisMode) {\r
1752 vtxESD = (AliESDVertex*)esdEvent->GetPrimaryVertexTPC();\r
1753 }\r
1754 else if(GetAnalysisMode() == kTPCITSAnalysisMode) {\r
1755 vtxESD = (AliESDVertex*)esdEvent->GetPrimaryVertexTracks();\r
1756 }\r
1757 else {\r
1758 return;\r
1759 }\r
1760 if(!vtxESD) return;\r
1761\r
1762 Bool_t isEventOK = evtCuts->AcceptEvent(esdEvent,mcEvent,vtxESD); \r
1763 //printf("isEventOK %d, isEventTriggered %d \n",isEventOK, isEventTriggered);\r
1764 //printf("GetAnalysisMode() %d \n",GetAnalysisMode());\r
1765\r
1766\r
1767 // check event cuts\r
1768 if(isEventOK && isEventTriggered)\r
1769 {\r
1770 Double_t vert[3] = {0}; \r
1771 vert[0] = vtxESD->GetXv();\r
1772 vert[1] = vtxESD->GetYv();\r
1773 vert[2] = vtxESD->GetZv();\r
1774 Int_t mult = vtxESD->GetNContributors();\r
1775 Double_t bz = esdEvent->GetMagneticField();\r
1776 Double_t runNumber = esdEvent->GetRunNumber();\r
1777 Double_t evtTimeStamp = esdEvent->GetTimeStamp();\r
1778 Int_t evtNumberInFile = esdEvent->GetEventNumberInFile();\r
1779\r
1780 // large dEdx\r
1781 for (Int_t iTrack = 0; iTrack < esdEvent->GetNumberOfTracks(); iTrack++)\r
1782 {\r
1783 AliESDtrack *track = esdEvent->GetTrack(iTrack);\r
1784 if(!track) continue;\r
1785 if(track->Charge()==0) continue;\r
1786 if(!esdTrackCuts->AcceptTrack(track)) continue;\r
1787 if(!accCuts->AcceptTrack(track)) continue;\r
1788\r
1789 if(!IsHighDeDxParticle(track)) continue;\r
04fe1e30 1790 TObjString triggerClass = esdEvent->GetFiredTriggerClasses().Data();\r
1791\r
1059f583 1792 if(!fTreeSRedirector) return;\r
1793 (*fTreeSRedirector)<<"dEdx"<<\r
1794 "fileName.="<<&fileName<<\r
1795 "runNumber="<<runNumber<<\r
1796 "evtTimeStamp="<<evtTimeStamp<<\r
1797 "evtNumberInFile="<<evtNumberInFile<<\r
04fe1e30 1798 "triggerClass="<<&triggerClass<< // trigger\r
1059f583 1799 "Bz="<<bz<<\r
1800 "vtxESD.="<<vtxESD<<\r
1801 "mult="<<mult<<\r
1802 "esdTrack.="<<track<<\r
1803 "\n";\r
1804 }\r
1805 }\r
1806}\r
1807\r
1808//_____________________________________________________________________________\r
1809Int_t AliAnalysisTaskFilteredTree::GetKFParticle(AliESDv0 *const v0, AliESDEvent * const event, AliKFParticle & kfparticle)\r
1810{\r
1811 //\r
1812 // Create KF particle in case the V0 fullfill selection criteria\r
1813 //\r
1814 // Selection criteria\r
1815 // 0. algorithm cut\r
1816 // 1. track cut\r
1817 // 3. chi2 cut\r
1818 // 4. rough mass cut\r
1819 // 5. Normalized pointing angle cut\r
1820 //\r
1821 const Double_t cutMass=0.2;\r
1822 const Double_t kSigmaDCACut=3;\r
1823 //\r
1824 // 0.) algo cut - accept only on the fly\r
1825 //\r
1826 if (v0->GetOnFlyStatus() ==kFALSE) return 0; \r
1827 //\r
1828 // 1.) track cut\r
1829 // \r
1830 AliESDtrack * track0 = event->GetTrack(v0->GetIndex(0));\r
1831 AliESDtrack * track1 = event->GetTrack(v0->GetIndex(1));\r
1832 /*\r
1833 TCut cutD="abs(track0.fD/sqrt(track0.fCdd))>2&&abs(track1.fD/sqrt(track1.fCdd))>2";\r
1834 TCut cutTheta="abs(track0.fP[3])<1&&abs(track1.fP[3])<1";\r
1835 TCut cutNcl="track0.GetTPCClusterInfo(2,1)>100&&track1.GetTPCClusterInfo(2,1)>100";\r
1836 */ \r
1837 if (TMath::Abs(track0->GetTgl())>1) return 0;\r
1838 if (TMath::Abs(track1->GetTgl())>1) return 0;\r
1839 if ((track0->GetTPCClusterInfo(2,1))<100) return 0;\r
1840 if ((track1->GetTPCClusterInfo(2,1))<100) return 0;\r
1841 //if ((track0->GetITSclusters(0))<2) return 0;\r
1842 //if ((track1->GetITSclusters(0))<2) return 0; \r
1843 Float_t pos0[2]={0}, cov0[3]={0};\r
1844 Float_t pos1[2]={0}, cov1[3]={0};\r
1845 track0->GetImpactParameters(pos0,cov0);\r
1846 track0->GetImpactParameters(pos1,cov1);\r
1847 //\r
1848 if (TMath::Abs(pos0[0])<kSigmaDCACut*TMath::Sqrt(cov0[0])) return 0;\r
1849 if (TMath::Abs(pos1[0])<kSigmaDCACut*TMath::Sqrt(cov1[0])) return 0;\r
1850 // \r
1851 //\r
1852 // 3.) Chi2 cut\r
1853 //\r
1854 Double_t chi2KF = v0->GetKFInfo(2,2,2);\r
1855 if (chi2KF>25) return 0;\r
1856 //\r
1857 // 4.) Rough mass cut - 0.200 GeV\r
1858 //\r
1859 static Double_t masses[2]={-1};\r
1860 if (masses[0]<0){\r
1861 masses[0] = TDatabasePDG::Instance()->GetParticle("K_S0")->Mass();\r
1862 masses[1] = TDatabasePDG::Instance()->GetParticle("Lambda0")->Mass();\r
1863 }\r
1864 Double_t mass00= v0->GetEffMass(0,0);\r
1865 Double_t mass22= v0->GetEffMass(2,2);\r
1866 Double_t mass42= v0->GetEffMass(4,2);\r
1867 Double_t mass24= v0->GetEffMass(2,4);\r
1868 Bool_t massOK=kFALSE;\r
1869 Int_t type=0;\r
1870 Int_t ptype=0;\r
1871 Double_t dmass=1;\r
1872 Int_t p1=0, p2=0;\r
1873 if (TMath::Abs(mass00-0)<cutMass) {\r
1874 massOK=kTRUE; type+=1; \r
1875 if (TMath::Abs(mass00-0)<dmass) {\r
1876 ptype=1;\r
1877 dmass=TMath::Abs(mass00-0); \r
1878 p1=0; p2=0;\r
1879 } \r
1880 }\r
1881 if (TMath::Abs(mass24-masses[1])<cutMass) {\r
1882 massOK=kTRUE; type+=2; \r
1883 if (TMath::Abs(mass24-masses[1])<dmass){\r
1884 dmass = TMath::Abs(mass24-masses[1]);\r
1885 ptype=2;\r
1886 p1=2; p2=4;\r
1887 }\r
1888 }\r
1889 if (TMath::Abs(mass42-masses[1])<cutMass) {\r
1890 massOK=kTRUE; type+=4;\r
1891 if (TMath::Abs(mass42-masses[1])<dmass){\r
1892 dmass = TMath::Abs(mass42-masses[1]);\r
1893 ptype=4;\r
1894 p1=4; p2=2;\r
1895 }\r
1896 }\r
1897 if (TMath::Abs(mass22-masses[0])<cutMass) {\r
1898 massOK=kTRUE; type+=8;\r
1899 if (TMath::Abs(mass22-masses[0])<dmass){\r
1900 dmass = TMath::Abs(mass22-masses[0]);\r
1901 ptype=8;\r
1902 p1=2; p2=2;\r
1903 }\r
1904 }\r
1905 if (type==0) return 0;\r
1906 //\r
1907 const Int_t spdg[5]={kPositron,kMuonPlus,kPiPlus, kKPlus, kProton};\r
1908 const AliExternalTrackParam *paramP = v0->GetParamP();\r
1909 const AliExternalTrackParam *paramN = v0->GetParamN();\r
1910 if (paramP->GetSign()<0){\r
1911 paramP=v0->GetParamP();\r
1912 paramN=v0->GetParamN();\r
1913 }\r
1914 //Double_t *pparam1 = (Double_t*)paramP->GetParameter();\r
1915 //Double_t *pparam2 = (Double_t*)paramN->GetParameter();\r
1916 //\r
1917 AliKFParticle kfp1( *paramP, spdg[p1] );\r
1918 AliKFParticle kfp2( *paramN, -1 *spdg[p2] );\r
1919 AliKFParticle V0KF;\r
1920 (V0KF)+=kfp1;\r
1921 (V0KF)+=kfp2;\r
1922 kfparticle=V0KF;\r
1923 //\r
1924 // Pointing angle\r
1925 //\r
1926 Double_t errPhi = V0KF.GetErrPhi();\r
1927 Double_t pointAngle= TMath::ACos(v0->GetV0CosineOfPointingAngle());\r
1928 if (pointAngle/errPhi>10) return 0; \r
1929 //\r
1930 return ptype; \r
1931}\r
1932\r
1933//_____________________________________________________________________________\r
1934Bool_t AliAnalysisTaskFilteredTree::IsV0Downscaled(AliESDv0 *const v0)\r
1935{\r
1936 //\r
1937 // Downscale randomly low pt V0\r
1938 //\r
1939 //return kFALSE;\r
1940 Double_t maxPt= TMath::Max(v0->GetParamP()->Pt(), v0->GetParamN()->Pt());\r
1941 Double_t scalempt= TMath::Min(maxPt,10.);\r
1942 Double_t downscaleF = gRandom->Rndm();\r
1943 downscaleF *= fLowPtV0DownscaligF;\r
1944 \r
1945 //printf("V0 TMath::Exp(2*scalempt) %e, downscaleF %e \n",TMath::Exp(2*scalempt), downscaleF);\r
1946 if (TMath::Exp(2*scalempt)<downscaleF) return kTRUE;\r
1947 return kFALSE;\r
1948\r
1949 /*\r
1950 TH1F his1("his1","his1",100,0,10);\r
1951 TH1F his2("his2","his2",100,0,10);\r
1952 {for (Int_t i=0; i<10000; i++){\r
1953 Double_t rnd=gRandom->Exp(1);\r
1954 Bool_t isDownscaled =TMath::Exp(rnd)<100*gRandom->Rndm();\r
1955 his1->Fill(rnd); \r
1956 if (!isDownscaled) his2->Fill(rnd); \r
1957 }}\r
1958\r
1959 */\r
1960\r
1961}\r
1962\r
1963\r
1964\r
1965//_____________________________________________________________________________\r
1966Bool_t AliAnalysisTaskFilteredTree::ConstrainTPCInner(AliExternalTrackParam *const tpcInnerC, const AliESDVertex* vtx, Double_t b[3])\r
1967{\r
1968 // Constrain TPC inner params constrained\r
1969 //\r
1970 if(!tpcInnerC) return kFALSE; \r
1971 if(!vtx) return kFALSE;\r
1972\r
1973 Double_t dz[2],cov[3];\r
1974 //AliESDVertex *vtx= (AliESDVertex *)esdEvent->GetPrimaryVertex();\r
1975 //if(!tpcInnerC->PropagateToDCA(vtx, esdEvent->GetMagneticField(), 3, dz, cov)) return kFALSE; \r
1976 //if(!tpcInnerC->PropagateToDCA(vtx, Bz, 3, dz, cov)) return kFALSE; \r
1977 if(!tpcInnerC->PropagateToDCABxByBz(vtx, b, 3, dz, cov)) return kFALSE; \r
1978\r
1979\r
1980 Double_t covar[6]; vtx->GetCovMatrix(covar);\r
1981 Double_t p[2]={tpcInnerC->GetParameter()[0]-dz[0],tpcInnerC->GetParameter()[1]-dz[1]};\r
1982 Double_t c[3]={covar[2],0.,covar[5]};\r
1983 Double_t chi2C=tpcInnerC->GetPredictedChi2(p,c);\r
1984 if (chi2C>kVeryBig) return kFALSE; \r
1985\r
1986 if(!tpcInnerC->Update(p,c)) return kFALSE;\r
1987\r
1988 return kTRUE;\r
1989}\r
1990\r
1991//_____________________________________________________________________________\r
1992Bool_t AliAnalysisTaskFilteredTree::ConstrainTrackInner(AliExternalTrackParam *const trackInnerC, const AliESDVertex* vtx, Double_t mass, Double_t b[3])\r
1993{\r
1994 // Constrain TPC inner params constrained\r
1995 //\r
1996 if(!trackInnerC) return kFALSE; \r
1997 if(!vtx) return kFALSE;\r
1998\r
1999 const Double_t kRadius = 2.8; \r
2000 const Double_t kMaxStep = 1.0; \r
2001\r
2002 Double_t dz[2],cov[3];\r
2003\r
2004 //AliESDVertex *vtx= (AliESDVertex *)esdEvent->GetPrimaryVertex();\r
2005 //if(!trackInnerC->PropagateToDCA(vtx, esdEvent->GetMagneticField(), 3, dz, cov)) return kFALSE; \r
2006 //if(!trackInnerC->PropagateToDCA(vtx, Bz, 3, dz, cov)) return kFALSE; \r
2007\r
2008 if(!AliTracker::PropagateTrackToBxByBz(trackInnerC,kRadius,mass,kMaxStep,kFALSE)) return kFALSE;\r
2009 if(!trackInnerC->PropagateToDCABxByBz(vtx, b, 3, dz, cov)) return kFALSE; \r
2010\r
2011 //\r
2012 Double_t covar[6]; vtx->GetCovMatrix(covar);\r
2013 Double_t p[2]={trackInnerC->GetParameter()[0]-dz[0],trackInnerC->GetParameter()[1]-dz[1]};\r
2014 Double_t c[3]={covar[2],0.,covar[5]};\r
2015 Double_t chi2C=trackInnerC->GetPredictedChi2(p,c);\r
2016 if (chi2C>kVeryBig) return kFALSE; \r
2017\r
2018 if(!trackInnerC->Update(p,c)) return kFALSE;\r
2019\r
2020 return kTRUE;\r
2021}\r
2022\r
2023\r
2024//_____________________________________________________________________________\r
2025TParticle *AliAnalysisTaskFilteredTree::GetMother(TParticle *const particle, AliStack *const stack) \r
2026{\r
2027 if(!particle) return NULL;\r
2028 if(!stack) return NULL;\r
2029\r
2030 Int_t motherLabel = TMath::Abs(particle->GetMother(0)); \r
2031 TParticle* mother = NULL; \r
2032 mother = stack->Particle(motherLabel); \r
2033\r
2034return mother;\r
2035}\r
2036\r
2037//_____________________________________________________________________________\r
2038Bool_t AliAnalysisTaskFilteredTree::IsFromConversion(const Int_t label, AliStack *const stack) \r
2039{\r
2040 Bool_t isFromConversion = kFALSE;\r
2041\r
2042 if(stack) {\r
2043 TParticle* particle = stack->Particle(label);\r
2044\r
2045 if(particle && particle->GetPDG() && particle->GetPDG()->Charge()!=0) \r
2046 {\r
2047 Int_t mech = particle->GetUniqueID(); // production mechanism \r
2048 Bool_t isPrim = stack->IsPhysicalPrimary(label);\r
2049\r
2050 Int_t motherLabel = TMath::Abs(particle->GetMother(0)); \r
2051 TParticle* mother = stack->Particle(motherLabel); \r
2052 if(mother) {\r
2053 Int_t motherPdg = mother->GetPdgCode();\r
2054\r
2055 if(!isPrim && mech==5 && motherPdg==kGamma) { \r
2056 isFromConversion=kTRUE; \r
2057 }\r
2058 }\r
2059 } \r
2060 } \r
2061\r
2062 return isFromConversion;\r
2063}\r
2064\r
2065//_____________________________________________________________________________\r
2066Bool_t AliAnalysisTaskFilteredTree::IsFromMaterial(const Int_t label, AliStack *const stack) \r
2067{\r
2068 Bool_t isFromMaterial = kFALSE;\r
2069\r
2070 if(stack) {\r
2071 TParticle* particle = stack->Particle(label);\r
2072\r
2073 if(particle && particle->GetPDG() && particle->GetPDG()->Charge()!=0) \r
2074 {\r
2075 Int_t mech = particle->GetUniqueID(); // production mechanism \r
2076 Bool_t isPrim = stack->IsPhysicalPrimary(label);\r
2077\r
2078 Int_t motherLabel = TMath::Abs(particle->GetMother(0)); \r
2079 TParticle* mother = stack->Particle(motherLabel); \r
2080 if(mother) {\r
2081 if(!isPrim && mech==13) { \r
2082 isFromMaterial=kTRUE; \r
2083 }\r
2084 }\r
2085 } \r
2086 } \r
2087\r
2088 return isFromMaterial;\r
2089}\r
2090\r
2091//_____________________________________________________________________________\r
2092Bool_t AliAnalysisTaskFilteredTree::IsFromStrangeness(const Int_t label, AliStack *const stack) \r
2093{\r
2094 Bool_t isFromStrangeness = kFALSE;\r
2095\r
2096 if(stack) {\r
2097 TParticle* particle = stack->Particle(label);\r
2098\r
2099 if(particle && particle->GetPDG() && particle->GetPDG()->Charge()!=0) \r
2100 {\r
2101 Int_t mech = particle->GetUniqueID(); // production mechanism \r
2102 Bool_t isPrim = stack->IsPhysicalPrimary(label);\r
2103\r
2104 Int_t motherLabel = TMath::Abs(particle->GetMother(0)); \r
2105 TParticle* mother = stack->Particle(motherLabel); \r
2106 if(mother) {\r
2107 Int_t motherPdg = mother->GetPdgCode();\r
2108\r
2109 // K+-, lambda, antilambda, K0s decays\r
2110 if(!isPrim && mech==4 && \r
2111 (TMath::Abs(motherPdg)==kKPlus || TMath::Abs(motherPdg)==kLambda0 || motherPdg==kK0Short))\r
2112 {\r
2113 isFromStrangeness = kTRUE;\r
2114 } \r
2115 }\r
2116 } \r
2117 } \r
2118\r
2119 return isFromStrangeness;\r
2120}\r
2121\r
2122\r
2123//_____________________________________________________________________________\r
2124void AliAnalysisTaskFilteredTree::FinishTaskOutput() \r
2125{\r
2126 //\r
2127 // Called one at the end \r
2128 // locally on working node\r
2129 //\r
2130\r
a1419baa 2131 //// must be deleted to store trees\r
2132 //if(fTreeSRedirector) delete fTreeSRedirector; fTreeSRedirector=0;\r
2133\r
2134 //// open temporary file and copy trees to the ouptut container\r
2135\r
2136 //TChain* chain = 0;\r
2137 ////\r
2138 //chain = new TChain("highPt");\r
2139 //if(chain) { \r
2140 // chain->Add("jotwinow_Temp_Trees.root");\r
2141 // fHighPtTree = chain->CopyTree("1");\r
2142 // delete chain; chain=0; \r
2143 //}\r
2144 //if(fHighPtTree) fHighPtTree->Print();\r
2145\r
2146 ////\r
2147 //chain = new TChain("V0s");\r
2148 //if(chain) { \r
2149 // chain->Add("jotwinow_Temp_Trees.root");\r
2150 // fV0Tree = chain->CopyTree("1");\r
2151 // delete chain; chain=0; \r
2152 //}\r
2153 //if(fV0Tree) fV0Tree->Print();\r
2154\r
2155 ////\r
2156 //chain = new TChain("dEdx");\r
2157 //if(chain) { \r
2158 // chain->Add("jotwinow_Temp_Trees.root");\r
2159 // fdEdxTree = chain->CopyTree("1");\r
2160 // delete chain; chain=0; \r
2161 //}\r
2162 //if(fdEdxTree) fdEdxTree->Print();\r
2163\r
2164 ////\r
2165 //chain = new TChain("Laser");\r
2166 //if(chain) { \r
2167 // chain->Add("jotwinow_Temp_Trees.root");\r
2168 // fLaserTree = chain->CopyTree("1");\r
2169 // delete chain; chain=0; \r
2170 //}\r
2171 //if(fLaserTree) fLaserTree->Print();\r
2172\r
2173 ////\r
2174 //chain = new TChain("MCEffTree");\r
2175 //if(chain) { \r
2176 // chain->Add("jotwinow_Temp_Trees.root");\r
2177 // fMCEffTree = chain->CopyTree("1");\r
2178 // delete chain; chain=0; \r
2179 //}\r
2180 //if(fMCEffTree) fMCEffTree->Print();\r
2181\r
2182 ////\r
2183 //chain = new TChain("CosmicPairs");\r
2184 //if(chain) { \r
2185 // chain->Add("jotwinow_Temp_Trees.root");\r
2186 // fCosmicPairsTree = chain->CopyTree("1");\r
2187 // delete chain; chain=0; \r
2188 //}\r
2189 //if(fCosmicPairsTree) fCosmicPairsTree->Print(); \r
2190\r
2191\r
2192 //OpenFile(1);\r
1059f583 2193\r
2194 // Post output data.\r
a1419baa 2195 //PostData(1, fHighPtTree);\r
2196 //PostData(2, fV0Tree);\r
2197 //PostData(3, fdEdxTree);\r
2198 //PostData(4, fLaserTree);\r
2199 //PostData(5, fMCEffTree);\r
2200 //PostData(6, fCosmicPairsTree);\r
1059f583 2201}\r
2202\r
2203//_____________________________________________________________________________\r
2204void AliAnalysisTaskFilteredTree::Terminate(Option_t *) \r
2205{\r
2206 // Called one at the end \r
2207 /*\r
2208 fOutputSummary = dynamic_cast<TTree*> (GetOutputData(1));\r
2209 if(fOutputSummary) delete fOutputSummary; fOutputSummary=0;\r
2210 TChain* chain = new TChain("highPt");\r
2211 if(!chain) return;\r
2212 chain->Add("jotwinow_HighPt_TrackAndV0_Trees.root");\r
2213 TTree *tree = chain->CopyTree("1");\r
2214 if (chain) { delete chain; chain=0; }\r
2215 if(!tree) return;\r
2216 tree->Print();\r
2217 fOutputSummary = tree;\r
2218\r
2219 if (!fOutputSummary) {\r
2220 Printf("ERROR: AliAnalysisTaskFilteredTree::Terminate(): Output data not avaiable %p \n", GetOutputData(1));\r
2221 return;\r
2222 }\r
2223 */\r
2224}\r
2225\r
2226//_____________________________________________________________________________\r
2227Int_t AliAnalysisTaskFilteredTree::GetMCTrueTrackMult(AliMCEvent *const mcEvent, AliFilteredTreeEventCuts *const evtCuts, AliFilteredTreeAcceptanceCuts *const accCuts)\r
2228{\r
2229 //\r
2230 // calculate mc event true track multiplicity\r
2231 //\r
2232 if(!mcEvent) return 0;\r
2233\r
2234 AliStack* stack = 0;\r
2235 Int_t mult = 0;\r
2236\r
2237 // MC particle stack\r
2238 stack = mcEvent->Stack();\r
2239 if (!stack) return 0;\r
2240\r
2241 //\r
2242 //printf("minZv %f, maxZv %f \n", evtCuts->GetMinZv(), evtCuts->GetMaxZv());\r
2243 //\r
2244\r
2245 Bool_t isEventOK = evtCuts->AcceptMCEvent(mcEvent);\r
2246 if(!isEventOK) return 0; \r
2247\r
2248 Int_t nPart = stack->GetNtrack();\r
2249 for (Int_t iMc = 0; iMc < nPart; ++iMc) \r
2250 {\r
2251 TParticle* particle = stack->Particle(iMc);\r
2252 if (!particle)\r
2253 continue;\r
2254\r
2255 // only charged particles\r
2256 if(!particle->GetPDG()) continue;\r
2257 Double_t charge = particle->GetPDG()->Charge()/3.;\r
2258 if (TMath::Abs(charge) < 0.001)\r
2259 continue;\r
2260 \r
2261 // physical primary\r
2262 Bool_t prim = stack->IsPhysicalPrimary(iMc);\r
2263 if(!prim) continue;\r
2264\r
2265 // checked accepted without pt cut\r
2266 //if(accCuts->AcceptTrack(particle)) \r
2267 if( particle->Eta() > accCuts->GetMinEta() && particle->Eta() < accCuts->GetMaxEta() ) \r
2268 {\r
2269 mult++;\r
2270 }\r
2271 }\r
2272\r
2273return mult; \r
2274}\r
2275\r
2276\r