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