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