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       AliAODHeader * header = dynamic_cast<AliAODHeader*>(source.GetHeader());
 270       if(!header) AliFatal("Not a standard AOD");
 271       *fHeader = *(header);
 272     }
 273
 274   fVertices->Clear("C");
 275
 276   fNuclei->Clear("C");
 277
 278   fSecondaryVerices->Clear("C");
 279
 280   fDaughterTracks->Clear("C");
 281
 282   //----------------------------------
 283
 284   //retrive MC infos
 285
 286   TClonesArray *arrayMC = 0;
 287   AliAODMCHeader *mcHeader=0;
 288   Int_t mumpdg=-100;
 289   //  Int_t mumpdgNeg=-100;
 290
 291
 292   arrayMC = (TClonesArray*) source.GetList()->FindObject(AliAODMCParticle::StdBranchName());
 293   if (!arrayMC) {
 294     Printf("Error: MC particles branch not found!\n");
 295     return;
 296   }
 297
 298   // if(arrayMC)
 299   //   cout<<"Ho caricato array mc"<<endl;
 300
 301   mcHeader =  (AliAODMCHeader*)source.GetList()->FindObject(AliAODMCHeader::StdBranchName());
 302   if(!mcHeader) {
 303     printf("AliAnalysisTaskSEDs::UserExec: MC header branch not found!\n");
 304     return;
 305   }
 306
 307   // if(mcHeader)
 308   //   cout<<"Ho caricato MC header"<<endl;
 309
 310   //  cout<<"Centrality AOD source: "<<source.GetHeader()->GetCentrality()<<endl;
 311
 312   Int_t nsv(0);
 313   //  Int_t nnuclei(0);
 314   Int_t ntracksD(0);
 315
 316   Int_t input(0);
 317   Double_t xdummy,ydummy;
 318
 319   AliAODRecoDecayLF2Prong *io2Prong  = 0;
 320
 321   TObjArray *twoTrackArray    = new TObjArray(2);
 322   Double_t dispersion;
 323
 324   // cout<<"Qui"<<endl;
 325   // cout<<source.GetMagneticField()<<endl;
 326
 327   AliAODVertex *vtx = source.GetPrimaryVertex();
 328
 329   // cout<<"Source "<<source<<endl;
 330   // cout<<"vtx: "<<vtx<<endl;
 331
 332   // A Set of very loose cut for a weak strange decay
 333
 334   fCosMin       = 0.99;
 335   fDCAtracksMin = 1;
 336   fRmax         = 200.;
 337   fRmin         = 0.1;
 338   fDNmin        = 0.05;
 339   fDPmin        = 0.05;
 340
 341   //----------------------------------------------------------
 342
 343   //  Int_t nindices=0;
 344   UShort_t *indices = 0;
 345   const Int_t entries = source.GetNumberOfTracks();
 346
 347   Double_t pos[3],cov[6];
 348   vtx->GetXYZ(pos);
 349   vtx->GetCovarianceMatrix(cov);
 350   fV1 = new AliESDVertex(pos,cov,100.,100,vtx->GetName());
 351   cout<<"fV1 pricipal loop: "<<fV1<<endl;
 352
 353   if(entries<=0) return;
 354   indices = new UShort_t[entries];
 355   memset(indices,0,sizeof(UShort_t)*entries);
 356   fAODMapSize = 100000;
 357   fAODMap = new Int_t[fAODMapSize];
 358   memset(fAODMap,0,sizeof(Int_t)*fAODMapSize);
 359   //  cent=((AliAODEvent&)source)->GetCentrality();
 360
 361   //-------------------------------------------------------------
 362
 363   if(vtx->GetNContributors()<1) {
 364
 365     // SPD vertex cut
 366     vtx =source.GetPrimaryVertexSPD();
 367
 368     if(vtx->GetNContributors()<1) {
 369       Info("AliAnalysisTaskHelium3Pi","No good vertex, skip event");
 370       return; // NO GOOD VERTEX, SKIP EVENT
 371     }
 372   }
 373
 374   Double_t xPrimaryVertex=0.,yPrimaryVertex=0.;
 375   xPrimaryVertex=vtx->GetX();
 376   yPrimaryVertex=vtx->GetY();
 377
 378   fBzkG=source.GetMagneticField();
 379   fVertexerTracks=new AliVertexerTracks(fBzkG);
 380
 381   Double_t TrackNumber = source.GetNumberOfTracks();
 382   Int_t label =-1;
 383
 384   //Tracks arrays
 385
 386   TArrayI Track0(TrackNumber);        //Pions
 387   Int_t nTrack0=0;
 388
 389   TArrayI Track1(TrackNumber);        //Helium3
 390   Int_t nTrack1=0;
 391
 392   for(Int_t j=0; j<TrackNumber; j++){
 393
 394     //    cout<<"Inside loop tracks"<<endl;
 395
 396
 397     AliVTrack *track = (AliVTrack*)source.GetTrack(j);
 398
 399     AliAODTrack *aodtrack =(AliAODTrack*)track;
 400
 401     //-----------------------------------------------------------
 402     //Track cuts
 403     if(aodtrack->GetTPCNcls() < 70 )continue;
 404     if(aodtrack->Chi2perNDF() > 4 )continue;
 405
 406     if (!aodtrack->IsOn(AliAODTrack::kTPCrefit)) continue;
 407     if (!aodtrack->IsOn(AliAODTrack::kTPCin)) continue;
 408     if (aodtrack->IsOn(AliAODTrack::kITSpureSA)) continue;
 409
 410     //---------------------------------------------------------------
 411
 412     Double_t mom = aodtrack->P();
 413
 414     if(mom<0.150)continue;
 415
 416     label = TMath::Abs(aodtrack->GetLabel());
 417
 418     AliAODMCParticle *part = (AliAODMCParticle*) arrayMC->At(label);
 419
 420     Int_t PDGCode=part->GetPdgCode();
 421
 422     Int_t mumid = part->GetMother();
 423
 424     if(mumid>-1){
 425       AliAODMCParticle *mother = (AliAODMCParticle*) arrayMC->At(mumid);
 426       mumpdg = mother->GetPdgCode();
 427     }
 428
 429     //    if(mumpdg == 1010010030 ||mumpdg == -1010010030 ){
 430
 431     if(mumpdg == -1010010030){
 432
 433       if(PDGCode==-211 || PDGCode==+211){
 434         //      if(PDGCode==-211){
 435
 436         //      if(aodtrack->Charge()>0)continue;
 437
 438         Track0[nTrack0++]=j;
 439       }
 440
 441       if(PDGCode==1000020030 ||PDGCode==-1000020030 ){
 442       //if(PDGCode==1000020030){
 443
 444         //      if(aodtrack->Charge()<0)continue;
 445
 446         Track1[nTrack1++]=j;
 447         //      new((*fNuclei)[nnuclei++]) AliAODTrack(*aodtrack);
 448       }
 449     }
 450   }
 451
 452   //Set Track Daughters
 453   Track0.Set(nTrack0);
 454   Track1.Set(nTrack1);
 455
 456
 457   // cout<<"Track loops..."<<endl;
 458   // cout<<"npos "<<nTrack1<<endl;
 459   // cout<<"nneg "<<nTrack0<<endl;
 460
 461
 462   AliAODTrack *postrackAOD = 0;
 463   AliAODTrack *negtrackAOD = 0;
 464
 465   AliESDtrack *postrack = 0;
 466   AliESDtrack *negtrack = 0;
 467
 468   Bool_t isOk=kFALSE;
 469
 470   for (Int_t i=0; i<nTrack1; i++){                            //! He Tracks Loop
 471
 472     Int_t Track1idx=Track1[i];
 473
 474     AliVTrack *trackpos = (AliVTrack*)source.GetTrack(Track1idx);
 475
 476     postrackAOD = (AliAODTrack*)trackpos;
 477     postrack = new AliESDtrack(trackpos);
 478
 479     //--- MC infos
 480
 481     Int_t labelpos = TMath::Abs(postrack->GetLabel());
 482     AliAODMCParticle *partPos = (AliAODMCParticle*) arrayMC->At(labelpos);
 483     Int_t mumidPos = partPos->GetMother();
 484
 485     //    cout<<"\n\nmumidPos: "<<mumidPos <<endl;
 486
 487     //------------------------------
 488
 489     for (Int_t k=0; k <nTrack0 ; k++) {                           //! Pion Tracks Loop
 490
 491       Int_t Track0idx=Track0[k];
 492
 493       AliVTrack *trackneg = (AliVTrack*)source.GetTrack(Track0idx);
 494       negtrackAOD =(AliAODTrack*)trackneg;
 495       negtrack = new AliESDtrack(trackneg);
 496
 497       //--- MC infos
 498
 499       Int_t labelneg = TMath::Abs(negtrack->GetLabel());
 500       AliAODMCParticle *partNeg = (AliAODMCParticle*) arrayMC->At(labelneg);
 501       Int_t mumidNeg = partNeg->GetMother();
 502
 503       //      cout<<"mumidNeg: "<<mumidNeg <<endl;
 504
 505       // AliAODMCParticle *motherNeg = (AliAODMCParticle*) arrayMC->At(mumidNeg);
 506       // mumpdgNeg = motherNeg->GetPdgCode();
 507
 508       //------------------------------
 509       //  if(mumidPos == mumidNeg && mumidNeg > 0){
 510       isOk=kFALSE;
 511
 512       // if(mumidPos == mumidNeg && mumidNeg > 0 && mumpdgNeg  == 1010010030)
 513       if(mumidPos == mumidNeg && mumidNeg > 0)
 514         isOk = kTRUE;
 515
 516       twoTrackArray->AddAt(negtrack,0);
 517       twoTrackArray->AddAt(postrack,1);
 518
 519       Double_t dcap1n1 = postrack->GetDCA(negtrack,fBzkG,xdummy,ydummy);
 520
 521       Double_t dcap1toPV = TMath::Abs(postrack->GetD(xPrimaryVertex, yPrimaryVertex,fBzkG));
 522       Double_t dcan1toPV = TMath::Abs(negtrack->GetD(xPrimaryVertex, yPrimaryVertex,fBzkG));
 523
 524       if(dcap1n1>fDCAtracksMin)continue;
 525       if((xdummy+ydummy)>fRmax )continue;
 526       if((xdummy+ydummy)< fRmin)continue;
 527
 528       if ( dcan1toPV < fDNmin)
 529         if ( dcap1toPV < fDPmin) continue;
 530
 531       //      cout<<"dcap1n1: "<<dcap1n1<<endl;
 532
 533       AliAODVertex *vertexp1n1 = ReconstructSecondaryVertex(twoTrackArray,dispersion,kTRUE);
 534
 535       if(!vertexp1n1) {
 536
 537         twoTrackArray->Clear();
 538         delete negtrack;
 539         negtrack=NULL;
 540         continue;
 541       }
 542
 543       io2Prong = Make2Prong(twoTrackArray,source,vertexp1n1,dcap1n1);
 544
 545       if(io2Prong->CosPointingAngle()<fCosMin)continue;
 546
 547       //      cout<<"CPA: "<<io2Prong->CosPointingAngle()<<endl;
 548
 549       // AliAODTrack *trk0 = (AliAODTrack*)io2Prong->GetDaughter(0);
 550       // AliAODTrack *trk1 = (AliAODTrack*)io2Prong->GetDaughter(1);
 551
 552       // cout<<"**********************************************"<<endl;
 553       // cout<<trk0/*->GetID()*/<<" "<<negtrackAOD->GetID()<<endl;
 554       // cout<<trk1/*->GetID()*/<<" "<<postrackAOD->GetID()<<endl;
 555       // cout<<"d0 io2Prong: "<<io2Prong->GetProngID(1)<<endl;
 556       // cout<<"d1 io2Prong: "<<io2Prong->GetProngID(0)<<endl;
 557       // cout<<"**********************************************"<<endl;
 558
 559       //      rd =  new((*fSecondaryVerices)[nsv++]) AliAODRecoDecayLF2Prong(*io2Prong);
 560
 561       //      cout<<"is ok??? "<<isOk<<endl;
 562       if(isOk){
 563         new((*fSecondaryVerices)[nsv++]) AliAODRecoDecayLF2Prong(*io2Prong);
 564         new((*fDaughterTracks)[ntracksD++]) AliAODTrack(*negtrackAOD);
 565         new((*fDaughterTracks)[ntracksD++]) AliAODTrack(*postrackAOD);
 566
 567       }
 568       // rd->SetSecondaryVtx(vertexp1n1);
 569       // vertexp1n1->SetParent(rd);
 570       // AddRefs(vertexp1n1,rd,source,twoTrackArray);
 571
 572       delete negtrack;
 573       negtrack = NULL;
 574
 575       delete vertexp1n1;
 576       vertexp1n1= NULL;
 577       continue;
 578       //   }
 579     }
 580
 581     delete postrack;
 582     postrack = NULL;
 583
 584   }
 585
 586   //----------------------------------------------------------
 587
 588   assert(fVertices!=0x0);
 589   fVertices->Clear("C");
 590   TIter nextV(source.GetVertices());
 591   AliAODVertex* v;
 592   Int_t nvertices(0);
 593
 594   while ( ( v = static_cast<AliAODVertex*>(nextV()) ) )
 595     {
 596       if ( v->GetType() == AliAODVertex::kPrimary     ||
 597            v->GetType() == AliAODVertex::kMainSPD     ||
 598            v->GetType() == AliAODVertex::kPileupSPD   ||
 599            v->GetType() == AliAODVertex::kPileupTracks||
 600            v->GetType() == AliAODVertex::kMainTPC  )
 601         {
 602
 603           AliAODVertex* tmp = v->CloneWithoutRefs();
 604           //AliAODVertex* copiedVertex = new((*fVertices)[nvertices++]) AliAODVertex(*tmp);
 605           new((*fVertices)[nvertices++]) AliAODVertex(*tmp);
 606
 607           // to insure the main vertex retains the ncontributors information
 608           // (which is otherwise computed dynamically from
 609           // references to tracks, which we do not keep in muon aods...)
 610           // we set it here
 611
 612           //copiedVertex->SetNContributors(v->GetNContributors());
 613
 614           //  fVertices->Delete();
 615           //      delete copiedVertex;
 616           delete tmp;
 617           //      printf("....Prendo il vertice primario...\n");
 618         }
 619       //       printf("....Prendo il vertice primario...\n");
 620     }
 621
 622   //printf("....Done NuclEx Replicator...\n");
 623
 624   AliDebug(1,Form("input mu tracks=%d tracks=%d vertices=%d nnuclei=%d",
 625                   input,fSecondaryVerices->GetEntries(),fVertices->GetEntries(),fNuclei->GetEntries()));
 626
 627   // cout<<"Delete..."<<endl;
 628   // delete foPion;
 629   // foPion = NULL;
 630   //cout<<"Delete 1"<<  endl;
 631
 632   if(io2Prong) {delete io2Prong; io2Prong=NULL;}
 633   //cout<<"Delete 2"<<  endl;
 634   twoTrackArray->Delete();  delete twoTrackArray;
 635   //  cout<<"Delete 3"<<  endl;
 636   // vtx->Delete();  delete vtx;
 637   //  cout<<"Delete 4"<<  endl;
 638   if(fV1) { delete fV1; fV1=NULL; }
 639   //  cout<<"Delete 5"<<  endl;
 640   if(fAODMap) { delete [] fAODMap; fAODMap=NULL; }
 641   delete []indices;
 642   //  cout<<"Delete 6"<<  endl;
 643   delete fVertexerTracks;
 644
 645   //  cout<<"Mi rompo alla fine. Perche???"<<endl;
 646
 647 }
 648
 649 //-----------------------------------------------------------------------------
 650
 651 AliAODVertex *AliAODMCNuclExReplicator::ReconstructSecondaryVertex(TObjArray *trkArray,
 652                                                                 Double_t &dispersion,Bool_t useTRefArray) const
 653
 654 {
 655   // Secondary vertex reconstruction with AliVertexerTracks or AliKFParticle
 656   //AliCodeTimerAuto("",0);
 657
 658   AliESDVertex *vertexESD = 0;
 659   AliAODVertex *vertexAOD = 0;
 660
 661   if(!fSecVtxWithKF) { // AliVertexerTracks
 662
 663     fVertexerTracks->SetVtxStart(fV1);
 664     vertexESD = (AliESDVertex*)fVertexerTracks->VertexForSelectedESDTracks(trkArray);
 665
 666     if(!vertexESD) return vertexAOD;
 667
 668
 669     Double_t vertRadius2=vertexESD->GetX()*vertexESD->GetX()+vertexESD->GetY()*vertexESD->GetY();
 670     if(vertRadius2>200.){
 671       // vertex outside beam pipe, reject candidate to avoid propagation through material
 672       delete vertexESD; vertexESD=NULL;
 673       return vertexAOD;
 674     }
 675
 676   } else { // Kalman Filter vertexer (AliKFParticle)
 677
 678     AliKFParticle::SetField(fBzkG);
 679
 680     AliKFVertex vertexKF;
 681
 682     Int_t nTrks = trkArray->GetEntriesFast();
 683     for(Int_t i=0; i<nTrks; i++) {
 684       AliESDtrack *esdTrack = (AliESDtrack*)trkArray->At(i);
 685       AliKFParticle daughterKF(*esdTrack,211);
 686       vertexKF.AddDaughter(daughterKF);
 687     }
 688     vertexESD = new AliESDVertex(vertexKF.Parameters(),
 689                                  vertexKF.CovarianceMatrix(),
 690                                  vertexKF.GetChi2(),
 691                                  vertexKF.GetNContributors());
 692
 693   }
 694   // convert to AliAODVertex
 695   Double_t pos[3],cov[6],chi2perNDF;
 696   vertexESD->GetXYZ(pos); // position
 697   vertexESD->GetCovMatrix(cov); //covariance matrix
 698   chi2perNDF = vertexESD->GetChi2toNDF();
 699   dispersion = vertexESD->GetDispersion();
 700   delete vertexESD;
 701   vertexESD=NULL;
 702
 703   Int_t nprongs= (useTRefArray ? 0 : trkArray->GetEntriesFast());
 704   vertexAOD = new AliAODVertex(pos,cov,chi2perNDF,0x0,-1,AliAODVertex::kUndef,nprongs);
 705
 706   //  cout<<"------------------> Reconstruct vertexAOD: "<<vertexAOD<<endl;
 707
 708   return vertexAOD;
 709
 710
 711 }
 712
 713 //-----------------------------------------------------------------------------
 714
 715 AliAODRecoDecayLF2Prong* AliAODMCNuclExReplicator::Make2Prong(TObjArray *twoTrackArray,const AliAODEvent &evento,
 716                                                            AliAODVertex *secVert,Double_t dca)
 717
 718
 719 {
 720
 721   //cout<<"Inside make 2 prong"<<endl;
 722
 723   Int_t charge[2];
 724   charge[0]=1; //it was -1 //Ramona
 725   charge[1]=2;
 726
 727   // From AliAnalysisVertexingLF.cxx Method:Make2Prongs
 728
 729   fBzkG = evento.GetMagneticField();
 730
 731   Double_t px[2],py[2],pz[2],d0[2],d0err[2];
 732   // Also this has been changed //Ramona
 733   AliESDtrack *negtrack = (AliESDtrack*)twoTrackArray->UncheckedAt(0);
 734   AliESDtrack *postrack = (AliESDtrack*)twoTrackArray->UncheckedAt(1);
 735
 736   //cout<<"negtrack: "<<negtrack<<" postrack: "<<postrack<<endl;
 737   //cout<<"kVeryBig: "<<kVeryBig<<endl;
 738   //cout<<"secVert: "<<secVert<<endl;
 739
 740   // // propagate tracks to secondary vertex, to compute inv. mass
 741
 742   postrack->PropagateToDCA(secVert,fBzkG,kVeryBig);
 743   negtrack->PropagateToDCA(secVert,fBzkG,kVeryBig);
 744
 745   Double_t momentum[3];
 746
 747   negtrack->GetPxPyPz(momentum);
 748   px[0] = charge[0]*momentum[0];
 749   py[0] = charge[0]*momentum[1];
 750   pz[0] = charge[0]*momentum[2];
 751
 752   postrack->GetPxPyPz(momentum);
 753   px[1] = charge[1]*momentum[0];
 754   py[1] = charge[1]*momentum[1];
 755   pz[1] = charge[1]*momentum[2];
 756
 757   //cout<< px[0] <<" "<< " "<< py[0] << " "<< pz[0]<<endl;
 758   //  px[0] = momentum[0]; py[0] = momentum[1]; pz[0] = momentum[2];
 759   //cout<< px[1] <<" "<< " "<< py[1] << " "<< pz[1]<<endl;
 760   //px[1] = momentum[0]; py[1] = momentum[1]; pz[1] = momentum[2];
 761
 762
 763   // primary vertex to be used by this candidate
 764   AliAODVertex *primVertexAOD  = evento.GetPrimaryVertex();
 765   //cout<<"primVertexAOD "<<primVertexAOD<<endl;
 766   if(!primVertexAOD) return 0x0;
 767
 768   Double_t d0z0[2],covd0z0[3];
 769
 770   d0z0[0] = -999.;
 771   d0z0[1] = -999.;
 772   covd0z0[0] = -999.;
 773   covd0z0[1] = -999.;
 774   covd0z0[2] = -999.;
 775
 776   d0[0] = d0z0[0];
 777   d0err[0] = TMath::Sqrt(TMath::Abs(covd0z0[0]));
 778
 779   d0[1] = d0z0[0];
 780   d0err[1] = TMath::Sqrt(TMath::Abs(covd0z0[0]));
 781
 782   // create the object AliAODRecoDecayLF2Prong
 783   //  AliAODRecoDecayLF2Prong *the2Prong = new AliAODRecoDecayLF2Prong(secVert,px,py,pz,d0,d0err,dcap1n1);
 784   AliAODRecoDecayLF2Prong *the2Prong = new AliAODRecoDecayLF2Prong(secVert,px,py,pz,d0,d0err,dca);
 785   the2Prong->SetOwnPrimaryVtx(primVertexAOD);
 786
 787   //  the2Prong->SetProngIDs(2,{-999,999});
 788   // UShort_t id0[2]={99999,999999};
 789   // the2Prong->SetProngIDs(2,id0);
 790
 791   UShort_t id[2]={(UShort_t)postrack->GetID(),(UShort_t)negtrack->GetID()};
 792   the2Prong->SetProngIDs(2,id);
 793
 794   // cout<<"\n\n\nMake 2 Prong: id[0]"<<id[0]<<" id[1]: "<<id[1]<<endl;
 795   // cout<<"Get: 1 "<<the2Prong->GetProngID(0)<<" 2 "<<the2Prong->GetProngID(1)<<endl;
 796   // cout<<"Get: 3 "<<the2Prong->GetProngID(2)<<" 4 "<<the2Prong->GetProngID(3)<<endl<<endl<<endl;
 797   //delete primVertexAOD; primVertexAOD=NULL;
 798
 799   if(postrack->Charge()!=0 && negtrack->Charge()!=0) {
 800
 801     AddDaughterRefs(secVert,(AliAODEvent&)evento,twoTrackArray);
 802     //    AddDaughterRefs(secVert,(AliAODEvent*)evento,twoTrackArray);
 803
 804   }
 805
 806   return the2Prong;
 807
 808   delete the2Prong;
 809 }
 810
 811 //----------------------------------------------------------------------------
 812 void AliAODMCNuclExReplicator::AddDaughterRefs(AliAODVertex *v,
 813                                             const AliAODEvent &event,
 814                                             const TObjArray *trkArray) const
 815
 816 {
 817   // Add the AOD tracks as daughters of the vertex (TRef)
 818   // AliCodeTimerAuto("",0);
 819   // cout<<"Inside  AddDaughterRefs "<<endl;
 820
 821   Int_t nDg = v->GetNDaughters();
 822
 823   TObject *dg = 0;
 824   if(nDg) dg = v->GetDaughter(0);
 825   if(dg) return; // daughters already added
 826
 827   Int_t nTrks = trkArray->GetEntriesFast();
 828
 829   AliExternalTrackParam *track = 0;
 830   AliAODTrack *aodTrack = 0;
 831   Int_t id;
 832
 833   for(Int_t i=0; i<nTrks; i++) {
 834     track = (AliExternalTrackParam*)trkArray->UncheckedAt(i);
 835
 836     id = (Int_t)track->GetID();
 837     //    printf("---> %d\n",id);
 838     if(id<0) continue; // this track is a AliAODRecoDecay
 839
 840     aodTrack = (AliAODTrack*)event.GetTrack(fAODMap[id]);
 841     v->AddDaughter(aodTrack);
 842   }
 843   return;
 844 }
 845 //----------------------------------------------------------------------------
 846
 847 void AliAODMCNuclExReplicator::AddRefs(AliAODVertex *v,AliAODRecoDecayLF *rd,
 848                                     const AliAODEvent &event,
 849                                     const TObjArray *trkArray) const
 850
 851 {
 852   // Add the AOD tracks as daughters of the vertex (TRef)
 853   // Also add the references to the primary vertex and to the cuts
 854   //AliCodeTimerAuto("",0);
 855
 856   AddDaughterRefs(v,event,trkArray);
 857   rd->SetPrimaryVtxRef((AliAODVertex*)event.GetPrimaryVertex());
 858   return;
 859 }
 860
 861 //---------------------------------------------------------------------------
 862
 863 void AliAODMCNuclExReplicator::Terminate(){
 864
 865 }