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