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