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