fix a bug where the tree remains empty for data
[u/mrichter/AliRoot.git] / PWGPP / AliAnalysisTaskFilteredTree.cxx
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 "AliVTrack.h"\r
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
63 using namespace std;\r
64 \r
65 ClassImp(AliAnalysisTaskFilteredTree)\r
66 \r
67 //_____________________________________________________________________________\r
68 AliAnalysisTaskFilteredTree::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
77   , fReducePileUp(kTRUE)\r
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
109 AliAnalysisTaskFilteredTree::~AliAnalysisTaskFilteredTree()\r
110 {\r
111   //the output trees not to be deleted in case of proof analysis\r
112   //Bool_t deleteTrees=kTRUE;\r
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
120 \r
121   delete fFilteredTreeEventCuts;\r
122   delete fFilteredTreeAcceptanceCuts;\r
123   delete fFilteredTreeRecAcceptanceCuts;\r
124   delete fEsdTrackCuts;\r
125 }\r
126 \r
127 //____________________________________________________________________________\r
128 Bool_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
138 return kTRUE;\r
139 }\r
140 \r
141 //_____________________________________________________________________________\r
142 void AliAnalysisTaskFilteredTree::UserCreateOutputObjects()\r
143 {\r
144   // Create histograms\r
145   // Called once\r
146 \r
147   //\r
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
151  \r
152   //\r
153   // Create trees\r
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
163   PostData(3,fdEdxTree);\r
164   PostData(4,fLaserTree);\r
165   PostData(5,fMCEffTree);\r
166   PostData(6,fCosmicPairsTree);\r
167 }\r
168 \r
169 //_____________________________________________________________________________\r
170 void AliAnalysisTaskFilteredTree::UserExec(Option_t *) \r
171 {\r
172   //\r
173   // Called for each event\r
174   //\r
175 \r
176   // ESD event\r
177   fESD = (AliESDEvent*) (InputEvent());\r
178   if (!fESD) {\r
179     Printf("ERROR: ESD event not available");\r
180     return;\r
181   }\r
182 \r
183   //// MC event\r
184   //if(fUseMCInfo) {\r
185   //  fMC = MCEvent();\r
186   //  if (!fMC) {\r
187   //    Printf("ERROR: MC event not available");\r
188   //    return;\r
189   //  }\r
190   //}\r
191   \r
192   //if MC info available - use it.\r
193   fMC = MCEvent();\r
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
216   if(fMC) { ProcessMCEff(fESD,fMC,fESDfriend);}\r
217 }\r
218 \r
219 //_____________________________________________________________________________\r
220 void 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
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
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
362             "triggerClass="<<&triggerClass<<      //  trigger\r
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
378 void 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
449   if(mcEvent)\r
450   {\r
451     // get MC event header\r
452     header = mcEvent->Header();\r
453     if (!header) {\r
454       AliDebug(AliLog::kError, "Header not available");\r
455       return;\r
456     }\r
457     // MC particle stack\r
458     stack = mcEvent->Stack();\r
459     if (!stack) {\r
460       AliDebug(AliLog::kError, "Stack not available");\r
461       return;\r
462     }\r
463 \r
464     // get MC vertex\r
465     genHeader = header->GenEventHeader();\r
466     if (!genHeader) {\r
467       AliDebug(AliLog::kError, "Could not retrieve genHeader from Header");\r
468       return;\r
469     }\r
470     genHeader->PrimaryVertex(vtxMC);\r
471 \r
472     // multipliticy of all MC primary tracks\r
473     // in Zv, pt and eta ranges)\r
474     multMCTrueTracks = GetMCTrueTrackMult(mcEvent,evtCuts,accCuts);\r
475   } // 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
495   Int_t ntracks = esdEvent->GetNumberOfTracks();\r
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
520     Float_t bz = esdEvent->GetMagneticField();\r
521     Int_t runNumber = esdEvent->GetRunNumber();\r
522     Int_t evtTimeStamp = esdEvent->GetTimeStamp();\r
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
553       TObjString triggerClass = esdEvent->GetFiredTriggerClasses().Data();\r
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
573       if(!fTreeSRedirector) return;\r
574       (*fTreeSRedirector)<<"highPt"<<\r
575         "fileName.="<<&fileName<<            \r
576         "runNumber="<<runNumber<<\r
577         "evtTimeStamp="<<evtTimeStamp<<\r
578         "evtNumberInFile="<<evtNumberInFile<<\r
579         "triggerClass="<<&triggerClass<<      //  trigger\r
580         "Bz="<<bz<<                           //  magnetic field\r
581         "vtxESD.="<<vtxESD<<\r
582         "ntracksESD="<<ntracks<<              // number of tracks in the ESD\r
583       //  "vtxESDx="<<vtxX<<\r
584       //  "vtxESDy="<<vtxY<<\r
585       //  "vtxESDz="<<vtxZ<<\r
586         "IRtot="<<ir1<<                      // interaction record history info\r
587         "IRint2="<<ir2<<\r
588         "mult="<<mult<<                      // multiplicity of tracks pointing to the primary vertex\r
589         "esdTrack.="<<track<<\r
590         "centralityF="<<centralityF<<       \r
591         "\n";\r
592     }\r
593   }\r
594   \r
595 }\r
596 \r
597 \r
598 //_____________________________________________________________________________\r
599 void 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
632       Int_t runNumber = esdEvent->GetRunNumber();\r
633       Int_t evtTimeStamp = esdEvent->GetTimeStamp();\r
634       Int_t evtNumberInFile = esdEvent->GetEventNumberInFile();\r
635       Float_t bz = esdEvent->GetMagneticField();\r
636       TObjString triggerClass = esdEvent->GetFiredTriggerClasses().Data();\r
637 \r
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
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
654         "triggerClass="<<&triggerClass<<      //  trigger\r
655         "Bz="<<bz<<\r
656         "multTPCtracks="<<countLaserTracks<<\r
657         "\n";\r
658     }\r
659   }\r
660 }\r
661 \r
662 //_____________________________________________________________________________\r
663 void AliAnalysisTaskFilteredTree::ProcessAll(AliESDEvent *const esdEvent, AliMCEvent * const mcEvent, AliESDfriend *const esdFriend)\r
664 {\r
665   //\r
666   // Select real events with high-pT tracks\r
667   // Calculate and stor in the output tree:\r
668   //  TPC constrained tracks\r
669   //  InnerParams constrained tracks\r
670   //  TPC-ITS tracks\r
671   //  ITSout-InnerParams tracks\r
672   //  chi2 distance between TPC constrained and TPC-ITS tracks\r
673   //  chi2 distance between TPC refitted constrained and TPC-ITS tracks\r
674   //  chi2 distance between ITSout and InnerParams tracks\r
675   //  MC information: \r
676   //   track references at ITSin, TPCin; InnerParam at first TPC track reference, \r
677   //   particle ID, mother ID, production mechanism ...\r
678   // \r
679 \r
680   if(!esdEvent) {\r
681     AliDebug(AliLog::kError, "esdEvent not available");\r
682     return;\r
683   }\r
684 \r
685   // get selection cuts\r
686   AliFilteredTreeEventCuts *evtCuts = GetEventCuts(); \r
687   AliFilteredTreeAcceptanceCuts *accCuts = GetAcceptanceCuts(); \r
688   AliESDtrackCuts *esdTrackCuts = GetTrackCuts(); \r
689 \r
690   if(!evtCuts || !accCuts  || !esdTrackCuts) {\r
691     AliDebug(AliLog::kError, "cuts not available");\r
692     return;\r
693   }\r
694 \r
695   // trigger selection\r
696   Bool_t isEventTriggered = kTRUE;\r
697   AliPhysicsSelection *physicsSelection = NULL;\r
698   AliTriggerAnalysis* triggerAnalysis = NULL;\r
699 \r
700   // \r
701   AliInputEventHandler* inputHandler = (AliInputEventHandler*) AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler();\r
702   if (!inputHandler)\r
703   {\r
704     Printf("ERROR: Could not receive input handler");\r
705     return;\r
706   }\r
707    \r
708   // get file name\r
709   TTree *chain = (TChain*)GetInputData(0);\r
710   if(!chain) { \r
711     Printf("ERROR: Could not receive input chain");\r
712     return;\r
713   }\r
714   TObjString fileName(chain->GetCurrentFile()->GetName());\r
715 \r
716   // trigger\r
717   if(evtCuts->IsTriggerRequired())  \r
718   {\r
719     // always MB\r
720     isEventTriggered = inputHandler->IsEventSelected() & AliVEvent::kMB;\r
721 \r
722     physicsSelection = static_cast<AliPhysicsSelection*> (inputHandler->GetEventSelection());\r
723     if(!physicsSelection) {AliInfo("no physics selection"); return;}\r
724     //SetPhysicsTriggerSelection(physicsSelection);\r
725 \r
726     if (isEventTriggered && (GetTrigger() == AliTriggerAnalysis::kV0AND)) {\r
727       // set trigger (V0AND)\r
728       triggerAnalysis = physicsSelection->GetTriggerAnalysis();\r
729       if(!triggerAnalysis) {AliInfo("no trigger analysis");return;}\r
730       isEventTriggered = triggerAnalysis->IsOfflineTriggerFired(esdEvent, GetTrigger());\r
731     }\r
732   }\r
733 \r
734   // centrality determination\r
735   Float_t centralityF = -1;\r
736   AliCentrality *esdCentrality = esdEvent->GetCentrality();\r
737   centralityF = esdCentrality->GetCentralityPercentile(fCentralityEstimator.Data());\r
738 \r
739   // use MC information\r
740   AliHeader* header = 0;\r
741   AliGenEventHeader* genHeader = 0;\r
742   AliStack* stack = 0;\r
743   TArrayF vtxMC(3);\r
744 \r
745   Int_t multMCTrueTracks = 0;\r
746   if(mcEvent)\r
747   {\r
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
785     AliInfo("no ESD vertex");\r
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
804     if(!esdHeader) {AliInfo("no esdHeader");return;}\r
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
817     Float_t bz = esdEvent->GetMagneticField();\r
818     Int_t runNumber = esdEvent->GetRunNumber();\r
819     Int_t evtTimeStamp = esdEvent->GetTimeStamp();\r
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
1010       AliTrackReference *refITS = NULL;\r
1011 \r
1012       Bool_t isOKtrackInnerC3 = kFALSE;\r
1013       AliExternalTrackParam *trackInnerC3 = new AliExternalTrackParam(*(track->GetInnerParam()));\r
1014 \r
1015       if(mcEvent) \r
1016       {\r
1017         if(!stack) return;\r
1018 \r
1019         //\r
1020         // global track\r
1021         //\r
1022         Int_t label = TMath::Abs(track->GetLabel()); \r
1023         particle = stack->Particle(label);\r
1024         if(particle && particle->GetPDG() && particle->GetPDG()->Charge()!=0.)\r
1025         {\r
1026           particleMother = GetMother(particle,stack);\r
1027           mech = particle->GetUniqueID();\r
1028           isPrim = stack->IsPhysicalPrimary(label);\r
1029           isFromStrangess  = IsFromStrangeness(label,stack);\r
1030           isFromConversion = IsFromConversion(label,stack);\r
1031           isFromMaterial   = IsFromMaterial(label,stack);\r
1032         }\r
1033 \r
1034         //\r
1035         // TPC track\r
1036         //\r
1037         Int_t labelTPC = TMath::Abs(track->GetTPCLabel()); \r
1038         particleTPC = stack->Particle(labelTPC);\r
1039         if(particleTPC && particleTPC->GetPDG() && particleTPC->GetPDG()->Charge()!=0.)\r
1040         {\r
1041           particleMotherTPC = GetMother(particleTPC,stack);\r
1042           mechTPC = particleTPC->GetUniqueID();\r
1043           isPrimTPC = stack->IsPhysicalPrimary(labelTPC);\r
1044           isFromStrangessTPC  = IsFromStrangeness(labelTPC,stack);\r
1045           isFromConversionTPC = IsFromConversion(labelTPC,stack);\r
1046           isFromMaterialTPC   = IsFromMaterial(labelTPC,stack);\r
1047         }\r
1048 \r
1049         //\r
1050         // store first track reference\r
1051         // for TPC track\r
1052         //\r
1053         TParticle *part=0;\r
1054         TClonesArray *trefs=0;\r
1055         Int_t status = mcEvent->GetParticleAndTR(track->GetTPCLabel(), part, trefs);\r
1056 \r
1057         if(status>0 && part && trefs && part->GetPDG() && part->GetPDG()->Charge()!=0.) \r
1058         {\r
1059           Int_t nTrackRef = trefs->GetEntries();\r
1060           //printf("nTrackRef %d \n",nTrackRef);\r
1061 \r
1062           Int_t countITS = 0;\r
1063           for (Int_t iref = 0; iref < nTrackRef; iref++) \r
1064           {\r
1065             AliTrackReference *ref = (AliTrackReference *)trefs->At(iref);\r
1066 \r
1067              // Ref. in the middle ITS \r
1068             if(ref && ref->DetectorId()==AliTrackReference::kITS)\r
1069             {\r
1070               if(!refITS && countITS==2) {\r
1071                 refITS = ref;\r
1072                 //printf("refITS %p \n",refITS);\r
1073               }\r
1074               countITS++;\r
1075             }\r
1076 \r
1077             // TPC\r
1078             if(ref && ref->DetectorId()==AliTrackReference::kTPC)\r
1079             {\r
1080               if(!refTPCIn) {\r
1081                 refTPCIn = ref;\r
1082                 //printf("refTPCIn %p \n",refTPCIn);\r
1083                 //break;\r
1084               }\r
1085             }\r
1086           }\r
1087 \r
1088           // transform inner params to TrackRef\r
1089           // reference frame\r
1090           if(refTPCIn && trackInnerC3) \r
1091           {\r
1092             Double_t kRefPhi = TMath::ATan2(refTPCIn->Y(),refTPCIn->X());\r
1093             isOKtrackInnerC3 = trackInnerC3->Rotate(kRefPhi);\r
1094             isOKtrackInnerC3 = AliTracker::PropagateTrackToBxByBz(trackInnerC3,refTPCIn->R(),track->GetMass(),kStep,kFALSE);\r
1095           }\r
1096         }\r
1097 \r
1098         //\r
1099         // ITS track\r
1100         //\r
1101         Int_t labelITS = TMath::Abs(track->GetITSLabel()); \r
1102         particleITS = stack->Particle(labelITS);\r
1103         if(particleITS && particleITS->GetPDG() && particleITS->GetPDG()->Charge()!=0.)\r
1104         {\r
1105           particleMotherITS = GetMother(particleITS,stack);\r
1106           mechITS = particleITS->GetUniqueID();\r
1107           isPrimITS = stack->IsPhysicalPrimary(labelITS);\r
1108           isFromStrangessITS  = IsFromStrangeness(labelITS,stack);\r
1109           isFromConversionITS = IsFromConversion(labelITS,stack);\r
1110           isFromMaterialITS   = IsFromMaterial(labelITS,stack);\r
1111         }\r
1112       }\r
1113 \r
1114       //\r
1115       Bool_t dumpToTree=kFALSE;\r
1116       \r
1117       if(isOKtpcInnerC  && isOKtrackInnerC) dumpToTree = kTRUE;\r
1118       if(fUseESDfriends && isOKtrackInnerC2 && isOKouterITSc) dumpToTree = kTRUE;\r
1119       if(mcEvent && isOKtrackInnerC3) dumpToTree = kTRUE;\r
1120       TObjString triggerClass = esdEvent->GetFiredTriggerClasses().Data();\r
1121       if (fReducePileUp){  \r
1122         //\r
1123         // 18.03 - Reduce pile-up chunks, done outside of the ESDTrackCuts for 2012/2013 data pile-up about 95 % of tracks\r
1124         // Done only in case no MC info \r
1125         //\r
1126         Float_t dcaTPC[2];\r
1127         track->GetImpactParametersTPC(dcaTPC[0],dcaTPC[1]);\r
1128         Bool_t isRoughPrimary = TMath::Abs(dcaTPC[1])<10;\r
1129         Bool_t hasOuter=(track->IsOn(AliVTrack::kITSin))||(track->IsOn(AliVTrack::kTOFout))||(track->IsOn(AliVTrack::kTRDin));\r
1130         Bool_t keepPileUp=gRandom->Rndm()<0.05;\r
1131         if ( (!hasOuter) && (!isRoughPrimary) && (!keepPileUp)){\r
1132           dumpToTree=kFALSE;\r
1133         }\r
1134       }\r
1135       /////////////////\r
1136       //book keeping of created dummy objects (to avoid NULL in trees)\r
1137       Bool_t newvtxESD=kFALSE;\r
1138       Bool_t newtrack=kFALSE;\r
1139       Bool_t newtpcInnerC=kFALSE;\r
1140       Bool_t newtrackInnerC=kFALSE;\r
1141       Bool_t newtrackInnerC2=kFALSE;\r
1142       Bool_t newouterITSc=kFALSE;\r
1143       Bool_t newtrackInnerC3=kFALSE;\r
1144       Bool_t newrefTPCIn=kFALSE;\r
1145       Bool_t newrefITS=kFALSE;\r
1146       Bool_t newparticle=kFALSE;\r
1147       Bool_t newparticleMother=kFALSE;\r
1148       Bool_t newparticleTPC=kFALSE;\r
1149       Bool_t newparticleMotherTPC=kFALSE;\r
1150       Bool_t newparticleITS=kFALSE;\r
1151       Bool_t newparticleMotherITS=kFALSE;\r
1152       \r
1153       //check that the vertex is there and that it is OK, \r
1154       //i.e. no null member arrays, otherwise a problem with merging\r
1155       //later on.\r
1156       //this is a very ugly hack!\r
1157       if (!vtxESD)\r
1158       {\r
1159         AliInfo("fixing the ESD vertex for streaming");\r
1160         vtxESD=new AliESDVertex(); \r
1161         //vtxESD->SetNContributors(1);\r
1162         //UShort_t pindices[1]; pindices[0]=0;\r
1163         //vtxESD->SetIndices(1,pindices);\r
1164         //vtxESD->SetNContributors(0);\r
1165         newvtxESD=kTRUE;\r
1166       }\r
1167       //\r
1168       if (!track) {track=new AliESDtrack();newtrack=kTRUE;}\r
1169       if (!tpcInnerC) {tpcInnerC=new AliExternalTrackParam();newtpcInnerC=kTRUE;}\r
1170       if (!trackInnerC) {trackInnerC=new AliExternalTrackParam();newtrackInnerC=kTRUE;}\r
1171       if (!trackInnerC2) {trackInnerC2=new AliExternalTrackParam();newtrackInnerC2=kTRUE;}\r
1172       if (!outerITSc) {outerITSc=new AliExternalTrackParam();newouterITSc=kTRUE;}\r
1173       if (!trackInnerC3) {trackInnerC3=new AliExternalTrackParam();newtrackInnerC3=kTRUE;}\r
1174       if (mcEvent)\r
1175       {\r
1176         if (!refTPCIn) {refTPCIn=new AliTrackReference(); newrefTPCIn=kTRUE;}\r
1177         if (!refITS) {refITS=new AliTrackReference();newrefITS=kTRUE;}\r
1178         if (!particle) {particle=new TParticle(); newparticle=kTRUE;}\r
1179         if (!particleMother) {particleMother=new TParticle();newparticleMother=kTRUE;}\r
1180         if (!particleTPC) {particleTPC=new TParticle();newparticleTPC=kTRUE;}\r
1181         if (!particleMotherTPC) {particleMotherTPC=new TParticle();newparticleMotherTPC=kTRUE;}\r
1182         if (!particleITS) {particleITS=new TParticle();newparticleITS=kTRUE;}\r
1183         if (!particleMotherITS) {particleMotherITS=new TParticle();newparticleMotherITS=kTRUE;}\r
1184       }\r
1185       /////////////////////////\r
1186 \r
1187       //Double_t vtxX=vtxESD->GetX();\r
1188       //Double_t vtxY=vtxESD->GetY();\r
1189       //Double_t vtxZ=vtxESD->GetZ();\r
1190 \r
1191       AliESDVertex* pvtxESD = (AliESDVertex*)vtxESD->Clone();\r
1192       AliESDtrack* ptrack=(AliESDtrack*)track->Clone();\r
1193       AliExternalTrackParam* ptpcInnerC = (AliExternalTrackParam*)tpcInnerC->Clone();\r
1194       AliExternalTrackParam* ptrackInnerC = (AliExternalTrackParam*)trackInnerC->Clone();\r
1195       AliExternalTrackParam* ptrackInnerC2 = (AliExternalTrackParam*)trackInnerC2->Clone();\r
1196       AliExternalTrackParam* pouterITSc = (AliExternalTrackParam*)outerITSc->Clone();\r
1197       AliExternalTrackParam* ptrackInnerC3 = (AliExternalTrackParam*)trackInnerC3->Clone();\r
1198       Int_t ntracks = esdEvent->GetNumberOfTracks();\r
1199       \r
1200       if(fTreeSRedirector && dumpToTree) \r
1201       {\r
1202 \r
1203         (*fTreeSRedirector)<<"highPt"<<\r
1204           "fileName.="<<&fileName<<                // name of the chunk file (hopefully full)\r
1205           "runNumber="<<runNumber<<                // runNumber\r
1206           "evtTimeStamp="<<evtTimeStamp<<          // time stamp of event (in seconds)\r
1207           "evtNumberInFile="<<evtNumberInFile<<    // event number\r
1208           "triggerClass="<<&triggerClass<<         // trigger class as a string\r
1209           "Bz="<<bz<<                              // solenoid magnetic field in the z direction (in kGaus)\r
1210           "vtxESD.="<<pvtxESD<<                    // vertexer ESD tracks (can be biased by TPC pileup tracks)\r
1211           //"vtxESDx="<<vtxX<<\r
1212           //"vtxESDy="<<vtxY<<\r
1213           //"vtxESDz="<<vtxZ<<\r
1214           "IRtot="<<ir1<<                         // interaction record (trigger) counters - coutner 1\r
1215           "IRint2="<<ir2<<                        // interaction record (trigger) coutners - counter 2\r
1216           "mult="<<mult<<                         // multiplicity of tracks pointing to the primary vertex\r
1217           "ntracks="<<ntracks<<                   // number of the esd tracks (to take into account the pileup in the TPC)\r
1218           "esdTrack.="<<ptrack<<                  // esdTrack as used in the physical analysis\r
1219           "extTPCInnerC.="<<ptpcInnerC<<          // ??? \r
1220           "extInnerParamC.="<<ptrackInnerC<<      // ???\r
1221           "extInnerParam.="<<ptrackInnerC2<<      // ???\r
1222           "extOuterITS.="<<pouterITSc<<           // ???\r
1223           "extInnerParamRef.="<<ptrackInnerC3<<   // ???\r
1224           "chi2TPCInnerC="<<chi2(0,0)<<           // chi2   of tracks ???\r
1225           "chi2InnerC="<<chi2trackC(0,0)<<        // chi2s  of tracks TPCinner to the combined\r
1226           "chi2OuterITS="<<chi2OuterITS(0,0)<<    // chi2s  of tracks TPC at inner wall to the ITSout\r
1227           "centralityF="<<centralityF;\r
1228         if (mcEvent)\r
1229         {\r
1230           (*fTreeSRedirector)<<"highPt"<<\r
1231           "refTPCIn.="<<refTPCIn<<\r
1232           "refITS.="<<refITS<<\r
1233           "particle.="<<particle<<\r
1234           "particleMother.="<<particleMother<<\r
1235           "mech="<<mech<<\r
1236           "isPrim="<<isPrim<<\r
1237           "isFromStrangess="<<isFromStrangess<<\r
1238           "isFromConversion="<<isFromConversion<<\r
1239           "isFromMaterial="<<isFromMaterial<<\r
1240           "particleTPC.="<<particleTPC<<\r
1241           "particleMotherTPC.="<<particleMotherTPC<<\r
1242           "mechTPC="<<mechTPC<<\r
1243           "isPrimTPC="<<isPrimTPC<<\r
1244           "isFromStrangessTPC="<<isFromStrangessTPC<<\r
1245           "isFromConversionTPC="<<isFromConversionTPC<<\r
1246           "isFromMaterialTPC="<<isFromMaterialTPC<<\r
1247           "particleITS.="<<particleITS<<\r
1248           "particleMotherITS.="<<particleMotherITS<<\r
1249           "mechITS="<<mechITS<<\r
1250           "isPrimITS="<<isPrimITS<<\r
1251           "isFromStrangessITS="<<isFromStrangessITS<<\r
1252           "isFromConversionITS="<<isFromConversionITS<<\r
1253           "isFromMaterialITS="<<isFromMaterialITS;\r
1254         }\r
1255         //finish writing the entry\r
1256         AliInfo("writing tree highPt");\r
1257         (*fTreeSRedirector)<<"highPt"<<"\n";\r
1258       }\r
1259 \r
1260       delete pvtxESD;\r
1261       delete ptrack;\r
1262       delete ptpcInnerC;\r
1263       delete ptrackInnerC;\r
1264       delete ptrackInnerC2;\r
1265       delete pouterITSc;\r
1266       delete ptrackInnerC3;\r
1267 \r
1268       ////////////////////\r
1269       //delete the dummy objects we might have created.\r
1270       if (newvtxESD) {delete vtxESD; vtxESD=NULL;}\r
1271       if (newtrack) {delete track; track=NULL;}\r
1272       if (newtpcInnerC) {delete tpcInnerC; tpcInnerC=NULL;}\r
1273       if (newtrackInnerC) {delete trackInnerC; trackInnerC=NULL;}\r
1274       if (newtrackInnerC2) {delete trackInnerC2; trackInnerC2=NULL;}\r
1275       if (newouterITSc) {delete outerITSc; outerITSc=NULL;}\r
1276       if (newtrackInnerC3) {delete trackInnerC3; trackInnerC3=NULL;}\r
1277       if (newrefTPCIn) {delete refTPCIn; refTPCIn=NULL;}\r
1278       if (newrefITS) {delete refITS; refITS=NULL;}\r
1279       if (newparticle) {delete particle; particle=NULL;}\r
1280       if (newparticleMother) {delete particleMother; particleMother=NULL;}\r
1281       if (newparticleTPC) {delete particleTPC; particleTPC=NULL;}\r
1282       if (newparticleMotherTPC) {delete particleMotherTPC; particleMotherTPC=NULL;}\r
1283       if (newparticleITS) {delete particleITS; particleITS=NULL;}\r
1284       if (newparticleMotherITS) {delete particleMotherITS; particleMotherITS=NULL;}\r
1285       ///////////////\r
1286 \r
1287       delete tpcInnerC;\r
1288       delete trackInnerC;\r
1289       delete trackInnerC2;\r
1290       delete outerITSc;\r
1291       delete trackInnerC3;\r
1292     }\r
1293   }\r
1294 \r
1295 }\r
1296 \r
1297 \r
1298 //_____________________________________________________________________________\r
1299 void AliAnalysisTaskFilteredTree::ProcessMCEff(AliESDEvent *const esdEvent, AliMCEvent * const mcEvent, AliESDfriend *const /*esdFriend*/)\r
1300 {\r
1301   //\r
1302   // Fill tree for efficiency studies MC only\r
1303   AliInfo("we start!");\r
1304 \r
1305   if(!esdEvent) {\r
1306     AliDebug(AliLog::kError, "esdEvent not available");\r
1307     return;\r
1308   }\r
1309 \r
1310   if(!mcEvent) {\r
1311     AliDebug(AliLog::kError, "mcEvent not available");\r
1312     return;\r
1313   }\r
1314 \r
1315   // get selection cuts\r
1316   AliFilteredTreeEventCuts *evtCuts = GetEventCuts(); \r
1317   AliFilteredTreeAcceptanceCuts *accCuts = GetAcceptanceCuts(); \r
1318   AliESDtrackCuts *esdTrackCuts = GetTrackCuts(); \r
1319 \r
1320   if(!evtCuts || !accCuts  || !esdTrackCuts) {\r
1321     AliDebug(AliLog::kError, "cuts not available");\r
1322     return;\r
1323   }\r
1324 \r
1325   // trigger selection\r
1326   Bool_t isEventTriggered = kTRUE;\r
1327   AliPhysicsSelection *physicsSelection = NULL;\r
1328   AliTriggerAnalysis* triggerAnalysis = NULL;\r
1329 \r
1330   // \r
1331   AliInputEventHandler* inputHandler = (AliInputEventHandler*) AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler();\r
1332   if (!inputHandler)\r
1333   {\r
1334     Printf("ERROR: Could not receive input handler");\r
1335     return;\r
1336   }\r
1337 \r
1338   // get file name\r
1339   TTree *chain = (TChain*)GetInputData(0);\r
1340   if(!chain) { \r
1341     Printf("ERROR: Could not receive input chain");\r
1342     return;\r
1343   }\r
1344   TObjString fileName(chain->GetCurrentFile()->GetName());\r
1345 \r
1346   // trigger\r
1347   if(evtCuts->IsTriggerRequired())  \r
1348   {\r
1349     // always MB\r
1350     isEventTriggered = inputHandler->IsEventSelected() & AliVEvent::kMB;\r
1351 \r
1352     physicsSelection = static_cast<AliPhysicsSelection*> (inputHandler->GetEventSelection());\r
1353     if(!physicsSelection) return;\r
1354 \r
1355     if (isEventTriggered && (GetTrigger() == AliTriggerAnalysis::kV0AND)) {\r
1356       // set trigger (V0AND)\r
1357       triggerAnalysis = physicsSelection->GetTriggerAnalysis();\r
1358       if(!triggerAnalysis) return;\r
1359       isEventTriggered = triggerAnalysis->IsOfflineTriggerFired(esdEvent, GetTrigger());\r
1360     }\r
1361   }\r
1362 \r
1363   // centrality determination\r
1364   Float_t centralityF = -1;\r
1365   AliCentrality *esdCentrality = esdEvent->GetCentrality();\r
1366   centralityF = esdCentrality->GetCentralityPercentile(fCentralityEstimator.Data());\r
1367 \r
1368   // use MC information\r
1369   AliHeader* header = 0;\r
1370   AliGenEventHeader* genHeader = 0;\r
1371   AliStack* stack = 0;\r
1372   TArrayF vtxMC(3);\r
1373 \r
1374   Int_t multMCTrueTracks = 0;\r
1375   //\r
1376   if(!mcEvent) {\r
1377     AliDebug(AliLog::kError, "mcEvent not available");\r
1378     return;\r
1379   }\r
1380   // get MC event header\r
1381   header = mcEvent->Header();\r
1382   if (!header) {\r
1383     AliDebug(AliLog::kError, "Header not available");\r
1384     return;\r
1385   }\r
1386   // MC particle stack\r
1387   stack = mcEvent->Stack();\r
1388   if (!stack) {\r
1389     AliDebug(AliLog::kError, "Stack not available");\r
1390     return;\r
1391   }\r
1392 \r
1393   // get MC vertex\r
1394   genHeader = header->GenEventHeader();\r
1395   if (!genHeader) {\r
1396     AliDebug(AliLog::kError, "Could not retrieve genHeader from Header");\r
1397     return;\r
1398   }\r
1399   genHeader->PrimaryVertex(vtxMC);\r
1400 \r
1401   // multipliticy of all MC primary tracks\r
1402   // in Zv, pt and eta ranges)\r
1403   multMCTrueTracks = GetMCTrueTrackMult(mcEvent,evtCuts,accCuts);\r
1404 \r
1405 \r
1406   // get reconstructed vertex  \r
1407   //const AliESDVertex* vtxESD = 0; \r
1408   AliESDVertex* vtxESD = 0; \r
1409   if(GetAnalysisMode() == kTPCAnalysisMode) {\r
1410     vtxESD = (AliESDVertex*)esdEvent->GetPrimaryVertexTPC();\r
1411   }\r
1412   else if(GetAnalysisMode() == kTPCITSAnalysisMode) {\r
1413     vtxESD = (AliESDVertex*)esdEvent->GetPrimaryVertexTracks();\r
1414   }\r
1415   else {\r
1416     return;\r
1417   }\r
1418 \r
1419   if(!vtxESD) return;\r
1420 \r
1421   Bool_t isEventOK = evtCuts->AcceptEvent(esdEvent,mcEvent,vtxESD); \r
1422   //printf("isEventOK %d, isEventTriggered %d \n",isEventOK, isEventTriggered);\r
1423   //printf("GetAnalysisMode() %d \n",GetAnalysisMode());\r
1424 \r
1425   TObjString triggerClass = esdEvent->GetFiredTriggerClasses().Data();\r
1426 \r
1427   // check event cuts\r
1428   if(isEventOK && isEventTriggered)\r
1429   {\r
1430     if(!stack) return;\r
1431 \r
1432     //\r
1433     // MC info\r
1434     //\r
1435     TParticle *particle=NULL;\r
1436     TParticle *particleMother=NULL;\r
1437     Int_t mech=-1;\r
1438 \r
1439     // reco event info\r
1440     Double_t vert[3] = {0}; \r
1441     vert[0] = vtxESD->GetXv();\r
1442     vert[1] = vtxESD->GetYv();\r
1443     vert[2] = vtxESD->GetZv();\r
1444     Int_t mult = vtxESD->GetNContributors();\r
1445     Double_t bz = esdEvent->GetMagneticField();\r
1446     Double_t runNumber = esdEvent->GetRunNumber();\r
1447     Double_t evtTimeStamp = esdEvent->GetTimeStamp();\r
1448     Int_t evtNumberInFile = esdEvent->GetEventNumberInFile();\r
1449 \r
1450     // loop over MC stack\r
1451     for (Int_t iMc = 0; iMc < stack->GetNtrack(); ++iMc) \r
1452     {\r
1453       particle = stack->Particle(iMc);\r
1454       if (!particle)\r
1455         continue;\r
1456 \r
1457       // only charged particles\r
1458       if(!particle->GetPDG()) continue;\r
1459       Double_t charge = particle->GetPDG()->Charge()/3.;\r
1460       if (TMath::Abs(charge) < 0.001)\r
1461         continue;\r
1462 \r
1463       // only primary particles\r
1464       Bool_t prim = stack->IsPhysicalPrimary(iMc);\r
1465       if(!prim) continue;\r
1466 \r
1467       // downscale low-pT particles\r
1468       Double_t scalempt= TMath::Min(particle->Pt(),10.);\r
1469       Double_t downscaleF = gRandom->Rndm();\r
1470       downscaleF *= fLowPtTrackDownscaligF;\r
1471       if(TMath::Exp(2*scalempt)<downscaleF) continue;\r
1472 \r
1473       // is particle in acceptance\r
1474       if(!accCuts->AcceptTrack(particle)) continue;\r
1475 \r
1476       // check if particle reconstructed\r
1477       Bool_t isRec = kFALSE;\r
1478       Int_t  trackIndex = -1;\r
1479       for (Int_t iTrack = 0; iTrack < esdEvent->GetNumberOfTracks(); iTrack++)\r
1480       {\r
1481 \r
1482         AliESDtrack *track = esdEvent->GetTrack(iTrack);\r
1483         if(!track) continue;\r
1484         if(track->Charge()==0) continue;\r
1485         if(esdTrackCuts->AcceptTrack(track) && accCuts->AcceptTrack(track)) \r
1486         {\r
1487           Int_t label =  TMath::Abs(track->GetLabel());\r
1488           if(label == iMc) {\r
1489             isRec = kTRUE;\r
1490             trackIndex = iTrack;\r
1491             break;\r
1492           }\r
1493         } \r
1494       }\r
1495 \r
1496       // Store information in the output tree\r
1497       AliESDtrack *recTrack = NULL; \r
1498       if(trackIndex>-1)  { \r
1499         recTrack = esdEvent->GetTrack(trackIndex); \r
1500       } else {\r
1501         recTrack = new AliESDtrack(); \r
1502       } \r
1503 \r
1504       particleMother = GetMother(particle,stack);\r
1505       mech = particle->GetUniqueID();\r
1506 \r
1507       //MC particle track length\r
1508       Double_t tpcTrackLength = 0.;\r
1509       AliMCParticle *mcParticle = (AliMCParticle*) mcEvent->GetTrack(iMc);\r
1510       if(mcParticle) {\r
1511         Int_t counter;\r
1512         tpcTrackLength = mcParticle->GetTPCTrackLength(bz,0.05,counter,3.0);\r
1513       } \r
1514 \r
1515 \r
1516       //\r
1517       if(fTreeSRedirector) {\r
1518         (*fTreeSRedirector)<<"MCEffTree"<<\r
1519           "fileName.="<<&fileName<<\r
1520           "triggerClass.="<<&triggerClass<<\r
1521           "runNumber="<<runNumber<<\r
1522           "evtTimeStamp="<<evtTimeStamp<<\r
1523           "evtNumberInFile="<<evtNumberInFile<<\r
1524           "Bz="<<bz<<\r
1525           "vtxESD.="<<vtxESD<<\r
1526           "mult="<<mult<<\r
1527           "esdTrack.="<<recTrack<<\r
1528           "isRec="<<isRec<<\r
1529           "tpcTrackLength="<<tpcTrackLength<<\r
1530           "particle.="<<particle<<\r
1531           "particleMother.="<<particleMother<<\r
1532           "mech="<<mech<<\r
1533           "\n";\r
1534       }\r
1535 \r
1536       if(trackIndex <0 && recTrack) delete recTrack; recTrack=0;\r
1537     }\r
1538   }\r
1539 \r
1540 }\r
1541 \r
1542 //_____________________________________________________________________________\r
1543 Bool_t AliAnalysisTaskFilteredTree::IsHighDeDxParticle(AliESDtrack * track) {\r
1544   //\r
1545   // check if particle is Z > 1 \r
1546   //\r
1547   if (track->GetTPCNcls() < 60) return kFALSE;\r
1548   Double_t mom = track->GetInnerParam()->GetP();\r
1549   if (mom < 0.2) return kFALSE; // protection against unexpected behavior of Aleph parameterization\r
1550   Float_t dca[2], bCov[3];\r
1551   track->GetImpactParameters(dca,bCov);\r
1552   //\r
1553 \r
1554   Double_t triggerDeDx = 4*AliExternalTrackParam::BetheBlochAleph((mom*2)/(0.938*3),1.0288,31.9806,5.04114e-11,2.13096,2.38541);\r
1555 \r
1556   if (track->GetTPCsignal() > triggerDeDx && track->GetTPCsignal()<1000 && TMath::Abs(dca[0])<3.) return kTRUE;\r
1557 \r
1558   return kFALSE;\r
1559 }\r
1560 \r
1561 //_____________________________________________________________________________\r
1562 void AliAnalysisTaskFilteredTree::ProcessV0(AliESDEvent *const esdEvent, AliMCEvent * const mcEvent, AliESDfriend *const /*esdFriend*/)\r
1563 {\r
1564   //\r
1565   // Select real events with V0 (K0s and Lambda and Gamma) high-pT candidates\r
1566   //\r
1567   if(!esdEvent) {\r
1568     AliDebug(AliLog::kError, "esdEvent not available");\r
1569     return;\r
1570   }\r
1571 \r
1572   // get selection cuts\r
1573   AliFilteredTreeEventCuts *evtCuts = GetEventCuts(); \r
1574   AliFilteredTreeAcceptanceCuts *accCuts = GetAcceptanceCuts(); \r
1575   AliESDtrackCuts *esdTrackCuts = GetTrackCuts(); \r
1576 \r
1577   if(!evtCuts || !accCuts  || !esdTrackCuts) {\r
1578     AliDebug(AliLog::kError, "cuts not available");\r
1579     return;\r
1580   }\r
1581 \r
1582   // trigger selection\r
1583   Bool_t isEventTriggered = kTRUE;\r
1584   AliPhysicsSelection *physicsSelection = NULL;\r
1585   AliTriggerAnalysis* triggerAnalysis = NULL;\r
1586 \r
1587   // \r
1588   AliInputEventHandler* inputHandler = (AliInputEventHandler*) AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler();\r
1589   if (!inputHandler)\r
1590   {\r
1591     Printf("ERROR: Could not receive input handler");\r
1592     return;\r
1593   }\r
1594    \r
1595   // get file name\r
1596   TTree *chain = (TChain*)GetInputData(0);\r
1597   if(!chain) { \r
1598     Printf("ERROR: Could not receive input chain");\r
1599     return;\r
1600   }\r
1601   TObjString fileName(chain->GetCurrentFile()->GetName());\r
1602 \r
1603   // trigger\r
1604   if(evtCuts->IsTriggerRequired())  \r
1605   {\r
1606     // always MB\r
1607     isEventTriggered = inputHandler->IsEventSelected() & AliVEvent::kMB;\r
1608 \r
1609     physicsSelection = static_cast<AliPhysicsSelection*> (inputHandler->GetEventSelection());\r
1610     if(!physicsSelection) return;\r
1611     //SetPhysicsTriggerSelection(physicsSelection);\r
1612 \r
1613     if (isEventTriggered && (GetTrigger() == AliTriggerAnalysis::kV0AND)) {\r
1614       // set trigger (V0AND)\r
1615       triggerAnalysis = physicsSelection->GetTriggerAnalysis();\r
1616       if(!triggerAnalysis) return;\r
1617       isEventTriggered = triggerAnalysis->IsOfflineTriggerFired(esdEvent, GetTrigger());\r
1618     }\r
1619   }\r
1620 \r
1621   // centrality determination\r
1622   Float_t centralityF = -1;\r
1623   AliCentrality *esdCentrality = esdEvent->GetCentrality();\r
1624   centralityF = esdCentrality->GetCentralityPercentile(fCentralityEstimator.Data());\r
1625 \r
1626 \r
1627   // get reconstructed vertex  \r
1628   //const AliESDVertex* vtxESD = 0; \r
1629   AliESDVertex* vtxESD = 0; \r
1630   if(GetAnalysisMode() == kTPCAnalysisMode) {\r
1631         vtxESD = (AliESDVertex*)esdEvent->GetPrimaryVertexTPC();\r
1632   }\r
1633   else if(GetAnalysisMode() == kTPCITSAnalysisMode) {\r
1634      vtxESD = (AliESDVertex*)esdEvent->GetPrimaryVertexTracks();\r
1635   }\r
1636   else {\r
1637         return;\r
1638   }\r
1639 \r
1640   if(!vtxESD) return;\r
1641 \r
1642   Bool_t isEventOK = evtCuts->AcceptEvent(esdEvent,mcEvent,vtxESD); \r
1643   //printf("isEventOK %d, isEventTriggered %d \n",isEventOK, isEventTriggered);\r
1644   //printf("GetAnalysisMode() %d \n",GetAnalysisMode());\r
1645 \r
1646   // check event cuts\r
1647   if(isEventOK && isEventTriggered) {\r
1648   //\r
1649   // Dump the pt downscaled V0 into the tree\r
1650   // \r
1651   Int_t ntracks = esdEvent->GetNumberOfTracks();\r
1652   Int_t nV0s = esdEvent->GetNumberOfV0s();\r
1653   Int_t run = esdEvent->GetRunNumber();\r
1654   Int_t time= esdEvent->GetTimeStamp();\r
1655   Int_t evNr=esdEvent->GetEventNumberInFile();\r
1656   Double_t bz = esdEvent->GetMagneticField();\r
1657 \r
1658 \r
1659   for (Int_t iv0=0; iv0<nV0s; iv0++){\r
1660     AliESDv0 * v0 = esdEvent->GetV0(iv0);\r
1661     if (!v0) continue;\r
1662     AliESDtrack * track0 = esdEvent->GetTrack(v0->GetIndex(0));\r
1663     AliESDtrack * track1 = esdEvent->GetTrack(v0->GetIndex(1));\r
1664     if (!track0) continue;\r
1665     if (!track1) continue;\r
1666     if (track0->GetSign()<0) {\r
1667       track1 = esdEvent->GetTrack(v0->GetIndex(0));\r
1668       track0 = esdEvent->GetTrack(v0->GetIndex(1));\r
1669     }\r
1670     //\r
1671     Bool_t isDownscaled = IsV0Downscaled(v0);\r
1672     if (isDownscaled) continue;\r
1673     AliKFParticle kfparticle; //\r
1674     Int_t type=GetKFParticle(v0,esdEvent,kfparticle);\r
1675     if (type==0) continue;   \r
1676     TObjString triggerClass = esdEvent->GetFiredTriggerClasses().Data();\r
1677 \r
1678     if(!fTreeSRedirector) return;\r
1679     (*fTreeSRedirector)<<"V0s"<<\r
1680       "isDownscaled="<<isDownscaled<<\r
1681       "triggerClass="<<&triggerClass<<      //  trigger\r
1682       "Bz="<<bz<<\r
1683       "fileName.="<<&fileName<<\r
1684       "runNumber="<<run<<\r
1685       "evtTimeStamp="<<time<<\r
1686       "evtNumberInFile="<<evNr<<\r
1687       "type="<<type<<\r
1688       "ntracks="<<ntracks<<\r
1689       "v0.="<<v0<<\r
1690       "kf.="<<&kfparticle<<\r
1691       "track0.="<<track0<<\r
1692       "track1.="<<track1<<\r
1693       "centralityF="<<centralityF<<\r
1694       "\n";\r
1695   }\r
1696   }\r
1697 }\r
1698 \r
1699 //_____________________________________________________________________________\r
1700 void AliAnalysisTaskFilteredTree::ProcessdEdx(AliESDEvent *const esdEvent, AliMCEvent * const mcEvent, AliESDfriend *const /*esdFriend*/)\r
1701 {\r
1702   //\r
1703   // Select real events with large TPC dEdx signal\r
1704   //\r
1705   if(!esdEvent) {\r
1706     AliDebug(AliLog::kError, "esdEvent not available");\r
1707     return;\r
1708   }\r
1709 \r
1710   // get selection cuts\r
1711   AliFilteredTreeEventCuts *evtCuts = GetEventCuts(); \r
1712   AliFilteredTreeAcceptanceCuts *accCuts = GetAcceptanceCuts(); \r
1713   AliESDtrackCuts *esdTrackCuts = GetTrackCuts(); \r
1714 \r
1715   if(!evtCuts || !accCuts  || !esdTrackCuts) {\r
1716     AliDebug(AliLog::kError, "cuts not available");\r
1717     return;\r
1718   }\r
1719 \r
1720   // get file name\r
1721   TTree *chain = (TChain*)GetInputData(0);\r
1722   if(!chain) { \r
1723     Printf("ERROR: Could not receive input chain");\r
1724     return;\r
1725   }\r
1726   TObjString fileName(chain->GetCurrentFile()->GetName());\r
1727 \r
1728   // trigger\r
1729   Bool_t isEventTriggered = kTRUE;\r
1730   AliPhysicsSelection *physicsSelection = NULL;\r
1731   AliTriggerAnalysis* triggerAnalysis = NULL;\r
1732 \r
1733   // \r
1734   AliInputEventHandler* inputHandler = (AliInputEventHandler*) AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler();\r
1735   if (!inputHandler)\r
1736   {\r
1737     Printf("ERROR: Could not receive input handler");\r
1738     return;\r
1739   }\r
1740 \r
1741   if(evtCuts->IsTriggerRequired())  \r
1742   {\r
1743     // always MB\r
1744     isEventTriggered = inputHandler->IsEventSelected() & AliVEvent::kMB;\r
1745 \r
1746     physicsSelection = static_cast<AliPhysicsSelection*> (inputHandler->GetEventSelection());\r
1747     if(!physicsSelection) return;\r
1748 \r
1749     if (isEventTriggered && (GetTrigger() == AliTriggerAnalysis::kV0AND)) {\r
1750       // set trigger (V0AND)\r
1751       triggerAnalysis = physicsSelection->GetTriggerAnalysis();\r
1752       if(!triggerAnalysis) return;\r
1753       isEventTriggered = triggerAnalysis->IsOfflineTriggerFired(esdEvent, GetTrigger());\r
1754     }\r
1755   }\r
1756 \r
1757   // get reconstructed vertex  \r
1758   AliESDVertex* vtxESD = 0; \r
1759   if(GetAnalysisMode() == kTPCAnalysisMode) {\r
1760         vtxESD = (AliESDVertex*)esdEvent->GetPrimaryVertexTPC();\r
1761   }\r
1762   else if(GetAnalysisMode() == kTPCITSAnalysisMode) {\r
1763      vtxESD = (AliESDVertex*)esdEvent->GetPrimaryVertexTracks();\r
1764   }\r
1765   else {\r
1766         return;\r
1767   }\r
1768   if(!vtxESD) return;\r
1769 \r
1770   Bool_t isEventOK = evtCuts->AcceptEvent(esdEvent,mcEvent,vtxESD); \r
1771   //printf("isEventOK %d, isEventTriggered %d \n",isEventOK, isEventTriggered);\r
1772   //printf("GetAnalysisMode() %d \n",GetAnalysisMode());\r
1773 \r
1774 \r
1775   // check event cuts\r
1776   if(isEventOK && isEventTriggered)\r
1777   {\r
1778     Double_t vert[3] = {0}; \r
1779     vert[0] = vtxESD->GetXv();\r
1780     vert[1] = vtxESD->GetYv();\r
1781     vert[2] = vtxESD->GetZv();\r
1782     Int_t mult = vtxESD->GetNContributors();\r
1783     Double_t bz = esdEvent->GetMagneticField();\r
1784     Double_t runNumber = esdEvent->GetRunNumber();\r
1785     Double_t evtTimeStamp = esdEvent->GetTimeStamp();\r
1786     Int_t evtNumberInFile = esdEvent->GetEventNumberInFile();\r
1787 \r
1788     // large dEdx\r
1789     for (Int_t iTrack = 0; iTrack < esdEvent->GetNumberOfTracks(); iTrack++)\r
1790     {\r
1791       AliESDtrack *track = esdEvent->GetTrack(iTrack);\r
1792       if(!track) continue;\r
1793       if(track->Charge()==0) continue;\r
1794       if(!esdTrackCuts->AcceptTrack(track)) continue;\r
1795       if(!accCuts->AcceptTrack(track)) continue;\r
1796 \r
1797       if(!IsHighDeDxParticle(track)) continue;\r
1798       TObjString triggerClass = esdEvent->GetFiredTriggerClasses().Data();\r
1799 \r
1800       if(!fTreeSRedirector) return;\r
1801       (*fTreeSRedirector)<<"dEdx"<<\r
1802       "fileName.="<<&fileName<<\r
1803       "runNumber="<<runNumber<<\r
1804       "evtTimeStamp="<<evtTimeStamp<<\r
1805       "evtNumberInFile="<<evtNumberInFile<<\r
1806         "triggerClass="<<&triggerClass<<      //  trigger\r
1807       "Bz="<<bz<<\r
1808       "vtxESD.="<<vtxESD<<\r
1809       "mult="<<mult<<\r
1810       "esdTrack.="<<track<<\r
1811       "\n";\r
1812     }\r
1813   }\r
1814 }\r
1815 \r
1816 //_____________________________________________________________________________\r
1817 Int_t   AliAnalysisTaskFilteredTree::GetKFParticle(AliESDv0 *const v0, AliESDEvent * const event, AliKFParticle & kfparticle)\r
1818 {\r
1819   //\r
1820   // Create KF particle in case the V0 fullfill selection criteria\r
1821   //\r
1822   // Selection criteria\r
1823   //  0. algorithm cut\r
1824   //  1. track cut\r
1825   //  3. chi2 cut\r
1826   //  4. rough mass cut\r
1827   //  5. Normalized pointing angle cut\r
1828   //\r
1829   const Double_t cutMass=0.2;\r
1830   const Double_t kSigmaDCACut=3;\r
1831   //\r
1832   // 0.) algo cut - accept only on the fly\r
1833   //\r
1834   if (v0->GetOnFlyStatus() ==kFALSE) return 0;     \r
1835   //\r
1836   // 1.) track cut\r
1837   // \r
1838   AliESDtrack * track0 = event->GetTrack(v0->GetIndex(0));\r
1839   AliESDtrack * track1 = event->GetTrack(v0->GetIndex(1));\r
1840   /*\r
1841     TCut cutD="abs(track0.fD/sqrt(track0.fCdd))>2&&abs(track1.fD/sqrt(track1.fCdd))>2";\r
1842     TCut cutTheta="abs(track0.fP[3])<1&&abs(track1.fP[3])<1";\r
1843     TCut cutNcl="track0.GetTPCClusterInfo(2,1)>100&&track1.GetTPCClusterInfo(2,1)>100";\r
1844   */  \r
1845   if (TMath::Abs(track0->GetTgl())>1) return 0;\r
1846   if (TMath::Abs(track1->GetTgl())>1) return 0;\r
1847   if ((track0->GetTPCClusterInfo(2,1))<100) return 0;\r
1848   if ((track1->GetTPCClusterInfo(2,1))<100) return 0;\r
1849   //if ((track0->GetITSclusters(0))<2) return 0;\r
1850   //if ((track1->GetITSclusters(0))<2) return 0; \r
1851   Float_t pos0[2]={0}, cov0[3]={0};\r
1852   Float_t pos1[2]={0}, cov1[3]={0};\r
1853   track0->GetImpactParameters(pos0,cov0);\r
1854   track0->GetImpactParameters(pos1,cov1);\r
1855   //\r
1856   if (TMath::Abs(pos0[0])<kSigmaDCACut*TMath::Sqrt(cov0[0])) return 0;\r
1857   if (TMath::Abs(pos1[0])<kSigmaDCACut*TMath::Sqrt(cov1[0])) return 0;\r
1858   // \r
1859   //\r
1860   // 3.) Chi2 cut\r
1861   //\r
1862   Double_t chi2KF = v0->GetKFInfo(2,2,2);\r
1863   if (chi2KF>25) return 0;\r
1864   //\r
1865   // 4.) Rough mass cut - 0.200 GeV\r
1866   //\r
1867   static Double_t masses[2]={-1};\r
1868   if (masses[0]<0){\r
1869     masses[0] = TDatabasePDG::Instance()->GetParticle("K_S0")->Mass();\r
1870     masses[1] = TDatabasePDG::Instance()->GetParticle("Lambda0")->Mass();\r
1871   }\r
1872   Double_t mass00=  v0->GetEffMass(0,0);\r
1873   Double_t mass22=  v0->GetEffMass(2,2);\r
1874   Double_t mass42=  v0->GetEffMass(4,2);\r
1875   Double_t mass24=  v0->GetEffMass(2,4);\r
1876   Bool_t massOK=kFALSE;\r
1877   Int_t type=0;\r
1878   Int_t ptype=0;\r
1879   Double_t dmass=1;\r
1880   Int_t p1=0, p2=0;\r
1881   if (TMath::Abs(mass00-0)<cutMass) {\r
1882     massOK=kTRUE; type+=1; \r
1883     if (TMath::Abs(mass00-0)<dmass) {\r
1884       ptype=1;\r
1885       dmass=TMath::Abs(mass00-0);      \r
1886       p1=0; p2=0;\r
1887     } \r
1888   }\r
1889   if (TMath::Abs(mass24-masses[1])<cutMass) {\r
1890     massOK=kTRUE; type+=2; \r
1891     if (TMath::Abs(mass24-masses[1])<dmass){\r
1892       dmass = TMath::Abs(mass24-masses[1]);\r
1893       ptype=2;\r
1894       p1=2; p2=4;\r
1895     }\r
1896   }\r
1897   if (TMath::Abs(mass42-masses[1])<cutMass) {\r
1898     massOK=kTRUE; type+=4;\r
1899     if (TMath::Abs(mass42-masses[1])<dmass){\r
1900       dmass = TMath::Abs(mass42-masses[1]);\r
1901       ptype=4;\r
1902       p1=4; p2=2;\r
1903     }\r
1904   }\r
1905   if (TMath::Abs(mass22-masses[0])<cutMass) {\r
1906     massOK=kTRUE; type+=8;\r
1907     if (TMath::Abs(mass22-masses[0])<dmass){\r
1908       dmass = TMath::Abs(mass22-masses[0]);\r
1909       ptype=8;\r
1910       p1=2; p2=2;\r
1911     }\r
1912   }\r
1913   if (type==0) return 0;\r
1914   //\r
1915   const Int_t spdg[5]={kPositron,kMuonPlus,kPiPlus, kKPlus, kProton};\r
1916   const AliExternalTrackParam *paramP = v0->GetParamP();\r
1917   const AliExternalTrackParam *paramN = v0->GetParamN();\r
1918   if (paramP->GetSign()<0){\r
1919     paramP=v0->GetParamP();\r
1920     paramN=v0->GetParamN();\r
1921   }\r
1922   //Double_t *pparam1 = (Double_t*)paramP->GetParameter();\r
1923   //Double_t *pparam2 = (Double_t*)paramN->GetParameter();\r
1924   //\r
1925   AliKFParticle kfp1( *paramP, spdg[p1]  );\r
1926   AliKFParticle kfp2( *paramN, -1 *spdg[p2]  );\r
1927   AliKFParticle V0KF;\r
1928   (V0KF)+=kfp1;\r
1929   (V0KF)+=kfp2;\r
1930   kfparticle=V0KF;\r
1931   //\r
1932   // Pointing angle\r
1933   //\r
1934   Double_t  errPhi    = V0KF.GetErrPhi();\r
1935   Double_t  pointAngle= TMath::ACos(v0->GetV0CosineOfPointingAngle());\r
1936   if (pointAngle/errPhi>10) return 0;  \r
1937   //\r
1938   return ptype;  \r
1939 }\r
1940 \r
1941 //_____________________________________________________________________________\r
1942 Bool_t AliAnalysisTaskFilteredTree::IsV0Downscaled(AliESDv0 *const v0)\r
1943 {\r
1944   //\r
1945   // Downscale randomly low pt V0\r
1946   //\r
1947   //return kFALSE;\r
1948   Double_t maxPt= TMath::Max(v0->GetParamP()->Pt(), v0->GetParamN()->Pt());\r
1949   Double_t scalempt= TMath::Min(maxPt,10.);\r
1950   Double_t downscaleF = gRandom->Rndm();\r
1951   downscaleF *= fLowPtV0DownscaligF;\r
1952   //\r
1953   // Special treatment of the gamma conversion pt spectra is softer - \r
1954   Double_t mass00=  v0->GetEffMass(0,0);\r
1955   const Double_t cutMass=0.2;\r
1956   if (TMath::Abs(mass00-0)<cutMass){\r
1957     downscaleF/=10.;  // 10 times smaller downscaling for the gamma concersion candidate\r
1958   }\r
1959   //printf("V0 TMath::Exp(2*scalempt) %e, downscaleF %e \n",TMath::Exp(2*scalempt), downscaleF);\r
1960   if (TMath::Exp(2*scalempt)<downscaleF) return kTRUE;\r
1961   return kFALSE;\r
1962 \r
1963   /*\r
1964     TH1F his1("his1","his1",100,0,10);\r
1965     TH1F his2("his2","his2",100,0,10);\r
1966     {for (Int_t i=0; i<10000; i++){\r
1967        Double_t rnd=gRandom->Exp(1);\r
1968        Bool_t isDownscaled =TMath::Exp(rnd)<100*gRandom->Rndm();\r
1969        his1->Fill(rnd); \r
1970        if (!isDownscaled) his2->Fill(rnd); \r
1971     }}\r
1972 \r
1973    */\r
1974 \r
1975 }\r
1976 \r
1977 \r
1978 \r
1979 //_____________________________________________________________________________\r
1980 Bool_t AliAnalysisTaskFilteredTree::ConstrainTPCInner(AliExternalTrackParam *const tpcInnerC, const AliESDVertex* vtx, Double_t b[3])\r
1981 {\r
1982  // Constrain TPC inner params constrained\r
1983  //\r
1984       if(!tpcInnerC) return kFALSE; \r
1985       if(!vtx) return kFALSE;\r
1986 \r
1987       Double_t dz[2],cov[3];\r
1988       //AliESDVertex *vtx= (AliESDVertex *)esdEvent->GetPrimaryVertex();\r
1989       //if(!tpcInnerC->PropagateToDCA(vtx, esdEvent->GetMagneticField(), 3, dz, cov)) return kFALSE; \r
1990       //if(!tpcInnerC->PropagateToDCA(vtx, Bz, 3, dz, cov)) return kFALSE; \r
1991       if(!tpcInnerC->PropagateToDCABxByBz(vtx, b, 3, dz, cov)) return kFALSE; \r
1992 \r
1993 \r
1994       Double_t covar[6]; vtx->GetCovMatrix(covar);\r
1995       Double_t p[2]={tpcInnerC->GetParameter()[0]-dz[0],tpcInnerC->GetParameter()[1]-dz[1]};\r
1996       Double_t c[3]={covar[2],0.,covar[5]};\r
1997       Double_t chi2C=tpcInnerC->GetPredictedChi2(p,c);\r
1998       if (chi2C>kVeryBig) return kFALSE; \r
1999 \r
2000       if(!tpcInnerC->Update(p,c)) return kFALSE;\r
2001 \r
2002   return kTRUE;\r
2003 }\r
2004 \r
2005 //_____________________________________________________________________________\r
2006 Bool_t AliAnalysisTaskFilteredTree::ConstrainTrackInner(AliExternalTrackParam *const trackInnerC, const AliESDVertex* vtx, Double_t mass, Double_t b[3])\r
2007 {\r
2008  // Constrain TPC inner params constrained\r
2009  //\r
2010       if(!trackInnerC) return kFALSE; \r
2011       if(!vtx) return kFALSE;\r
2012 \r
2013       const Double_t kRadius  = 2.8; \r
2014       const Double_t kMaxStep = 1.0; \r
2015 \r
2016       Double_t dz[2],cov[3];\r
2017 \r
2018       //AliESDVertex *vtx= (AliESDVertex *)esdEvent->GetPrimaryVertex();\r
2019       //if(!trackInnerC->PropagateToDCA(vtx, esdEvent->GetMagneticField(), 3, dz, cov)) return kFALSE; \r
2020       //if(!trackInnerC->PropagateToDCA(vtx, Bz, 3, dz, cov)) return kFALSE; \r
2021 \r
2022       if(!AliTracker::PropagateTrackToBxByBz(trackInnerC,kRadius,mass,kMaxStep,kFALSE)) return kFALSE;\r
2023       if(!trackInnerC->PropagateToDCABxByBz(vtx, b, 3, dz, cov)) return kFALSE; \r
2024 \r
2025       //\r
2026       Double_t covar[6]; vtx->GetCovMatrix(covar);\r
2027       Double_t p[2]={trackInnerC->GetParameter()[0]-dz[0],trackInnerC->GetParameter()[1]-dz[1]};\r
2028       Double_t c[3]={covar[2],0.,covar[5]};\r
2029       Double_t chi2C=trackInnerC->GetPredictedChi2(p,c);\r
2030       if (chi2C>kVeryBig) return kFALSE; \r
2031 \r
2032       if(!trackInnerC->Update(p,c)) return kFALSE;\r
2033 \r
2034   return kTRUE;\r
2035 }\r
2036 \r
2037 \r
2038 //_____________________________________________________________________________\r
2039 TParticle *AliAnalysisTaskFilteredTree::GetMother(TParticle *const particle, AliStack *const stack) \r
2040 {\r
2041   if(!particle) return NULL;\r
2042   if(!stack) return NULL;\r
2043 \r
2044   Int_t motherLabel = TMath::Abs(particle->GetMother(0));  \r
2045   TParticle* mother = NULL; \r
2046   mother = stack->Particle(motherLabel); \r
2047 \r
2048 return mother;\r
2049 }\r
2050 \r
2051 //_____________________________________________________________________________\r
2052 Bool_t AliAnalysisTaskFilteredTree::IsFromConversion(const Int_t label, AliStack *const stack) \r
2053 {\r
2054   Bool_t isFromConversion = kFALSE;\r
2055 \r
2056   if(stack) {\r
2057     TParticle* particle = stack->Particle(label);\r
2058 \r
2059     if(particle && particle->GetPDG() && particle->GetPDG()->Charge()!=0) \r
2060     {\r
2061        Int_t mech = particle->GetUniqueID(); // production mechanism \r
2062        Bool_t isPrim = stack->IsPhysicalPrimary(label);\r
2063 \r
2064        Int_t motherLabel = TMath::Abs(particle->GetMother(0));  \r
2065        TParticle* mother = stack->Particle(motherLabel); \r
2066        if(mother) {\r
2067           Int_t motherPdg = mother->GetPdgCode();\r
2068 \r
2069           if(!isPrim && mech==5 && motherPdg==kGamma) { \r
2070             isFromConversion=kTRUE; \r
2071           }\r
2072        }\r
2073     } \r
2074   } \r
2075 \r
2076   return isFromConversion;\r
2077 }\r
2078 \r
2079 //_____________________________________________________________________________\r
2080 Bool_t AliAnalysisTaskFilteredTree::IsFromMaterial(const Int_t label, AliStack *const stack) \r
2081 {\r
2082   Bool_t isFromMaterial = kFALSE;\r
2083 \r
2084   if(stack) {\r
2085     TParticle* particle = stack->Particle(label);\r
2086 \r
2087     if(particle && particle->GetPDG() && particle->GetPDG()->Charge()!=0) \r
2088     {\r
2089        Int_t mech = particle->GetUniqueID(); // production mechanism \r
2090        Bool_t isPrim = stack->IsPhysicalPrimary(label);\r
2091 \r
2092        Int_t motherLabel = TMath::Abs(particle->GetMother(0));  \r
2093        TParticle* mother = stack->Particle(motherLabel); \r
2094        if(mother) {\r
2095           if(!isPrim && mech==13) { \r
2096             isFromMaterial=kTRUE; \r
2097           }\r
2098        }\r
2099      } \r
2100   } \r
2101 \r
2102   return isFromMaterial;\r
2103 }\r
2104 \r
2105 //_____________________________________________________________________________\r
2106 Bool_t AliAnalysisTaskFilteredTree::IsFromStrangeness(const Int_t label, AliStack *const stack) \r
2107 {\r
2108   Bool_t isFromStrangeness = kFALSE;\r
2109 \r
2110   if(stack) {\r
2111     TParticle* particle = stack->Particle(label);\r
2112 \r
2113     if(particle && particle->GetPDG() && particle->GetPDG()->Charge()!=0) \r
2114     {\r
2115        Int_t mech = particle->GetUniqueID(); // production mechanism \r
2116        Bool_t isPrim = stack->IsPhysicalPrimary(label);\r
2117 \r
2118        Int_t motherLabel = TMath::Abs(particle->GetMother(0));  \r
2119        TParticle* mother = stack->Particle(motherLabel); \r
2120        if(mother) {\r
2121           Int_t motherPdg = mother->GetPdgCode();\r
2122 \r
2123           // K+-, lambda, antilambda, K0s decays\r
2124           if(!isPrim && mech==4 && \r
2125               (TMath::Abs(motherPdg)==kKPlus || TMath::Abs(motherPdg)==kLambda0 || motherPdg==kK0Short))\r
2126           {\r
2127             isFromStrangeness = kTRUE;\r
2128           } \r
2129        }\r
2130     } \r
2131   } \r
2132 \r
2133   return isFromStrangeness;\r
2134 }\r
2135 \r
2136 \r
2137 //_____________________________________________________________________________\r
2138 void AliAnalysisTaskFilteredTree::FinishTaskOutput() \r
2139 {\r
2140   //\r
2141   // Called one at the end \r
2142   // locally on working node\r
2143   //\r
2144 \r
2145   //// must be deleted to store trees\r
2146   //if(fTreeSRedirector)  delete fTreeSRedirector; fTreeSRedirector=0;\r
2147 \r
2148   //// open temporary file and copy trees to the ouptut container\r
2149 \r
2150   //TChain* chain = 0;\r
2151   ////\r
2152   //chain = new TChain("highPt");\r
2153   //if(chain) { \r
2154   //  chain->Add("jotwinow_Temp_Trees.root");\r
2155   //  fHighPtTree = chain->CopyTree("1");\r
2156   //  delete chain; chain=0; \r
2157   //}\r
2158   //if(fHighPtTree) fHighPtTree->Print();\r
2159 \r
2160   ////\r
2161   //chain = new TChain("V0s");\r
2162   //if(chain) { \r
2163   //  chain->Add("jotwinow_Temp_Trees.root");\r
2164   //  fV0Tree = chain->CopyTree("1");\r
2165   //  delete chain; chain=0; \r
2166   //}\r
2167   //if(fV0Tree) fV0Tree->Print();\r
2168 \r
2169   ////\r
2170   //chain = new TChain("dEdx");\r
2171   //if(chain) { \r
2172   //  chain->Add("jotwinow_Temp_Trees.root");\r
2173   //  fdEdxTree = chain->CopyTree("1");\r
2174   //  delete chain; chain=0; \r
2175   //}\r
2176   //if(fdEdxTree) fdEdxTree->Print();\r
2177 \r
2178   ////\r
2179   //chain = new TChain("Laser");\r
2180   //if(chain) { \r
2181   //  chain->Add("jotwinow_Temp_Trees.root");\r
2182   //  fLaserTree = chain->CopyTree("1");\r
2183   //  delete chain; chain=0; \r
2184   //}\r
2185   //if(fLaserTree) fLaserTree->Print();\r
2186 \r
2187   ////\r
2188   //chain = new TChain("MCEffTree");\r
2189   //if(chain) { \r
2190   //  chain->Add("jotwinow_Temp_Trees.root");\r
2191   //  fMCEffTree = chain->CopyTree("1");\r
2192   //  delete chain; chain=0; \r
2193   //}\r
2194   //if(fMCEffTree) fMCEffTree->Print();\r
2195 \r
2196   ////\r
2197   //chain = new TChain("CosmicPairs");\r
2198   //if(chain) { \r
2199   //  chain->Add("jotwinow_Temp_Trees.root");\r
2200   //  fCosmicPairsTree = chain->CopyTree("1");\r
2201   //  delete chain; chain=0; \r
2202   //}\r
2203   //if(fCosmicPairsTree) fCosmicPairsTree->Print();  \r
2204 \r
2205   Bool_t deleteTrees=kTRUE;\r
2206   if ((AliAnalysisManager::GetAnalysisManager()))\r
2207   {\r
2208     if (AliAnalysisManager::GetAnalysisManager()->GetAnalysisType() == \r
2209              AliAnalysisManager::kProofAnalysis)\r
2210       deleteTrees=kFALSE;\r
2211   }\r
2212   if (deleteTrees) delete fTreeSRedirector;\r
2213   fTreeSRedirector=NULL;\r
2214 \r
2215   //OpenFile(1);\r
2216 \r
2217   // Post output data.\r
2218   //PostData(1, fHighPtTree);\r
2219   //PostData(2, fV0Tree);\r
2220   //PostData(3, fdEdxTree);\r
2221   //PostData(4, fLaserTree);\r
2222   //PostData(5, fMCEffTree);\r
2223   //PostData(6, fCosmicPairsTree);\r
2224 }\r
2225 \r
2226 //_____________________________________________________________________________\r
2227 void AliAnalysisTaskFilteredTree::Terminate(Option_t *) \r
2228 {\r
2229   // Called one at the end \r
2230   /*\r
2231   fOutputSummary = dynamic_cast<TTree*> (GetOutputData(1));\r
2232   if(fOutputSummary) delete fOutputSummary; fOutputSummary=0;\r
2233   TChain* chain = new TChain("highPt");\r
2234   if(!chain) return;\r
2235   chain->Add("jotwinow_HighPt_TrackAndV0_Trees.root");\r
2236   TTree *tree = chain->CopyTree("1");\r
2237   if (chain) { delete chain; chain=0; }\r
2238   if(!tree) return;\r
2239   tree->Print();\r
2240   fOutputSummary = tree;\r
2241 \r
2242   if (!fOutputSummary) {\r
2243     Printf("ERROR: AliAnalysisTaskFilteredTree::Terminate(): Output data not avaiable %p \n", GetOutputData(1));\r
2244     return;\r
2245   }\r
2246   */\r
2247   \r
2248 }\r
2249 \r
2250 //_____________________________________________________________________________\r
2251 Int_t AliAnalysisTaskFilteredTree::GetMCTrueTrackMult(AliMCEvent *const mcEvent, AliFilteredTreeEventCuts *const evtCuts, AliFilteredTreeAcceptanceCuts *const accCuts)\r
2252 {\r
2253   //\r
2254   // calculate mc event true track multiplicity\r
2255   //\r
2256   if(!mcEvent) return 0;\r
2257 \r
2258   AliStack* stack = 0;\r
2259   Int_t mult = 0;\r
2260 \r
2261   // MC particle stack\r
2262   stack = mcEvent->Stack();\r
2263   if (!stack) return 0;\r
2264 \r
2265   //\r
2266   //printf("minZv %f, maxZv %f \n", evtCuts->GetMinZv(), evtCuts->GetMaxZv());\r
2267   //\r
2268 \r
2269   Bool_t isEventOK = evtCuts->AcceptMCEvent(mcEvent);\r
2270   if(!isEventOK) return 0; \r
2271 \r
2272   Int_t nPart  = stack->GetNtrack();\r
2273   for (Int_t iMc = 0; iMc < nPart; ++iMc) \r
2274   {\r
2275      TParticle* particle = stack->Particle(iMc);\r
2276      if (!particle)\r
2277      continue;\r
2278 \r
2279      // only charged particles\r
2280      if(!particle->GetPDG()) continue;\r
2281      Double_t charge = particle->GetPDG()->Charge()/3.;\r
2282      if (TMath::Abs(charge) < 0.001)\r
2283      continue;\r
2284       \r
2285      // physical primary\r
2286      Bool_t prim = stack->IsPhysicalPrimary(iMc);\r
2287      if(!prim) continue;\r
2288 \r
2289      // checked accepted without pt cut\r
2290      //if(accCuts->AcceptTrack(particle)) \r
2291      if( particle->Eta() > accCuts->GetMinEta() && particle->Eta() < accCuts->GetMaxEta() ) \r
2292      {\r
2293        mult++;\r
2294      }\r
2295   }\r
2296 \r
2297 return mult;  \r
2298 }\r
2299 \r
2300 \r