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