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 }