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