93943dfccf29aa5bbc2e73635ab1ff65e1d374c1
[u/mrichter/AliRoot.git] / PWGLF / STRANGENESS / Hypernuclei / AliAODMCNuclExReplicator.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 //
17 // Implementation of a branch replicator 
18 // to produce aods with only few branches.
19 //
20 // This replicator is in charge of replicating the nuclei primary vertices
21 // tracks identified as nuclei with Z>=2, secondary vertices in form of 
22 // AliAODRecoDecayLF2Prong and their daughter tracks.
23 // These informations are stored into a reduced AODs (AliAOD.NuclEx.root) 
24 // 
25 // The vertices are filtered so that only the primary vertices make it
26 // to the output aods.
27 //
28 // The secondary vertices are recreated here, as a AliAODRecoDecayLF2Prong
29 // plus cuts that select secondary vericesoutside the primary vertex
30
31 // Authors: S. Bufalino (stefania.bufalino@cern.ch)
32 //          R. Lea      (ramona.lea@cern.ch)
33 // Based on AliAODMuonReplicator.cxx 
34
35 //NOTE : This is a test on MC : no PID response only MC truth + select only 3LH
36 // daughters : NOT INTENDED for any efficiency!!
37
38 class AliESDv0;
39 class AliESDVertex;
40 class AliAODVertex;
41 class AliAODRecoDecay;
42
43 #include "AliStack.h"
44 #include "AliAODEvent.h"
45 #include "AliAODMCHeader.h"
46 #include "AliAODMCParticle.h"
47 #include "AliAODTZERO.h"
48 #include "AliAODTrack.h"
49 #include "AliAODVZERO.h"
50 //#include "AliAnalysisCuts.h"
51 #include "TF1.h"
52 #include "AliExternalTrackParam.h"
53 #include "AliESDv0.h"
54 #include "AliAODv0.h"
55 //#include "AliPIDResponse.h"
56 #include <iostream>
57 #include <cassert>
58 #include "AliESDtrack.h"
59 #include "TObjArray.h"
60 #include "AliAnalysisFilter.h"
61 #include "AliAODRecoDecay.h"
62 #include "AliAODRecoDecayLF.h"
63 #include "AliAODRecoDecayLF2Prong.h"
64
65 #include <TFile.h>
66 #include <TDatabasePDG.h>
67 #include <TString.h>
68 #include <TList.h>
69 #include "AliLog.h"
70 #include "AliVEvent.h"
71 #include "AliVVertex.h"
72 #include "AliVTrack.h"
73 #include "AliVertexerTracks.h"
74 #include "AliKFVertex.h"
75 #include "AliESDEvent.h"
76 #include "AliESDVertex.h"
77 #include "AliESDtrackCuts.h"
78 #include "AliAODEvent.h"
79 #include "AliAnalysisFilter.h"
80 //#include "AliAnalysisVertexingLF.h"
81 #include "AliAnalysisManager.h"
82 #include "AliAODMCNuclExReplicator.h"
83 #include "TH1.h"
84 #include "TCanvas.h"
85 #include "AliInputEventHandler.h"
86
87 using std::cout;
88 using std::endl;
89
90 ClassImp(AliAODMCNuclExReplicator)
91
92 //_____________________________________________________________________________
93 AliAODMCNuclExReplicator::AliAODMCNuclExReplicator(const char* name, const char* title,
94                                                    Int_t mcMode, 
95                                                    Int_t nsigmaTrk1, Int_t partType1,
96                                                    Int_t nsigmaTrk2, Int_t partType2
97                                                    )
98   :AliAODBranchReplicator(name,title),
99   fBzkG(0.),
100   fCosMin(),
101   fDCAtracksMin(),
102   fRmax(),
103   fRmin(),
104   fDNmin(),
105   fDPmin(), 
106   fHeader(0x0),
107   fVertices(0x0), 
108   fNuclei(0x0),
109   fSecondaryVerices(0x0), 
110   fDaughterTracks(0x0),
111   fList(0x0),
112   fMCParticles(0x0),
113   fMCHeader(0x0),
114   fMCMode(mcMode),
115   fLabelMap(),
116   fParticleSelected(),
117   fReplicateHeader(kTRUE), //replicateHeader //to be fixed
118   fnSigmaTrk1(nsigmaTrk1),
119   fnSigmaTrk2(nsigmaTrk2),
120   fpartType1(partType1),
121   fpartType2(partType2),
122   fSecVtxWithKF(kFALSE),
123   fVertexerTracks(0x0),
124   fV1(0x0),
125   fAODMapSize(0),
126   fAODMap(0)
127   
128 {
129   // default ctor
130 }
131
132 //_____________________________________________________________________________
133 AliAODMCNuclExReplicator::~AliAODMCNuclExReplicator()
134 {
135   // destructor
136   // delete fTrackCut;
137   // delete fVertexCut;
138   delete fList;
139 }
140
141 //_____________________________________________________________________________
142 void AliAODMCNuclExReplicator::SelectParticle(Int_t i)
143 {
144   // taking the absolute values here, need to take care 
145   // of negative daughter and mother
146   // IDs when setting!
147   
148   if (!IsParticleSelected(TMath::Abs(i)))
149   {
150     fParticleSelected.Add(TMath::Abs(i),1);    
151   }
152 }
153
154 //_____________________________________________________________________________
155 Bool_t AliAODMCNuclExReplicator::IsParticleSelected(Int_t i)  
156 {
157   // taking the absolute values here, need to take 
158   // care with negative daughter and mother
159   // IDs when setting!
160   return (fParticleSelected.GetValue(TMath::Abs(i))==1);
161 }
162
163
164 //_____________________________________________________________________________
165 void AliAODMCNuclExReplicator::CreateLabelMap(const AliAODEvent& source)
166 {  
167   //
168   // this should be called once all selections are done 
169   //
170   
171   fLabelMap.Delete();
172   
173   TClonesArray* mcParticles = static_cast<TClonesArray*>(source.FindListObject(AliAODMCParticle::StdBranchName()));
174   
175   Int_t i(0);
176   Int_t j(0);
177   
178   TIter next(mcParticles);
179   
180   while ( next() ) 
181   {
182     if (IsParticleSelected(i))
183     {
184       fLabelMap.Add(i,j++);
185     }
186     ++i;
187   }
188 }
189
190 //_____________________________________________________________________________
191 Int_t AliAODMCNuclExReplicator::GetNewLabel(Int_t i) 
192 {
193   // Gets the label from the new created Map
194   // Call CreatLabelMap before
195   // otherwise only 0 returned
196   return fLabelMap.GetValue(TMath::Abs(i));
197 }
198
199
200 //_____________________________________________________________________________
201 TList* AliAODMCNuclExReplicator::GetList() const
202 {
203   // return (and build if not already done) our internal list of managed objects
204   if (!fList)
205     {
206     
207       if ( fReplicateHeader )
208         {
209           fHeader = new AliAODHeader;
210         }
211       
212       
213       fSecondaryVerices = new TClonesArray("AliAODRecoDecayLF2Prong",30);
214       fSecondaryVerices->SetName("SecondaryVertices");    
215     
216       fVertices = new TClonesArray("AliAODVertex",2);
217       fVertices->SetName("vertices");    
218     
219       fNuclei = new TClonesArray("AliAODTrack",30);
220       fNuclei->SetName("Nuclei");
221     
222       fDaughterTracks = new TClonesArray("AliAODTrack",30);
223       fDaughterTracks->SetName("DaughterTracks");
224
225     
226       fList = new TList;
227       fList->SetOwner(kTRUE);
228   
229       fList->Add(fHeader);
230       fList->Add(fVertices);
231       fList->Add(fNuclei);
232       fList->Add(fSecondaryVerices);
233       fList->Add(fDaughterTracks);
234   
235       
236       if ( fMCMode > 0 )
237         {
238           fMCHeader = new AliAODMCHeader;    
239           fMCParticles = new TClonesArray("AliAODMCParticle",1000);
240           fMCParticles->SetName(AliAODMCParticle::StdBranchName());
241           fList->Add(fMCHeader);
242           fList->Add(fMCParticles);
243         }
244     }
245   return fList;
246 }
247
248 //_____________________________________________________________________________
249 void AliAODMCNuclExReplicator::ReplicateAndFilter(const AliAODEvent& source)
250 {
251   // Replicate (and filter if filters are there) the relevant parts we're interested in AODEvent
252   
253   //  cout<<"-------------------->QUESTO"<<endl;
254
255   //-----------------------------------------------
256   // AliPIDResponse
257   
258   // AliAnalysisManager *man=AliAnalysisManager::GetAnalysisManager();
259   // AliInputEventHandler* inputHandler = (AliInputEventHandler*) (man->GetInputEventHandler());
260     
261   //--------------------------------------------------------
262
263   //  printf("Execute NuclEx Replicator\n");
264
265   //---------------------------------
266
267   if (fReplicateHeader)
268     {
269       *fHeader = *(source.GetHeader());
270     }
271     
272   fVertices->Clear("C");                        
273     
274   fNuclei->Clear("C");
275
276   fSecondaryVerices->Clear("C");
277
278   fDaughterTracks->Clear("C");
279   
280   //----------------------------------
281   
282   //retrive MC infos
283   
284   TClonesArray *arrayMC = 0;
285   AliAODMCHeader *mcHeader=0;
286   Int_t mumpdg=-100;
287   //  Int_t mumpdgNeg=-100;
288   
289
290   arrayMC = (TClonesArray*) source.GetList()->FindObject(AliAODMCParticle::StdBranchName());
291   if (!arrayMC) {
292     Printf("Error: MC particles branch not found!\n");
293     return;
294   }
295   
296   // if(arrayMC)
297   //   cout<<"Ho caricato array mc"<<endl;
298   
299   mcHeader =  (AliAODMCHeader*)source.GetList()->FindObject(AliAODMCHeader::StdBranchName());
300   if(!mcHeader) {
301     printf("AliAnalysisTaskSEDs::UserExec: MC header branch not found!\n");
302     return;
303   } 
304
305   // if(mcHeader)
306   //   cout<<"Ho caricato MC header"<<endl;
307   
308   //  cout<<"Centrality AOD source: "<<source.GetHeader()->GetCentrality()<<endl;
309
310   Int_t nsv(0);
311   //  Int_t nnuclei(0);
312   Int_t ntracksD(0);
313
314   Int_t input(0);
315   Double_t xdummy,ydummy;
316
317   AliAODRecoDecayLF2Prong *io2Prong  = 0;
318
319   TObjArray *twoTrackArray    = new TObjArray(2);
320   Double_t dispersion;
321
322   // cout<<"Qui"<<endl;
323   // cout<<source.GetMagneticField()<<endl;
324
325   AliAODVertex *vtx = source.GetPrimaryVertex();
326                                                 
327   // cout<<"Source "<<source<<endl;
328   // cout<<"vtx: "<<vtx<<endl;
329
330   // A Set of very loose cut for a weak strange decay
331   
332   fCosMin       = 0.99;
333   fDCAtracksMin = 1;
334   fRmax         = 200.;
335   fRmin         = 0.1;
336   fDNmin        = 0.05;
337   fDPmin        = 0.05;
338
339   //----------------------------------------------------------
340
341   //  Int_t nindices=0;
342   UShort_t *indices = 0;
343   const Int_t entries = source.GetNumberOfTracks();
344
345   Double_t pos[3],cov[6];
346   vtx->GetXYZ(pos);
347   vtx->GetCovarianceMatrix(cov);
348   fV1 = new AliESDVertex(pos,cov,100.,100,vtx->GetName());
349   cout<<"fV1 pricipal loop: "<<fV1<<endl;
350   
351   if(entries<=0) return;
352   indices = new UShort_t[entries];
353   memset(indices,0,sizeof(UShort_t)*entries);
354   fAODMapSize = 100000;
355   fAODMap = new Int_t[fAODMapSize];
356   memset(fAODMap,0,sizeof(Int_t)*fAODMapSize);
357   //  cent=((AliAODEvent&)source)->GetCentrality();
358   
359   //-------------------------------------------------------------
360
361   if(vtx->GetNContributors()<1) {
362     
363     // SPD vertex cut
364     vtx =source.GetPrimaryVertexSPD();
365     
366     if(vtx->GetNContributors()<1) {
367       Info("AliAnalysisTaskHelium3Pi","No good vertex, skip event");
368       return; // NO GOOD VERTEX, SKIP EVENT 
369     }
370   }
371   
372   Double_t xPrimaryVertex=0.,yPrimaryVertex=0.;
373   xPrimaryVertex=vtx->GetX();
374   yPrimaryVertex=vtx->GetY();
375   
376   fBzkG=source.GetMagneticField();
377   fVertexerTracks=new AliVertexerTracks(fBzkG);
378
379   Double_t TrackNumber = source.GetNumberOfTracks();
380   Int_t label =-1;
381
382   //Tracks arrays
383   
384   TArrayI Track0(TrackNumber);        //Pions                                                                          
385   Int_t nTrack0=0;
386   
387   TArrayI Track1(TrackNumber);        //Helium3
388   Int_t nTrack1=0;
389
390   for(Int_t j=0; j<TrackNumber; j++){
391
392     //    cout<<"Inside loop tracks"<<endl;
393
394   
395     AliVTrack *track = (AliVTrack*)source.GetTrack(j);
396     
397     AliAODTrack *aodtrack =(AliAODTrack*)track;
398
399     //-----------------------------------------------------------
400     //Track cuts 
401     if(aodtrack->GetTPCNcls() < 70 )continue;
402     if(aodtrack->Chi2perNDF() > 4 )continue;
403     
404     if (!aodtrack->IsOn(AliAODTrack::kTPCrefit)) continue;
405     if (!aodtrack->IsOn(AliAODTrack::kTPCin)) continue;
406     if (aodtrack->IsOn(AliAODTrack::kITSpureSA)) continue;
407
408     //---------------------------------------------------------------
409      
410     Double_t mom = aodtrack->P();
411     
412     if(mom<0.150)continue;
413
414     label = TMath::Abs(aodtrack->GetLabel());
415
416     AliAODMCParticle *part = (AliAODMCParticle*) arrayMC->At(label);
417     
418     Int_t PDGCode=part->GetPdgCode();
419     
420     Int_t mumid = part->GetMother();
421     
422     if(mumid>-1){
423       AliAODMCParticle *mother = (AliAODMCParticle*) arrayMC->At(mumid);
424       mumpdg = mother->GetPdgCode();
425     }
426     
427     //    if(mumpdg == 1010010030 ||mumpdg == -1010010030 ){
428
429     if(mumpdg == -1010010030){
430
431       if(PDGCode==-211 || PDGCode==+211){  
432         //      if(PDGCode==-211){  
433         
434         //      if(aodtrack->Charge()>0)continue;
435         
436         Track0[nTrack0++]=j;
437       }
438       
439       if(PDGCode==1000020030 ||PDGCode==-1000020030 ){
440       //if(PDGCode==1000020030){
441         
442         //      if(aodtrack->Charge()<0)continue;
443         
444         Track1[nTrack1++]=j;
445         //      new((*fNuclei)[nnuclei++]) AliAODTrack(*aodtrack);
446       }
447     }
448   }
449   
450   //Set Track Daughters
451   Track0.Set(nTrack0);
452   Track1.Set(nTrack1);
453   
454   
455   // cout<<"Track loops..."<<endl;
456   // cout<<"npos "<<nTrack1<<endl;
457   // cout<<"nneg "<<nTrack0<<endl;
458
459
460   AliAODTrack *postrackAOD = 0;
461   AliAODTrack *negtrackAOD = 0;
462  
463   AliESDtrack *postrack = 0;
464   AliESDtrack *negtrack = 0;
465
466   Bool_t isOk=kFALSE;
467
468   for (Int_t i=0; i<nTrack1; i++){                            //! He Tracks Loop
469     
470     Int_t Track1idx=Track1[i];
471
472     AliVTrack *trackpos = (AliVTrack*)source.GetTrack(Track1idx);
473
474     postrackAOD = (AliAODTrack*)trackpos;
475     postrack = new AliESDtrack(trackpos);
476
477     //--- MC infos
478
479     Int_t labelpos = TMath::Abs(postrack->GetLabel());
480     AliAODMCParticle *partPos = (AliAODMCParticle*) arrayMC->At(labelpos);
481     Int_t mumidPos = partPos->GetMother();
482
483     //    cout<<"\n\nmumidPos: "<<mumidPos <<endl;
484
485     //------------------------------
486     
487     for (Int_t k=0; k <nTrack0 ; k++) {                           //! Pion Tracks Loop
488       
489       Int_t Track0idx=Track0[k];
490       
491       AliVTrack *trackneg = (AliVTrack*)source.GetTrack(Track0idx);
492       negtrackAOD =(AliAODTrack*)trackneg;
493       negtrack = new AliESDtrack(trackneg); 
494       
495       //--- MC infos
496       
497       Int_t labelneg = TMath::Abs(negtrack->GetLabel());
498       AliAODMCParticle *partNeg = (AliAODMCParticle*) arrayMC->At(labelneg);
499       Int_t mumidNeg = partNeg->GetMother();
500       
501       //      cout<<"mumidNeg: "<<mumidNeg <<endl;
502
503       // AliAODMCParticle *motherNeg = (AliAODMCParticle*) arrayMC->At(mumidNeg);
504       // mumpdgNeg = motherNeg->GetPdgCode();
505
506       //------------------------------
507       //  if(mumidPos == mumidNeg && mumidNeg > 0){
508       isOk=kFALSE;
509       
510       // if(mumidPos == mumidNeg && mumidNeg > 0 && mumpdgNeg  == 1010010030)
511       if(mumidPos == mumidNeg && mumidNeg > 0)
512         isOk = kTRUE;
513       
514       twoTrackArray->AddAt(negtrack,0);
515       twoTrackArray->AddAt(postrack,1);
516       
517       Double_t dcap1n1 = postrack->GetDCA(negtrack,fBzkG,xdummy,ydummy);
518       
519       Double_t dcap1toPV = TMath::Abs(postrack->GetD(xPrimaryVertex, yPrimaryVertex,fBzkG));
520       Double_t dcan1toPV = TMath::Abs(negtrack->GetD(xPrimaryVertex, yPrimaryVertex,fBzkG));
521       
522       if(dcap1n1>fDCAtracksMin)continue;
523       if((xdummy+ydummy)>fRmax )continue;
524       if((xdummy+ydummy)< fRmin)continue;
525         
526       if ( dcan1toPV < fDNmin)               
527         if ( dcap1toPV < fDPmin) continue;   
528       
529       //      cout<<"dcap1n1: "<<dcap1n1<<endl;
530         
531       AliAODVertex *vertexp1n1 = ReconstructSecondaryVertex(twoTrackArray,dispersion,kTRUE);
532         
533       if(!vertexp1n1) {
534           
535         twoTrackArray->Clear();
536         delete negtrack;
537         negtrack=NULL; 
538         continue; 
539       }
540         
541       io2Prong = Make2Prong(twoTrackArray,source,vertexp1n1,dcap1n1);
542         
543       if(io2Prong->CosPointingAngle()<fCosMin)continue;
544         
545       //      cout<<"CPA: "<<io2Prong->CosPointingAngle()<<endl;
546       
547       // AliAODTrack *trk0 = (AliAODTrack*)io2Prong->GetDaughter(0);
548       // AliAODTrack *trk1 = (AliAODTrack*)io2Prong->GetDaughter(1);
549         
550       // cout<<"**********************************************"<<endl;
551       // cout<<trk0/*->GetID()*/<<" "<<negtrackAOD->GetID()<<endl;
552       // cout<<trk1/*->GetID()*/<<" "<<postrackAOD->GetID()<<endl;
553       // cout<<"d0 io2Prong: "<<io2Prong->GetProngID(1)<<endl;
554       // cout<<"d1 io2Prong: "<<io2Prong->GetProngID(0)<<endl;
555       // cout<<"**********************************************"<<endl;
556         
557       //      rd =  new((*fSecondaryVerices)[nsv++]) AliAODRecoDecayLF2Prong(*io2Prong);
558
559       //      cout<<"is ok??? "<<isOk<<endl;
560       if(isOk){
561         new((*fSecondaryVerices)[nsv++]) AliAODRecoDecayLF2Prong(*io2Prong);
562         new((*fDaughterTracks)[ntracksD++]) AliAODTrack(*negtrackAOD);
563         new((*fDaughterTracks)[ntracksD++]) AliAODTrack(*postrackAOD);
564         
565       }
566       // rd->SetSecondaryVtx(vertexp1n1);
567       // vertexp1n1->SetParent(rd);
568       // AddRefs(vertexp1n1,rd,source,twoTrackArray);
569       
570       delete negtrack; 
571       negtrack = NULL;
572       
573       delete vertexp1n1;
574       vertexp1n1= NULL;
575       continue;
576       //   }
577     }
578     
579     delete postrack; 
580     postrack = NULL;
581     
582   }
583   
584   //----------------------------------------------------------
585  
586   assert(fVertices!=0x0);
587   fVertices->Clear("C");
588   TIter nextV(source.GetVertices());
589   AliAODVertex* v;
590   Int_t nvertices(0);
591   
592   while ( ( v = static_cast<AliAODVertex*>(nextV()) ) )
593     {
594       if ( v->GetType() == AliAODVertex::kPrimary     ||
595            v->GetType() == AliAODVertex::kMainSPD     ||
596            v->GetType() == AliAODVertex::kPileupSPD   ||
597            v->GetType() == AliAODVertex::kPileupTracks||
598            v->GetType() == AliAODVertex::kMainTPC  ) 
599         {
600                   
601           AliAODVertex* tmp = v->CloneWithoutRefs();
602           //AliAODVertex* copiedVertex = new((*fVertices)[nvertices++]) AliAODVertex(*tmp);
603           new((*fVertices)[nvertices++]) AliAODVertex(*tmp);
604           
605           // to insure the main vertex retains the ncontributors information
606           // (which is otherwise computed dynamically from
607           // references to tracks, which we do not keep in muon aods...)
608           // we set it here
609           
610           //copiedVertex->SetNContributors(v->GetNContributors()); 
611           
612           //  fVertices->Delete();
613           //      delete copiedVertex;
614           delete tmp;
615           //      printf("....Prendo il vertice primario...\n");
616         }
617       //       printf("....Prendo il vertice primario...\n");
618     }
619   
620   //printf("....Done NuclEx Replicator...\n");
621   
622   AliDebug(1,Form("input mu tracks=%d tracks=%d vertices=%d nnuclei=%d",
623                   input,fSecondaryVerices->GetEntries(),fVertices->GetEntries(),fNuclei->GetEntries())); 
624   
625   // cout<<"Delete..."<<endl;
626   // delete foPion;
627   // foPion = NULL;
628   //cout<<"Delete 1"<<  endl;
629   
630   if(io2Prong) {delete io2Prong; io2Prong=NULL;}
631   //cout<<"Delete 2"<<  endl;
632   twoTrackArray->Delete();  delete twoTrackArray;
633   //  cout<<"Delete 3"<<  endl;
634   // vtx->Delete();  delete vtx;
635   //  cout<<"Delete 4"<<  endl;
636   if(fV1) { delete fV1; fV1=NULL; }
637   //  cout<<"Delete 5"<<  endl;
638   if(fAODMap) { delete [] fAODMap; fAODMap=NULL; }
639   delete []indices;
640   //  cout<<"Delete 6"<<  endl;
641   delete fVertexerTracks;
642
643   //  cout<<"Mi rompo alla fine. Perche???"<<endl;
644   
645
646
647 //-----------------------------------------------------------------------------
648
649 AliAODVertex *AliAODMCNuclExReplicator::ReconstructSecondaryVertex(TObjArray *trkArray,
650                                                                 Double_t &dispersion,Bool_t useTRefArray) const
651
652 {
653   // Secondary vertex reconstruction with AliVertexerTracks or AliKFParticle
654   //AliCodeTimerAuto("",0);
655
656   AliESDVertex *vertexESD = 0;
657   AliAODVertex *vertexAOD = 0;
658
659   if(!fSecVtxWithKF) { // AliVertexerTracks
660
661     fVertexerTracks->SetVtxStart(fV1);
662     vertexESD = (AliESDVertex*)fVertexerTracks->VertexForSelectedESDTracks(trkArray);
663
664     if(!vertexESD) return vertexAOD;
665
666
667     Double_t vertRadius2=vertexESD->GetXv()*vertexESD->GetXv()+vertexESD->GetYv()*vertexESD->GetYv();
668     if(vertRadius2>200.){
669       // vertex outside beam pipe, reject candidate to avoid propagation through material
670       delete vertexESD; vertexESD=NULL;
671       return vertexAOD;
672     }
673
674   } else { // Kalman Filter vertexer (AliKFParticle)
675
676     AliKFParticle::SetField(fBzkG);
677
678     AliKFVertex vertexKF;
679
680     Int_t nTrks = trkArray->GetEntriesFast();
681     for(Int_t i=0; i<nTrks; i++) {
682       AliESDtrack *esdTrack = (AliESDtrack*)trkArray->At(i);
683       AliKFParticle daughterKF(*esdTrack,211);
684       vertexKF.AddDaughter(daughterKF);
685     }
686     vertexESD = new AliESDVertex(vertexKF.Parameters(),
687                                  vertexKF.CovarianceMatrix(),
688                                  vertexKF.GetChi2(),
689                                  vertexKF.GetNContributors());
690
691   }
692   // convert to AliAODVertex
693   Double_t pos[3],cov[6],chi2perNDF;
694   vertexESD->GetXYZ(pos); // position
695   vertexESD->GetCovMatrix(cov); //covariance matrix
696   chi2perNDF = vertexESD->GetChi2toNDF();
697   dispersion = vertexESD->GetDispersion();
698   delete vertexESD; 
699   vertexESD=NULL;
700
701   Int_t nprongs= (useTRefArray ? 0 : trkArray->GetEntriesFast());
702   vertexAOD = new AliAODVertex(pos,cov,chi2perNDF,0x0,-1,AliAODVertex::kUndef,nprongs);
703
704   //  cout<<"------------------> Reconstruct vertexAOD: "<<vertexAOD<<endl;
705
706   return vertexAOD;
707
708  
709 }
710
711 //-----------------------------------------------------------------------------
712
713 AliAODRecoDecayLF2Prong* AliAODMCNuclExReplicator::Make2Prong(TObjArray *twoTrackArray,const AliAODEvent &evento,
714                                                            AliAODVertex *secVert,Double_t dca)
715
716
717
718
719   //cout<<"Inside make 2 prong"<<endl;
720
721   Int_t charge[2];
722   charge[0]=1; //it was -1 //Ramona
723   charge[1]=2;
724       
725   // From AliAnalysisVertexingLF.cxx Method:Make2Prongs
726   
727   fBzkG = evento.GetMagneticField();
728       
729   Double_t px[2],py[2],pz[2],d0[2],d0err[2];
730   // Also this has been changed //Ramona
731   AliESDtrack *negtrack = (AliESDtrack*)twoTrackArray->UncheckedAt(0);
732   AliESDtrack *postrack = (AliESDtrack*)twoTrackArray->UncheckedAt(1);
733
734   //cout<<"negtrack: "<<negtrack<<" postrack: "<<postrack<<endl;
735   //cout<<"kVeryBig: "<<kVeryBig<<endl;
736   //cout<<"secVert: "<<secVert<<endl;
737
738   // // propagate tracks to secondary vertex, to compute inv. mass
739   
740   postrack->PropagateToDCA(secVert,fBzkG,kVeryBig);
741   negtrack->PropagateToDCA(secVert,fBzkG,kVeryBig);
742   
743   Double_t momentum[3];
744   
745   negtrack->GetPxPyPz(momentum);
746   px[0] = charge[0]*momentum[0]; 
747   py[0] = charge[0]*momentum[1]; 
748   pz[0] = charge[0]*momentum[2]; 
749
750   postrack->GetPxPyPz(momentum);
751   px[1] = charge[1]*momentum[0]; 
752   py[1] = charge[1]*momentum[1]; 
753   pz[1] = charge[1]*momentum[2]; 
754   
755   //cout<< px[0] <<" "<< " "<< py[0] << " "<< pz[0]<<endl;
756   //  px[0] = momentum[0]; py[0] = momentum[1]; pz[0] = momentum[2]; 
757   //cout<< px[1] <<" "<< " "<< py[1] << " "<< pz[1]<<endl;
758   //px[1] = momentum[0]; py[1] = momentum[1]; pz[1] = momentum[2]; 
759   
760
761   // primary vertex to be used by this candidate
762   AliAODVertex *primVertexAOD  = evento.GetPrimaryVertex();
763   //cout<<"primVertexAOD "<<primVertexAOD<<endl;
764   if(!primVertexAOD) return 0x0;
765       
766   Double_t d0z0[2],covd0z0[3];
767
768   d0z0[0] = -999.;
769   d0z0[1] = -999.;
770   covd0z0[0] = -999.;
771   covd0z0[1] = -999.;
772   covd0z0[2] = -999.;
773
774   d0[0] = d0z0[0];
775   d0err[0] = TMath::Sqrt(TMath::Abs(covd0z0[0]));
776
777   d0[1] = d0z0[0];
778   d0err[1] = TMath::Sqrt(TMath::Abs(covd0z0[0]));
779   
780   // create the object AliAODRecoDecayLF2Prong
781   //  AliAODRecoDecayLF2Prong *the2Prong = new AliAODRecoDecayLF2Prong(secVert,px,py,pz,d0,d0err,dcap1n1);
782   AliAODRecoDecayLF2Prong *the2Prong = new AliAODRecoDecayLF2Prong(secVert,px,py,pz,d0,d0err,dca);
783   the2Prong->SetOwnPrimaryVtx(primVertexAOD);
784   
785   //  the2Prong->SetProngIDs(2,{-999,999});
786   // UShort_t id0[2]={99999,999999};
787   // the2Prong->SetProngIDs(2,id0);
788
789   UShort_t id[2]={(UShort_t)postrack->GetID(),(UShort_t)negtrack->GetID()};
790   the2Prong->SetProngIDs(2,id);
791
792   // cout<<"\n\n\nMake 2 Prong: id[0]"<<id[0]<<" id[1]: "<<id[1]<<endl;
793   // cout<<"Get: 1 "<<the2Prong->GetProngID(0)<<" 2 "<<the2Prong->GetProngID(1)<<endl;
794   // cout<<"Get: 3 "<<the2Prong->GetProngID(2)<<" 4 "<<the2Prong->GetProngID(3)<<endl<<endl<<endl;
795   //delete primVertexAOD; primVertexAOD=NULL;
796   
797   if(postrack->Charge()!=0 && negtrack->Charge()!=0) {
798       
799     AddDaughterRefs(secVert,(AliAODEvent&)evento,twoTrackArray);
800     //    AddDaughterRefs(secVert,(AliAODEvent*)evento,twoTrackArray);
801       
802   }
803   
804   return the2Prong;  
805
806   delete the2Prong;
807 }
808
809 //----------------------------------------------------------------------------
810 void AliAODMCNuclExReplicator::AddDaughterRefs(AliAODVertex *v,
811                                             const AliAODEvent &event,
812                                             const TObjArray *trkArray) const
813
814 {
815   // Add the AOD tracks as daughters of the vertex (TRef)
816   // AliCodeTimerAuto("",0);
817   // cout<<"Inside  AddDaughterRefs "<<endl;
818
819   Int_t nDg = v->GetNDaughters();
820   
821   TObject *dg = 0;
822   if(nDg) dg = v->GetDaughter(0);
823   if(dg) return; // daughters already added
824   
825   Int_t nTrks = trkArray->GetEntriesFast();
826   
827   AliExternalTrackParam *track = 0;
828   AliAODTrack *aodTrack = 0;
829   Int_t id;
830   
831   for(Int_t i=0; i<nTrks; i++) {
832     track = (AliExternalTrackParam*)trkArray->UncheckedAt(i);
833   
834     id = (Int_t)track->GetID();
835     //    printf("---> %d\n",id);
836     if(id<0) continue; // this track is a AliAODRecoDecay
837   
838     aodTrack = (AliAODTrack*)event.GetTrack(fAODMap[id]);
839     v->AddDaughter(aodTrack);
840   }
841   return;
842 }
843 //----------------------------------------------------------------------------
844         
845 void AliAODMCNuclExReplicator::AddRefs(AliAODVertex *v,AliAODRecoDecayLF *rd,
846                                     const AliAODEvent &event,
847                                     const TObjArray *trkArray) const
848
849 {
850   // Add the AOD tracks as daughters of the vertex (TRef)
851   // Also add the references to the primary vertex and to the cuts
852   //AliCodeTimerAuto("",0);
853   
854   AddDaughterRefs(v,event,trkArray);
855   rd->SetPrimaryVtxRef((AliAODVertex*)event.GetPrimaryVertex());
856   return;
857 }       
858
859 //---------------------------------------------------------------------------
860
861 void AliAODMCNuclExReplicator::Terminate(){
862
863 }