cd338d9062c266bd0311425edf783a53e4de6a07
[u/mrichter/AliRoot.git] / PWGLF / STRANGENESS / Hypernuclei / AliAODNuclExReplicator.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
36 class AliESDv0;
37 class AliESDVertex;
38 class AliAODVertex;
39 class AliAODRecoDecay;
40
41 #include "AliAODDimuon.h"
42 #include "AliAODEvent.h"
43 #include "AliAODMCHeader.h"
44 #include "AliAODMCParticle.h"
45 #include "AliAODTZERO.h"
46 #include "AliAODTrack.h"
47 #include "AliAODVZERO.h"
48 //#include "AliAnalysisCuts.h"
49 #include "TF1.h"
50 #include "AliExternalTrackParam.h"
51 #include "AliESDv0.h"
52 #include "AliAODv0.h"
53 #include "AliPIDResponse.h"
54 #include <iostream>
55 #include <cassert>
56 #include "AliESDtrack.h"
57 #include "TObjArray.h"
58 #include "AliAnalysisFilter.h"
59 #include "AliAODRecoDecay.h"
60 #include "AliAODRecoDecayLF.h"
61 #include "AliAODRecoDecayLF2Prong.h"
62
63 #include <TFile.h>
64 #include <TDatabasePDG.h>
65 #include <TString.h>
66 #include <TList.h>
67 #include "AliLog.h"
68 #include "AliVEvent.h"
69 #include "AliVVertex.h"
70 #include "AliVTrack.h"
71 #include "AliVertexerTracks.h"
72 #include "AliKFVertex.h"
73 #include "AliESDEvent.h"
74 #include "AliESDVertex.h"
75 #include "AliESDtrackCuts.h"
76 #include "AliAODEvent.h"
77 #include "AliAnalysisFilter.h"
78 //#include "AliAnalysisVertexingLF.h"
79 #include "AliAnalysisManager.h"
80 #include "AliAODNuclExReplicator.h"
81 #include "TH1.h"
82 #include "TCanvas.h"
83 #include "AliInputEventHandler.h"
84
85 using std::cout;
86 using std::endl;
87
88 ClassImp(AliAODNuclExReplicator)
89
90 //_____________________________________________________________________________
91 AliAODNuclExReplicator::AliAODNuclExReplicator(const char* name, const char* title,
92                                                // AliAnalysisCuts* trackCut,
93                                                // AliAnalysisCuts* vertexCut,
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 //fVertexCut(vertexCut), 
108   fVertices(0x0), 
109   fNuclei(0x0),
110 //  fTrackCut(trackCut), 
111   fSecondaryVerices(0x0), 
112   fDaughterTracks(0x0),
113   fList(0x0),
114   fMCParticles(0x0),
115   fMCHeader(0x0),
116   fMCMode(mcMode),
117   fLabelMap(),
118   fParticleSelected(),
119   fReplicateHeader(kTRUE), //replicateHeader //to be fixed
120   fnSigmaTrk1(nsigmaTrk1),
121   fnSigmaTrk2(nsigmaTrk2),
122   fpartType1(partType1),
123   fpartType2(partType2),
124   fSecVtxWithKF(kFALSE),
125   fVertexerTracks(0x0),
126   fV1(0x0),
127   fAODMapSize(0),
128   fAODMap(0)
129   
130 {
131   // default ctor
132 }
133
134 //_____________________________________________________________________________
135 AliAODNuclExReplicator::~AliAODNuclExReplicator()
136 {
137   // destructor
138   // delete fTrackCut;
139   // delete fVertexCut;
140   delete fList;
141 }
142
143 //_____________________________________________________________________________
144 void AliAODNuclExReplicator::SelectParticle(Int_t i)
145 {
146   // taking the absolute values here, need to take care 
147   // of negative daughter and mother
148   // IDs when setting!
149   
150   if (!IsParticleSelected(TMath::Abs(i)))
151   {
152     fParticleSelected.Add(TMath::Abs(i),1);    
153   }
154 }
155
156 //_____________________________________________________________________________
157 Bool_t AliAODNuclExReplicator::IsParticleSelected(Int_t i)  
158 {
159   // taking the absolute values here, need to take 
160   // care with negative daughter and mother
161   // IDs when setting!
162   return (fParticleSelected.GetValue(TMath::Abs(i))==1);
163 }
164
165
166 //_____________________________________________________________________________
167 void AliAODNuclExReplicator::CreateLabelMap(const AliAODEvent& source)
168 {  
169   //
170   // this should be called once all selections are done 
171   //
172   
173   fLabelMap.Delete();
174   
175   TClonesArray* mcParticles = static_cast<TClonesArray*>(source.FindListObject(AliAODMCParticle::StdBranchName()));
176   
177   Int_t i(0);
178   Int_t j(0);
179   
180   TIter next(mcParticles);
181   
182   while ( next() ) 
183   {
184     if (IsParticleSelected(i))
185     {
186       fLabelMap.Add(i,j++);
187     }
188     ++i;
189   }
190 }
191
192 //_____________________________________________________________________________
193 Int_t AliAODNuclExReplicator::GetNewLabel(Int_t i) 
194 {
195   // Gets the label from the new created Map
196   // Call CreatLabelMap before
197   // otherwise only 0 returned
198   return fLabelMap.GetValue(TMath::Abs(i));
199 }
200
201 //_____________________________________________________________________________
202 // void AliAODNuclExReplicator::FilterMC(const AliAODEvent& source)
203 // {
204 //   // Filter MC information
205
206 //   fMCHeader->Reset();
207 //   fMCParticles->Clear("C");
208
209 //   AliAODMCHeader* mcHeader(0x0);
210 //   TClonesArray* mcParticles(0x0);
211   
212 //   fParticleSelected.Delete();
213   
214 //   if ( fMCMode>=2 && !fTracks->GetEntries() ) return;
215 //   // for fMCMode==1 we only copy MC information for events where there's at least one muon track
216     
217 //   mcHeader = static_cast<AliAODMCHeader*>(source.FindListObject(AliAODMCHeader::StdBranchName()));
218   
219 //   if ( mcHeader ) 
220 //     {
221 //       *fMCHeader = *mcHeader;
222 //     }
223   
224 //   mcParticles = static_cast<TClonesArray*>(source.FindListObject(AliAODMCParticle::StdBranchName()));
225   
226 //   if ( mcParticles && fMCMode>=2 )
227 //     {
228 //       // loop on (kept) muon tracks to find their ancestors
229 //       TIter nextMT(fTracks);
230 //       AliAODTrack* mt;
231     
232 //       while ( ( mt = static_cast<AliAODTrack*>(nextMT()) ) )
233 //      {
234 //        Int_t label = mt->GetLabel();
235       
236 //        while ( label >= 0 ) 
237 //          {
238 //            SelectParticle(label);
239 //            AliAODMCParticle* mother = static_cast<AliAODMCParticle*>(mcParticles->UncheckedAt(label));
240 //            if (!mother)
241 //              {
242 //                AliError("Got a null mother ! Check that !");
243 //                label = -1;
244 //              }
245 //            else
246 //              {
247 //                label = mother->GetMother();
248 //              }
249 //          }
250 //      }
251     
252 //       CreateLabelMap(source);
253     
254 //       // Actual filtering and label remapping (shamelessly taken for the implementation of AliAODHandler::StoreMCParticles)
255 //       TIter nextMC(mcParticles);
256 //       AliAODMCParticle* p;
257 //       Int_t nmc(0);
258 //       Int_t nmcout(0);
259     
260 //       while ( ( p = static_cast<AliAODMCParticle*>(nextMC()) ) )
261 //      {
262 //        AliAODMCParticle c(*p);
263       
264 //        if ( IsParticleSelected(nmc) )
265 //          {
266 //            // 
267 //            Int_t d0 =  p->GetDaughter(0);
268 //            Int_t d1 =  p->GetDaughter(1);
269 //            Int_t m =   p->GetMother();
270         
271 //            // other than for the track labels, negative values mean
272 //            // no daughter/mother so preserve it
273         
274 //            if(d0<0 && d1<0)
275 //              {
276 //                // no first daughter -> no second daughter
277 //                // nothing to be done
278 //                // second condition not needed just for sanity check at the end
279 //                c.SetDaughter(0,d0);
280 //                c.SetDaughter(1,d1);
281 //              } else if(d1 < 0 && d0 >= 0) 
282 //              {
283 //                // Only one daughter
284 //                // second condition not needed just for sanity check at the end
285 //                if(IsParticleSelected(d0))
286 //                  {
287 //                    c.SetDaughter(0,GetNewLabel(d0));
288 //                  } else 
289 //                  {
290 //                    c.SetDaughter(0,-1);
291 //                  }
292 //                c.SetDaughter(1,d1);
293 //              }
294 //            else if (d0 > 0 && d1 > 0 )
295 //              {
296 //                // we have two or more daughters loop on the stack to see if they are
297 //                // selected
298 //                Int_t d0tmp = -1;
299 //                Int_t d1tmp = -1;
300 //                for (int id = d0; id<=d1;++id)
301 //                  {
302 //                    if (IsParticleSelected(id))
303 //                      {
304 //                        if(d0tmp==-1)
305 //                          {
306 //                            // first time
307 //                            d0tmp = GetNewLabel(id);
308 //                            d1tmp = d0tmp; // this is to have the same schema as on the stack i.e. with one daugther d0 and d1 are the same 
309 //                          }
310 //                        else d1tmp = GetNewLabel(id);
311 //                      }
312 //                  }
313 //                c.SetDaughter(0,d0tmp);
314 //                c.SetDaughter(1,d1tmp);
315 //              } else 
316 //              {
317 //                AliError(Form("Unxpected indices %d %d",d0,d1));
318 //              }
319         
320 //            if ( m < 0 )
321 //              {
322 //                c.SetMother(m);
323 //              } else 
324 //              {
325 //                if (IsParticleSelected(m)) 
326 //                  {
327 //                    c.SetMother(GetNewLabel(m));              
328 //                  }
329 //                else 
330 //                  {
331 //                    AliError(Form("PROBLEM Mother not selected %d", m));              
332 //                  }
333 //              }
334         
335 //            new ((*fMCParticles)[nmcout++]) AliAODMCParticle(c);
336 //          }
337       
338 //        ++nmc;        
339 //      }      
340     
341 //       // now remap the tracks...
342     
343 //       TIter nextTrack(fTracks);
344 //       AliAODTrack* t;
345     
346 //       while ( ( t = static_cast<AliAODTrack*>(nextTrack()) ) )
347 //      {
348 //        t->SetLabel(GetNewLabel(t->GetLabel()));
349 //      }
350     
351 //     }
352 //   else if ( mcParticles ) 
353 //     {
354 //       // simple copy of input MC particles to ouput MC particles
355 //       TIter nextMC(mcParticles);
356 //       AliAODMCParticle* p;
357 //       Int_t nmcout(0);
358     
359 //       while ( ( p = static_cast<AliAODMCParticle*>(nextMC()) ) )
360 //      {
361 //        new ((*fMCParticles)[nmcout++]) AliAODMCParticle(*p);
362 //      }
363 //     }
364   
365 //   AliDebug(1,Form("input mc %d output mc %d",
366 //                   mcParticles ? mcParticles->GetEntries() : 0,
367 //                   fMCParticles ? fMCParticles->GetEntries() : 0));
368   
369 // }
370
371
372 //_____________________________________________________________________________
373 TList* AliAODNuclExReplicator::GetList() const
374 {
375   // return (and build if not already done) our internal list of managed objects
376   if (!fList)
377     {
378     
379       if ( fReplicateHeader )
380         {
381           fHeader = new AliAODHeader;
382         }
383       
384       
385       fSecondaryVerices = new TClonesArray("AliAODRecoDecayLF2Prong",30);
386       fSecondaryVerices->SetName("SecondaryVertices");    
387     
388       fVertices = new TClonesArray("AliAODVertex",2);
389       fVertices->SetName("vertices");    
390     
391       fNuclei = new TClonesArray("AliAODTrack",30);
392       fNuclei->SetName("Nuclei");
393     
394       fDaughterTracks = new TClonesArray("AliAODTrack",30);
395       fDaughterTracks->SetName("DaughterTracks");
396
397     
398       fList = new TList;
399       fList->SetOwner(kTRUE);
400   
401       fList->Add(fHeader);
402       fList->Add(fVertices);
403       fList->Add(fNuclei);
404       fList->Add(fSecondaryVerices);
405       fList->Add(fDaughterTracks);
406   
407       
408       if ( fMCMode > 0 )
409         {
410           fMCHeader = new AliAODMCHeader;    
411           fMCParticles = new TClonesArray("AliAODMCParticle",1000);
412           fMCParticles->SetName(AliAODMCParticle::StdBranchName());
413           fList->Add(fMCHeader);
414           fList->Add(fMCParticles);
415         }
416     }
417   return fList;
418 }
419
420 //_____________________________________________________________________________
421 void AliAODNuclExReplicator::ReplicateAndFilter(const AliAODEvent& source)
422 {
423   // Replicate (and filter if filters are there) the relevant parts we're interested in AODEvent
424   
425   //-----------------------------------------------
426   // AliPIDResponse
427   
428   AliPIDResponse *fPIDResponse;
429   AliAnalysisManager *man=AliAnalysisManager::GetAnalysisManager();
430   AliInputEventHandler* inputHandler = (AliInputEventHandler*) (man->GetInputEventHandler());
431   fPIDResponse = inputHandler->GetPIDResponse();
432   
433   //--------------------------------------------------------
434
435   //  printf("Execute NuclEx Replicator\n");
436
437   //---------------------------------
438
439   if (fReplicateHeader)
440     {
441       *fHeader = *(source.GetHeader());
442     }
443     
444   fVertices->Clear("C");                        
445     
446   fNuclei->Clear("C");
447
448   fSecondaryVerices->Clear("C");
449
450   fDaughterTracks->Clear("C");
451   
452   //----------------------------------
453
454   //  cout<<"Centrality AOD source: "<<source.GetHeader()->GetCentrality()<<endl;
455
456   Int_t nsv(0);
457   Int_t nnuclei(0);
458   Int_t ntracksD(0);
459
460   Int_t input(0);
461   Double_t xdummy,ydummy;
462
463   AliAODRecoDecayLF2Prong *io2Prong  = 0;
464
465   TObjArray *twoTrackArray    = new TObjArray(2);
466   Double_t dispersion;
467
468   // cout<<"Qui"<<endl;
469   //cout<<source.GetMagneticField()<<endl;
470
471   AliAODVertex *vtx = source.GetPrimaryVertex();
472                                                 
473   //  cout<<"Source "<<source<<endl;
474   //cout<<"vtx: "<<vtx<<endl;
475
476   // A Set of very loose cut for a weak strange decay
477   
478   fCosMin       = 0.99;
479   fDCAtracksMin = 1;
480   fRmax         = 200.;
481   fRmin         = 0.1;
482   fDNmin        = 0.05;
483   fDPmin        = 0.05;
484
485   //----------------------------------------------------------
486
487   //  Int_t nindices=0;
488   UShort_t *indices = 0;
489   const Int_t entries = source.GetNumberOfTracks();
490
491   Double_t pos[3],cov[6];
492   vtx->GetXYZ(pos);
493   vtx->GetCovarianceMatrix(cov);
494   fV1 = new AliESDVertex(pos,cov,100.,100,vtx->GetName());
495   //  cout<<"fV1 pricipal loop: "<<fV1<<endl;
496   
497   if(entries<=0) return;
498   indices = new UShort_t[entries];
499   memset(indices,0,sizeof(UShort_t)*entries);
500   fAODMapSize = 100000;
501   fAODMap = new Int_t[fAODMapSize];
502   memset(fAODMap,0,sizeof(Int_t)*fAODMapSize);
503   //  cent=((AliAODEvent&)source)->GetCentrality();
504   
505   //-------------------------------------------------------------
506
507   //AliAODRecoDecayLF   *rd = 0;
508   //  AliAODRecoDecayLF2Prong*rd = 0;
509   
510
511   if(vtx->GetNContributors()<1) {
512     
513     // SPD vertex cut
514     vtx =source.GetPrimaryVertexSPD();
515     
516     if(vtx->GetNContributors()<1) {
517       Info("AliAnalysisTaskHelium3Pi","No good vertex, skip event");
518       return; // NO GOOD VERTEX, SKIP EVENT 
519     }
520   }
521   
522   Double_t xPrimaryVertex=0.,yPrimaryVertex=0.;
523   xPrimaryVertex=vtx->GetX();
524   yPrimaryVertex=vtx->GetY();
525   
526   fBzkG=source.GetMagneticField();
527   fVertexerTracks=new AliVertexerTracks(fBzkG);
528
529   // TF1 *foPion=new TF1("foPion", "[0]*([1]*TMath::Power(TMath::Sqrt(1 + (x/0.13957)*(x/0.13957))/(x/0.13957) , [3]) - 1 - TMath::Power(TMath::Sqrt(1 + (x/0.13957)*(x/0.13957))/(x/0.13957) , [3])*TMath::Log([2] + 1/TMath::Power((x/0.13957), [4])))",0.01,20);
530   // foPion->SetParameters(4.1,8.98482806165147636e+00,1.54000000000000005e-05,2.30445734159456084e+00,2.25624744086878559e+00);
531   
532   Double_t TrackNumber = source.GetNumberOfTracks();
533   
534   //Tracks arrays
535   
536   TArrayI Track0(TrackNumber);        //Pions                                                                          
537   Int_t nTrack0=0;
538   
539   TArrayI Track1(TrackNumber);        //Helium3
540   Int_t nTrack1=0;
541   
542   for(Int_t j=0; j<TrackNumber; j++){
543
544     //    cout<<"Inside loop tracks"<<endl;
545     
546     AliVTrack *track = (AliVTrack*)source.GetTrack(j);
547     
548     AliAODTrack *aodtrack =(AliAODTrack*)track;
549
550     //-----------------------------------------------------------
551     //Track cuts 
552     if(aodtrack->GetTPCNcls() < 70 )continue;
553     if(aodtrack->Chi2perNDF() > 4 )continue;
554
555     // cout<<"Filter Map: "<<aodtrack->GetFilterMap()<<endl;
556     //    if (!aodtrack->TestFilterMask(BIT(4))) continue;
557     // if(aodtrack->TestFilterMask(BIT(4))!=0)
558     //   cout<<"!!!!!!!!!!!!!!!!!!!!!!!!!!! Filter bit mask: !!!!!!!!!!!!!!!!!!!!!!!!!!!"<<endl;
559     
560     if (!aodtrack->IsOn(AliAODTrack::kTPCrefit)) continue;
561     if (!aodtrack->IsOn(AliAODTrack::kTPCin)) continue;
562     if (aodtrack->IsOn(AliAODTrack::kITSpureSA)) continue;
563
564     //---------------------------------------------------------------
565      
566     //    Double_t mom = aodtrack->GetDetPid()->GetTPCmomentum();
567     Double_t mom = aodtrack->P();
568     
569     if(mom<0.150)continue;
570    
571
572  
573     Double_t nSigmaNegPion = TMath::Abs(fPIDResponse->NumberOfSigmasTPC(aodtrack,(AliPID::EParticleType)fpartType1)); 
574     //    cout<<"%%% nsigma Pi: "<<nSigmaNegPion<<endl;
575     //    Double_t nSigmaNegPion = TMath::Abs((aodtrack->GetTPCsignal() - foPion->Eval(mom))/foPion->Eval(mom))/0.07;
576     
577     if ( nSigmaNegPion <= fnSigmaTrk1 /*&& aodtrack->Charge()==-1*/){
578       Track0[nTrack0++]=j;
579     }
580     
581     Double_t nSigmaNuclei = fPIDResponse->NumberOfSigmasTPC(aodtrack,(AliPID::EParticleType)fpartType2);  //To be changed
582     //cout<<"%%% nsigma Nuclei: "<<nSigmaNuclei<<endl;
583   
584     if(nSigmaNuclei>-fnSigmaTrk2 && aodtrack->GetTPCsignal()<1000 && mom>0.2){ 
585    
586     //Double_t triggerDeDx = 4*AliExternalTrackParam::BetheBlochAleph((mom*2)/(0.938*3),1.0288,31.9806,5.04114e-11,2.13096,2.38541);
587       //if(aodtrack->GetTPCsignal() > triggerDeDx && aodtrack->GetTPCsignal()<5000 && mom>0.2 /*&& aodtrack->Charge()==1*/){
588       
589       Track1[nTrack1++]=j;
590       new((*fNuclei)[nnuclei++]) AliAODTrack(*aodtrack);
591       //cout<<"@@@ n sigma NUCLEI"<<nSigmaNuclei<<endl;
592     }
593     
594   }
595   
596   //  cout<<"Delete AodTrack..."<<endl;
597   //  delete aodtrack;
598   //  aodtrack->Delete();
599
600   //Set Track Daughters
601   Track0.Set(nTrack0);
602   Track1.Set(nTrack1);
603   
604   //  cout<<"Track loops..."<<endl;
605   
606   AliAODTrack *postrackAOD = 0;
607   AliAODTrack *negtrackAOD = 0;
608  
609   AliESDtrack *postrack = 0;
610   AliESDtrack *negtrack = 0;
611
612   // cout <<"nTrack1 "<< nTrack1<<endl;
613   // cout <<"nTrack0 "<< nTrack0<<endl;
614
615   for (Int_t i=0; i<nTrack1; i++){                            //! He Tracks Loop
616
617     //    cout<<"Inside track loop"<<endl;
618     
619     Int_t Track1idx=Track1[i];
620
621     AliVTrack *trackpos = (AliVTrack*)source.GetTrack(Track1idx);
622
623     postrackAOD = (AliAODTrack*)trackpos;
624     postrack = new AliESDtrack(trackpos);
625     
626     for (Int_t k=0; k <nTrack0 ; k++) {                           //! Pion Tracks Loop
627
628       Int_t Track0idx=Track0[k];
629
630       AliVTrack *trackneg = (AliVTrack*)source.GetTrack(Track0idx);
631       negtrackAOD =(AliAODTrack*)trackneg;
632       negtrack = new AliESDtrack(trackneg);       
633    
634       twoTrackArray->AddAt(negtrack,0);
635       twoTrackArray->AddAt(postrack,1);
636       
637       Double_t dcap1n1 = postrack->GetDCA(negtrack,fBzkG,xdummy,ydummy);
638
639       Double_t dcap1toPV = TMath::Abs(postrack->GetD(xPrimaryVertex, yPrimaryVertex,fBzkG));
640       Double_t dcan1toPV = TMath::Abs(negtrack->GetD(xPrimaryVertex, yPrimaryVertex,fBzkG));
641
642       if(dcap1n1>fDCAtracksMin)continue;
643       if((xdummy+ydummy)>fRmax )continue;
644       if((xdummy+ydummy)< fRmin)continue;
645       
646       if ( dcan1toPV < fDNmin)               
647         if ( dcap1toPV < fDPmin) continue;   
648       
649       //      cout<<"dcap1n1: "<<dcap1n1<<endl;
650  
651       AliAODVertex *vertexp1n1 = ReconstructSecondaryVertex(twoTrackArray,dispersion,kTRUE);
652       
653       if(!vertexp1n1) {
654
655         twoTrackArray->Clear();
656         delete negtrack;
657         negtrack=NULL; 
658         continue; 
659       }
660       
661       //      cout<<"vertexp1n1:  "<<vertexp1n1<<endl;
662
663       io2Prong = Make2Prong(twoTrackArray,source,vertexp1n1,dcap1n1);
664   
665       //cout<<"Fuori dal loop e' quello che salvo \npr0: "<<io2Prong->GetProngID(0)<<" pr1: "<<io2Prong->GetProngID(1)<<" pr2"<<io2Prong->GetProngID(2)<<endl;
666
667       if(io2Prong->CosPointingAngle()<fCosMin)continue;
668
669       //      cout<<"if CPA \npr0: "<<io2Prong->GetProngID(0)<<" pr1: "<<io2Prong->GetProngID(1)<<" pr2"<<io2Prong->GetProngID(2)<<endl;
670       // cout<<"pointing angle "<<io2Prong->CosPointingAngle()<<endl;
671       
672       
673       // AliAODTrack *trk0 = (AliAODTrack*)io2Prong->GetDaughter(0);
674       // AliAODTrack *trk1 = (AliAODTrack*)io2Prong->GetDaughter(1);
675     
676       // cout<<"**********************************************"<<endl;
677       // cout<<trk0/*->GetID()*/<<" "<<negtrackAOD->GetID()<<endl;
678       // cout<<trk1/*->GetID()*/<<" "<<postrackAOD->GetID()<<endl;
679       // cout<<"d0 io2Prong: "<<io2Prong->GetProngID(1)<<endl;
680       // cout<<"d1 io2Prong: "<<io2Prong->GetProngID(0)<<endl;
681       // cout<<"**********************************************"<<endl;
682
683       //      rd =  new((*fSecondaryVerices)[nsv++]) AliAODRecoDecayLF2Prong(*io2Prong);
684       
685       new((*fSecondaryVerices)[nsv++]) AliAODRecoDecayLF2Prong(*io2Prong);
686       
687       //      cout<<"QUELLO CHE SALVo \npr0: "<<rd->GetProngID(0)<<" pr1: "<<rd->GetProngID(1)<<" pr2"<<rd->GetProngID(2)<<endl;
688
689       // rd->SetSecondaryVtx(vertexp1n1);
690       // vertexp1n1->SetParent(rd);
691       // AddRefs(vertexp1n1,rd,source,twoTrackArray);
692       
693         
694       new((*fDaughterTracks)[ntracksD++]) AliAODTrack(*negtrackAOD);
695       new((*fDaughterTracks)[ntracksD++]) AliAODTrack(*postrackAOD);
696       
697       delete negtrack; 
698       negtrack = NULL;
699     
700       delete vertexp1n1;
701       vertexp1n1= NULL;
702       continue;
703     }
704
705     delete postrack; 
706     postrack = NULL;
707    
708   }
709   
710   //----------------------------------------------------------
711  
712   assert(fVertices!=0x0);
713   fVertices->Clear("C");
714   TIter nextV(source.GetVertices());
715   AliAODVertex* v;
716   Int_t nvertices(0);
717   
718   while ( ( v = static_cast<AliAODVertex*>(nextV()) ) )
719     {
720        if (  v->GetType() == AliAODVertex::kPrimary     ||
721             v->GetType() == AliAODVertex::kMainSPD     ||
722             v->GetType() == AliAODVertex::kPileupSPD   ||
723             v->GetType() == AliAODVertex::kPileupTracks||
724             v->GetType() == AliAODVertex::kMainTPC  ) 
725         {
726           AliAODVertex* tmp = v->CloneWithoutRefs();
727           AliAODVertex* copiedVertex = new((*fVertices)[nvertices++]) AliAODVertex(*tmp);
728           
729           // to insure the main vertex retains the ncontributors information
730           // (which is otherwise computed dynamically from
731           // references to tracks, which we do not keep in muon aods...)
732           // we set it here
733           
734           copiedVertex->SetNContributors(v->GetNContributors()); 
735           
736           //  fVertices->Delete();
737           //      delete copiedVertex;
738           delete tmp;
739         }
740     }
741   
742   //  printf("....Done NuclEx Replicator...\n");
743   
744   AliDebug(1,Form("input mu tracks=%d tracks=%d vertices=%d nnuclei=%d",
745                   input,fSecondaryVerices->GetEntries(),fVertices->GetEntries(),fNuclei->GetEntries())); 
746   
747   // Finally, deal with MC information, if needed
748   
749   // if ( fMCMode > 0 )
750   //   {
751   //     FilterMC(source);      
752   //   }
753   
754   //cout<<"Delete..."<<endl;
755   // delete foPion;
756   // foPion = NULL;
757   //cout<<"Delete 1"<<  endl;
758   
759   if(io2Prong) {delete io2Prong; io2Prong=NULL;}
760   //cout<<"Delete 2"<<  endl;
761   twoTrackArray->Delete();  delete twoTrackArray;
762   //cout<<"Delete 3"<<  endl;
763   // vtx->Delete();  delete vtx;
764   //cout<<"Delete 4"<<  endl;
765   if(fV1) { delete fV1; fV1=NULL; }
766   //cout<<"Delete 5"<<  endl;
767   if(fAODMap) { delete [] fAODMap; fAODMap=NULL; }
768   delete []indices;
769   //cout<<"Delete 6"<<  endl;
770   delete fVertexerTracks;
771
772   
773
774
775 //-----------------------------------------------------------------------------
776
777 AliAODVertex *AliAODNuclExReplicator::ReconstructSecondaryVertex(TObjArray *trkArray,
778                                                                 Double_t &dispersion,Bool_t useTRefArray) const
779
780 {
781   // Secondary vertex reconstruction with AliVertexerTracks or AliKFParticle
782   //AliCodeTimerAuto("",0);
783
784   AliESDVertex *vertexESD = 0;
785   AliAODVertex *vertexAOD = 0;
786
787   if(!fSecVtxWithKF) { // AliVertexerTracks
788
789     fVertexerTracks->SetVtxStart(fV1);
790     vertexESD = (AliESDVertex*)fVertexerTracks->VertexForSelectedESDTracks(trkArray);
791
792     if(!vertexESD) return vertexAOD;
793
794
795     Double_t vertRadius2=vertexESD->GetX()*vertexESD->GetX()+vertexESD->GetY()*vertexESD->GetY();
796     if(vertRadius2>200.){
797       // vertex outside beam pipe, reject candidate to avoid propagation through material
798       delete vertexESD; vertexESD=NULL;
799       return vertexAOD;
800     }
801
802   } else { // Kalman Filter vertexer (AliKFParticle)
803
804     AliKFParticle::SetField(fBzkG);
805
806     AliKFVertex vertexKF;
807
808     Int_t nTrks = trkArray->GetEntriesFast();
809     for(Int_t i=0; i<nTrks; i++) {
810       AliESDtrack *esdTrack = (AliESDtrack*)trkArray->At(i);
811       AliKFParticle daughterKF(*esdTrack,211);
812       vertexKF.AddDaughter(daughterKF);
813     }
814     vertexESD = new AliESDVertex(vertexKF.Parameters(),
815                                  vertexKF.CovarianceMatrix(),
816                                  vertexKF.GetChi2(),
817                                  vertexKF.GetNContributors());
818
819   }
820   // convert to AliAODVertex
821   Double_t pos[3],cov[6],chi2perNDF;
822   vertexESD->GetXYZ(pos); // position
823   vertexESD->GetCovMatrix(cov); //covariance matrix
824   chi2perNDF = vertexESD->GetChi2toNDF();
825   dispersion = vertexESD->GetDispersion();
826   delete vertexESD; 
827   vertexESD=NULL;
828
829   Int_t nprongs= (useTRefArray ? 0 : trkArray->GetEntriesFast());
830   vertexAOD = new AliAODVertex(pos,cov,chi2perNDF,0x0,-1,AliAODVertex::kUndef,nprongs);
831
832   //cout<<"------------------> Reconstruct vertexAOD: "<<vertexAOD<<endl;
833
834   return vertexAOD;
835
836  
837 }
838
839 //-----------------------------------------------------------------------------
840
841 AliAODRecoDecayLF2Prong* AliAODNuclExReplicator::Make2Prong(TObjArray *twoTrackArray,const AliAODEvent &evento,
842                                                            AliAODVertex *secVert,Double_t dca)
843
844
845
846
847   //cout<<"Inside make 2 prong"<<endl;
848
849   Int_t charge[2];
850   charge[0]=1; //it was -1 //Ramona
851   charge[1]=2;
852       
853   // From AliAnalysisVertexingLF.cxx Method:Make2Prongs
854   
855   fBzkG = evento.GetMagneticField();
856       
857   Double_t px[2],py[2],pz[2],d0[2],d0err[2];
858   // Also this has been changed //Ramona
859   AliESDtrack *negtrack = (AliESDtrack*)twoTrackArray->UncheckedAt(0);
860   AliESDtrack *postrack = (AliESDtrack*)twoTrackArray->UncheckedAt(1);
861
862   //cout<<"negtrack: "<<negtrack<<" postrack: "<<postrack<<endl;
863   //cout<<"kVeryBig: "<<kVeryBig<<endl;
864   //cout<<"secVert: "<<secVert<<endl;
865
866   // // propagate tracks to secondary vertex, to compute inv. mass
867   
868   postrack->PropagateToDCA(secVert,fBzkG,kVeryBig);
869   negtrack->PropagateToDCA(secVert,fBzkG,kVeryBig);
870   
871   Double_t momentum[3];
872   
873   negtrack->GetPxPyPz(momentum);
874   px[0] = charge[0]*momentum[0]; 
875   py[0] = charge[0]*momentum[1]; 
876   pz[0] = charge[0]*momentum[2]; 
877
878   postrack->GetPxPyPz(momentum);
879   px[1] = charge[1]*momentum[0]; 
880   py[1] = charge[1]*momentum[1]; 
881   pz[1] = charge[1]*momentum[2]; 
882   
883   //cout<< px[0] <<" "<< " "<< py[0] << " "<< pz[0]<<endl;
884   //  px[0] = momentum[0]; py[0] = momentum[1]; pz[0] = momentum[2]; 
885   //cout<< px[1] <<" "<< " "<< py[1] << " "<< pz[1]<<endl;
886   //px[1] = momentum[0]; py[1] = momentum[1]; pz[1] = momentum[2]; 
887   
888
889   // primary vertex to be used by this candidate
890   AliAODVertex *primVertexAOD  = evento.GetPrimaryVertex();
891   //cout<<"primVertexAOD "<<primVertexAOD<<endl;
892   if(!primVertexAOD) return 0x0;
893       
894   Double_t d0z0[2],covd0z0[3];
895
896   d0z0[0] = -999.;
897   d0z0[1] = -999.;
898   covd0z0[0] = -999.;
899   covd0z0[1] = -999.;
900   covd0z0[2] = -999.;
901
902   d0[0] = d0z0[0];
903   d0err[0] = TMath::Sqrt(TMath::Abs(covd0z0[0]));
904
905   d0[1] = d0z0[0];
906   d0err[1] = TMath::Sqrt(TMath::Abs(covd0z0[0]));
907   
908   // create the object AliAODRecoDecayLF2Prong
909   //  AliAODRecoDecayLF2Prong *the2Prong = new AliAODRecoDecayLF2Prong(secVert,px,py,pz,d0,d0err,dcap1n1);
910   AliAODRecoDecayLF2Prong *the2Prong = new AliAODRecoDecayLF2Prong(secVert,px,py,pz,d0,d0err,dca);
911   the2Prong->SetOwnPrimaryVtx(primVertexAOD);
912   
913   //  the2Prong->SetProngIDs(2,{-999,999});
914   // UShort_t id0[2]={99999,999999};
915   // the2Prong->SetProngIDs(2,id0);
916
917   UShort_t id[2]={(UShort_t)postrack->GetID(),(UShort_t)negtrack->GetID()};
918   the2Prong->SetProngIDs(2,id);
919
920   // cout<<"\n\n\nMake 2 Prong: id[0]"<<id[0]<<" id[1]: "<<id[1]<<endl;
921   // cout<<"Get: 1 "<<the2Prong->GetProngID(0)<<" 2 "<<the2Prong->GetProngID(1)<<endl;
922   // cout<<"Get: 3 "<<the2Prong->GetProngID(2)<<" 4 "<<the2Prong->GetProngID(3)<<endl<<endl<<endl;
923   //delete primVertexAOD; primVertexAOD=NULL;
924   
925   if(postrack->Charge()!=0 && negtrack->Charge()!=0) {
926       
927     AddDaughterRefs(secVert,(AliAODEvent&)evento,twoTrackArray);
928     //    AddDaughterRefs(secVert,(AliAODEvent*)evento,twoTrackArray);
929       
930   }
931   
932   return the2Prong;  
933
934   delete the2Prong;
935 }
936
937 //----------------------------------------------------------------------------
938 void AliAODNuclExReplicator::AddDaughterRefs(AliAODVertex *v,
939                                             const AliAODEvent &event,
940                                             const TObjArray *trkArray) const
941
942 {
943   // Add the AOD tracks as daughters of the vertex (TRef)
944   // AliCodeTimerAuto("",0);
945   // cout<<"Inside  AddDaughterRefs "<<endl;
946
947   Int_t nDg = v->GetNDaughters();
948   
949   TObject *dg = 0;
950   if(nDg) dg = v->GetDaughter(0);
951   if(dg) return; // daughters already added
952   
953   Int_t nTrks = trkArray->GetEntriesFast();
954   
955   AliExternalTrackParam *track = 0;
956   AliAODTrack *aodTrack = 0;
957   Int_t id;
958   
959   for(Int_t i=0; i<nTrks; i++) {
960     track = (AliExternalTrackParam*)trkArray->UncheckedAt(i);
961   
962     id = (Int_t)track->GetID();
963     //    printf("---> %d\n",id);
964     if(id<0) continue; // this track is a AliAODRecoDecay
965   
966     aodTrack = (AliAODTrack*)event.GetTrack(fAODMap[id]);
967     v->AddDaughter(aodTrack);
968   }
969   return;
970 }
971 //----------------------------------------------------------------------------
972         
973 void AliAODNuclExReplicator::AddRefs(AliAODVertex *v,AliAODRecoDecayLF *rd,
974                                     const AliAODEvent &event,
975                                     const TObjArray *trkArray) const
976
977 {
978   // Add the AOD tracks as daughters of the vertex (TRef)
979   // Also add the references to the primary vertex and to the cuts
980   //AliCodeTimerAuto("",0);
981   
982   AddDaughterRefs(v,event,trkArray);
983   rd->SetPrimaryVtxRef((AliAODVertex*)event.GetPrimaryVertex());
984   return;
985 }       
986
987 //---------------------------------------------------------------------------
988
989 void AliAODNuclExReplicator::Terminate(){
990
991 }