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
 568     if(mom<0.150)continue;
 569
 570     Double_t nSigmaNegPion = TMath::Abs(fPIDResponse->NumberOfSigmasTPC(aodtrack,(AliPID::EParticleType)fpartType1));
 571     //    cout<<"%%% nsigma Pi: "<<nSigmaNegPion<<endl;
 572     //    Double_t nSigmaNegPion = TMath::Abs((aodtrack->GetTPCsignal() - foPion->Eval(mom))/foPion->Eval(mom))/0.07;
 573
 574     if ( nSigmaNegPion <= fnSigmaTrk1 /*&& aodtrack->Charge()==-1*/){
 575       Track0[nTrack0++]=j;
 576     }
 577
 578     Double_t nSigmaNuclei = fPIDResponse->NumberOfSigmasTPC(aodtrack,(AliPID::EParticleType)fpartType2);  //To be changed
 579     //cout<<"%%% nsigma Nuclei: "<<nSigmaNuclei<<endl;
 580
 581     if(nSigmaNuclei>-fnSigmaTrk2 && aodtrack->GetTPCsignal()<1000 && mom>0.2){
 582
 583       //Double_t triggerDeDx = 4*AliExternalTrackParam::BetheBlochAleph((mom*2)/(0.938*3),1.0288,31.9806,5.04114e-11,2.13096,2.38541);
 584       //if(aodtrack->GetTPCsignal() > triggerDeDx && aodtrack->GetTPCsignal()<5000 && mom>0.2 /*&& aodtrack->Charge()==1*/){
 585
 586       Track1[nTrack1++]=j;
 587       new((*fNuclei)[nnuclei++]) AliAODTrack(*aodtrack);
 588       //cout<<"@@@ n sigma NUCLEI"<<nSigmaNuclei<<endl;
 589     }
 590
 591   }
 592
 593   //  cout<<"Delete AodTrack..."<<endl;
 594   //  delete aodtrack;
 595   //  aodtrack->Delete();
 596
 597   //Set Track Daughters
 598   Track0.Set(nTrack0);
 599   Track1.Set(nTrack1);
 600
 601   //  cout<<"Track loops..."<<endl;
 602
 603   AliAODTrack *postrackAOD = 0;
 604   AliAODTrack *negtrackAOD = 0;
 605
 606   AliESDtrack *postrack = 0;
 607   AliESDtrack *negtrack = 0;
 608
 609   // cout <<"nTrack1 "<< nTrack1<<endl;
 610   // cout <<"nTrack0 "<< nTrack0<<endl;
 611
 612   for (Int_t i=0; i<nTrack1; i++){                            //! He Tracks Loop
 613
 614     //    cout<<"Inside track loop"<<endl;
 615
 616     Int_t Track1idx=Track1[i];
 617
 618     AliVTrack *trackpos = (AliVTrack*)source.GetTrack(Track1idx);
 619
 620     postrackAOD = (AliAODTrack*)trackpos;
 621     postrack = new AliESDtrack(trackpos);
 622
 623     for (Int_t k=0; k <nTrack0 ; k++) {                           //! Pion Tracks Loop
 624
 625       Int_t Track0idx=Track0[k];
 626
 627       AliVTrack *trackneg = (AliVTrack*)source.GetTrack(Track0idx);
 628       negtrackAOD =(AliAODTrack*)trackneg;
 629       negtrack = new AliESDtrack(trackneg);
 630
 631       twoTrackArray->AddAt(negtrack,0);
 632       twoTrackArray->AddAt(postrack,1);
 633
 634       Double_t dcap1n1 = postrack->GetDCA(negtrack,fBzkG,xdummy,ydummy);
 635
 636       Double_t dcap1toPV = TMath::Abs(postrack->GetD(xPrimaryVertex, yPrimaryVertex,fBzkG));
 637       Double_t dcan1toPV = TMath::Abs(negtrack->GetD(xPrimaryVertex, yPrimaryVertex,fBzkG));
 638
 639       if(dcap1n1>fDCAtracksMin)continue;
 640       if((xdummy+ydummy)>fRmax )continue;
 641       if((xdummy+ydummy)< fRmin)continue;
 642
 643       if ( dcan1toPV < fDNmin)
 644         if ( dcap1toPV < fDPmin) continue;
 645
 646       //      cout<<"dcap1n1: "<<dcap1n1<<endl;
 647
 648       AliAODVertex *vertexp1n1 = ReconstructSecondaryVertex(twoTrackArray,dispersion,kTRUE);
 649
 650       if(!vertexp1n1) {
 651
 652         twoTrackArray->Clear();
 653         delete negtrack;
 654         negtrack=NULL;
 655         continue;
 656       }
 657
 658       //      cout<<"vertexp1n1:  "<<vertexp1n1<<endl;
 659
 660       io2Prong = Make2Prong(twoTrackArray,source,vertexp1n1,dcap1n1);
 661
 662       //cout<<"Fuori dal loop e' quello che salvo \npr0: "<<io2Prong->GetProngID(0)<<" pr1: "<<io2Prong->GetProngID(1)<<" pr2"<<io2Prong->GetProngID(2)<<endl;
 663
 664       if(io2Prong->CosPointingAngle()<fCosMin)continue;
 665
 666       //      cout<<"if CPA \npr0: "<<io2Prong->GetProngID(0)<<" pr1: "<<io2Prong->GetProngID(1)<<" pr2"<<io2Prong->GetProngID(2)<<endl;
 667       // cout<<"pointing angle "<<io2Prong->CosPointingAngle()<<endl;
 668
 669
 670       AliAODTrack *trk0 = (AliAODTrack*)io2Prong->GetDaughter(0);
 671       AliAODTrack *trk1 = (AliAODTrack*)io2Prong->GetDaughter(1);
 672
 673       cout<<"**********************************************"<<endl;
 674       cout<<trk0/*->GetID()*/<<" "<<negtrackAOD->GetID()<<endl;
 675       cout<<trk1/*->GetID()*/<<" "<<postrackAOD->GetID()<<endl;
 676       cout<<"d0 io2Prong: "<<io2Prong->GetProngID(1)<<endl;
 677       cout<<"d1 io2Prong: "<<io2Prong->GetProngID(0)<<endl;
 678       cout<<"**********************************************"<<endl;
 679
 680       //      rd =  new((*fSecondaryVerices)[nsv++]) AliAODRecoDecayLF2Prong(*io2Prong);
 681
 682       new((*fSecondaryVerices)[nsv++]) AliAODRecoDecayLF2Prong(*io2Prong);
 683
 684       //      cout<<"QUELLO CHE SALVo \npr0: "<<rd->GetProngID(0)<<" pr1: "<<rd->GetProngID(1)<<" pr2"<<rd->GetProngID(2)<<endl;
 685
 686       // rd->SetSecondaryVtx(vertexp1n1);
 687       // vertexp1n1->SetParent(rd);
 688       // AddRefs(vertexp1n1,rd,source,twoTrackArray);
 689
 690
 691       new((*fDaughterTracks)[ntracksD++]) AliAODTrack(*negtrackAOD);
 692       new((*fDaughterTracks)[ntracksD++]) AliAODTrack(*postrackAOD);
 693
 694       delete negtrack;
 695       negtrack = NULL;
 696
 697       delete vertexp1n1;
 698       vertexp1n1= NULL;
 699       continue;
 700     }
 701
 702     delete postrack;
 703     postrack = NULL;
 704
 705   }
 706
 707   //----------------------------------------------------------
 708
 709   assert(fVertices!=0x0);
 710   fVertices->Clear("C");
 711   TIter nextV(source.GetVertices());
 712   AliAODVertex* v;
 713   Int_t nvertices(0);
 714
 715   while ( ( v = static_cast<AliAODVertex*>(nextV()) ) )
 716     {
 717       if (  v->GetType() == AliAODVertex::kPrimary     ||
 718             v->GetType() == AliAODVertex::kMainSPD     ||
 719             v->GetType() == AliAODVertex::kPileupSPD   ||
 720             v->GetType() == AliAODVertex::kPileupTracks||
 721             v->GetType() == AliAODVertex::kMainTPC  )
 722         {
 723           AliAODVertex* tmp = v->CloneWithoutRefs();
 724           AliAODVertex* copiedVertex = new((*fVertices)[nvertices++]) AliAODVertex(*tmp);
 725
 726           // to insure the main vertex retains the ncontributors information
 727           // (which is otherwise computed dynamically from
 728           // references to tracks, which we do not keep in muon aods...)
 729           // we set it here
 730
 731           copiedVertex->SetNContributors(v->GetNContributors());
 732
 733           //  fVertices->Delete();
 734           //      delete copiedVertex;
 735           delete tmp;
 736         }
 737     }
 738
 739   //  printf("....Done NuclEx Replicator...\n");
 740
 741   AliDebug(1,Form("input mu tracks=%d tracks=%d vertices=%d nnuclei=%d",
 742                   input,fSecondaryVerices->GetEntries(),fVertices->GetEntries(),fNuclei->GetEntries()));
 743
 744   // Finally, deal with MC information, if needed
 745
 746   // if ( fMCMode > 0 )
 747   //   {
 748   //     FilterMC(source);
 749   //   }
 750
 751   //cout<<"Delete..."<<endl;
 752   // delete foPion;
 753   // foPion = NULL;
 754   //cout<<"Delete 1"<<  endl;
 755
 756   if(io2Prong) {delete io2Prong; io2Prong=NULL;}
 757   //cout<<"Delete 2"<<  endl;
 758   twoTrackArray->Delete();  delete twoTrackArray;
 759   //cout<<"Delete 3"<<  endl;
 760   // vtx->Delete();  delete vtx;
 761   //cout<<"Delete 4"<<  endl;
 762   if(fV1) { delete fV1; fV1=NULL; }
 763   //cout<<"Delete 5"<<  endl;
 764   if(fAODMap) { delete [] fAODMap; fAODMap=NULL; }
 765   delete []indices;
 766   //cout<<"Delete 6"<<  endl;
 767   delete fVertexerTracks;
 768
 769
 770 }
 771
 772 //-----------------------------------------------------------------------------
 773
 774 AliAODVertex *AliAODNuclExReplicator::ReconstructSecondaryVertex(TObjArray *trkArray,
 775                                                                 Double_t &dispersion,Bool_t useTRefArray) const
 776
 777 {
 778   // Secondary vertex reconstruction with AliVertexerTracks or AliKFParticle
 779   //AliCodeTimerAuto("",0);
 780
 781   AliESDVertex *vertexESD = 0;
 782   AliAODVertex *vertexAOD = 0;
 783
 784   if(!fSecVtxWithKF) { // AliVertexerTracks
 785
 786     fVertexerTracks->SetVtxStart(fV1);
 787     vertexESD = (AliESDVertex*)fVertexerTracks->VertexForSelectedESDTracks(trkArray);
 788
 789     if(!vertexESD) return vertexAOD;
 790
 791
 792     Double_t vertRadius2=vertexESD->GetXv()*vertexESD->GetXv()+vertexESD->GetYv()*vertexESD->GetYv();
 793     if(vertRadius2>200.){
 794       // vertex outside beam pipe, reject candidate to avoid propagation through material
 795       delete vertexESD; vertexESD=NULL;
 796       return vertexAOD;
 797     }
 798
 799   } else { // Kalman Filter vertexer (AliKFParticle)
 800
 801     AliKFParticle::SetField(fBzkG);
 802
 803     AliKFVertex vertexKF;
 804
 805     Int_t nTrks = trkArray->GetEntriesFast();
 806     for(Int_t i=0; i<nTrks; i++) {
 807       AliESDtrack *esdTrack = (AliESDtrack*)trkArray->At(i);
 808       AliKFParticle daughterKF(*esdTrack,211);
 809       vertexKF.AddDaughter(daughterKF);
 810     }
 811     vertexESD = new AliESDVertex(vertexKF.Parameters(),
 812                                  vertexKF.CovarianceMatrix(),
 813                                  vertexKF.GetChi2(),
 814                                  vertexKF.GetNContributors());
 815
 816   }
 817   // convert to AliAODVertex
 818   Double_t pos[3],cov[6],chi2perNDF;
 819   vertexESD->GetXYZ(pos); // position
 820   vertexESD->GetCovMatrix(cov); //covariance matrix
 821   chi2perNDF = vertexESD->GetChi2toNDF();
 822   dispersion = vertexESD->GetDispersion();
 823   delete vertexESD;
 824   vertexESD=NULL;
 825
 826   Int_t nprongs= (useTRefArray ? 0 : trkArray->GetEntriesFast());
 827   vertexAOD = new AliAODVertex(pos,cov,chi2perNDF,0x0,-1,AliAODVertex::kUndef,nprongs);
 828
 829   //cout<<"------------------> Reconstruct vertexAOD: "<<vertexAOD<<endl;
 830
 831   return vertexAOD;
 832
 833
 834 }
 835
 836 //-----------------------------------------------------------------------------
 837
 838 AliAODRecoDecayLF2Prong* AliAODNuclExReplicator::Make2Prong(TObjArray *twoTrackArray,const AliAODEvent &evento,
 839                                                            AliAODVertex *secVert,Double_t dca)
 840
 841
 842 {
 843
 844   //cout<<"Inside make 2 prong"<<endl;
 845
 846   Int_t charge[2];
 847   charge[0]=1; //it was -1 //Ramona
 848   charge[1]=2;
 849
 850   // From AliAnalysisVertexingLF.cxx Method:Make2Prongs
 851
 852   fBzkG = evento.GetMagneticField();
 853
 854   Double_t px[2],py[2],pz[2],d0[2],d0err[2];
 855   // Also this has been changed //Ramona
 856   AliESDtrack *negtrack = (AliESDtrack*)twoTrackArray->UncheckedAt(0);
 857   AliESDtrack *postrack = (AliESDtrack*)twoTrackArray->UncheckedAt(1);
 858
 859   //cout<<"negtrack: "<<negtrack<<" postrack: "<<postrack<<endl;
 860   //cout<<"kVeryBig: "<<kVeryBig<<endl;
 861   //cout<<"secVert: "<<secVert<<endl;
 862
 863   // // propagate tracks to secondary vertex, to compute inv. mass
 864
 865   postrack->PropagateToDCA(secVert,fBzkG,kVeryBig);
 866   negtrack->PropagateToDCA(secVert,fBzkG,kVeryBig);
 867
 868   Double_t momentum[3];
 869
 870   negtrack->GetPxPyPz(momentum);
 871   px[0] = charge[0]*momentum[0];
 872   py[0] = charge[0]*momentum[1];
 873   pz[0] = charge[0]*momentum[2];
 874
 875   postrack->GetPxPyPz(momentum);
 876   px[1] = charge[1]*momentum[0];
 877   py[1] = charge[1]*momentum[1];
 878   pz[1] = charge[1]*momentum[2];
 879
 880   //cout<< px[0] <<" "<< " "<< py[0] << " "<< pz[0]<<endl;
 881   //  px[0] = momentum[0]; py[0] = momentum[1]; pz[0] = momentum[2];
 882   //cout<< px[1] <<" "<< " "<< py[1] << " "<< pz[1]<<endl;
 883   //px[1] = momentum[0]; py[1] = momentum[1]; pz[1] = momentum[2];
 884
 885
 886   // primary vertex to be used by this candidate
 887   AliAODVertex *primVertexAOD  = evento.GetPrimaryVertex();
 888   //cout<<"primVertexAOD "<<primVertexAOD<<endl;
 889   if(!primVertexAOD) return 0x0;
 890
 891   Double_t d0z0[2],covd0z0[3];
 892
 893   d0z0[0] = -999.;
 894   d0z0[1] = -999.;
 895   covd0z0[0] = -999.;
 896   covd0z0[1] = -999.;
 897   covd0z0[2] = -999.;
 898
 899   d0[0] = d0z0[0];
 900   d0err[0] = TMath::Sqrt(TMath::Abs(covd0z0[0]));
 901
 902   d0[1] = d0z0[0];
 903   d0err[1] = TMath::Sqrt(TMath::Abs(covd0z0[0]));
 904
 905   // create the object AliAODRecoDecayLF2Prong
 906   //  AliAODRecoDecayLF2Prong *the2Prong = new AliAODRecoDecayLF2Prong(secVert,px,py,pz,d0,d0err,dcap1n1);
 907   AliAODRecoDecayLF2Prong *the2Prong = new AliAODRecoDecayLF2Prong(secVert,px,py,pz,d0,d0err,dca);
 908   the2Prong->SetOwnPrimaryVtx(primVertexAOD);
 909
 910   //  the2Prong->SetProngIDs(2,{-999,999});
 911   // UShort_t id0[2]={99999,999999};
 912   // the2Prong->SetProngIDs(2,id0);
 913
 914   UShort_t id[2]={(UShort_t)postrack->GetID(),(UShort_t)negtrack->GetID()};
 915   the2Prong->SetProngIDs(2,id);
 916
 917   // cout<<"\n\n\nMake 2 Prong: id[0]"<<id[0]<<" id[1]: "<<id[1]<<endl;
 918   // cout<<"Get: 1 "<<the2Prong->GetProngID(0)<<" 2 "<<the2Prong->GetProngID(1)<<endl;
 919   // cout<<"Get: 3 "<<the2Prong->GetProngID(2)<<" 4 "<<the2Prong->GetProngID(3)<<endl<<endl<<endl;
 920   //delete primVertexAOD; primVertexAOD=NULL;
 921
 922   if(postrack->Charge()!=0 && negtrack->Charge()!=0) {
 923
 924     AddDaughterRefs(secVert,(AliAODEvent&)evento,twoTrackArray);
 925     //    AddDaughterRefs(secVert,(AliAODEvent*)evento,twoTrackArray);
 926
 927   }
 928
 929   return the2Prong;
 930
 931   delete the2Prong;
 932 }
 933
 934 //----------------------------------------------------------------------------
 935 void AliAODNuclExReplicator::AddDaughterRefs(AliAODVertex *v,
 936                                             const AliAODEvent &event,
 937                                             const TObjArray *trkArray) const
 938
 939 {
 940   // Add the AOD tracks as daughters of the vertex (TRef)
 941   // AliCodeTimerAuto("",0);
 942   // cout<<"Inside  AddDaughterRefs "<<endl;
 943
 944   Int_t nDg = v->GetNDaughters();
 945
 946   TObject *dg = 0;
 947   if(nDg) dg = v->GetDaughter(0);
 948   if(dg) return; // daughters already added
 949
 950   Int_t nTrks = trkArray->GetEntriesFast();
 951
 952   AliExternalTrackParam *track = 0;
 953   AliAODTrack *aodTrack = 0;
 954   Int_t id;
 955
 956   for(Int_t i=0; i<nTrks; i++) {
 957     track = (AliExternalTrackParam*)trkArray->UncheckedAt(i);
 958
 959     id = (Int_t)track->GetID();
 960     //    printf("---> %d\n",id);
 961     if(id<0) continue; // this track is a AliAODRecoDecay
 962
 963     aodTrack = (AliAODTrack*)event.GetTrack(fAODMap[id]);
 964     v->AddDaughter(aodTrack);
 965   }
 966   return;
 967 }
 968 //----------------------------------------------------------------------------
 969
 970 void AliAODNuclExReplicator::AddRefs(AliAODVertex *v,AliAODRecoDecayLF *rd,
 971                                     const AliAODEvent &event,
 972                                     const TObjArray *trkArray) const
 973
 974 {
 975   // Add the AOD tracks as daughters of the vertex (TRef)
 976   // Also add the references to the primary vertex and to the cuts
 977   //AliCodeTimerAuto("",0);
 978
 979   AddDaughterRefs(v,event,trkArray);
 980   rd->SetPrimaryVtxRef((AliAODVertex*)event.GetPrimaryVertex());
 981   return;
 982 }
 983
 984 //---------------------------------------------------------------------------
 985
 986 void AliAODNuclExReplicator::Terminate(){
 987
 988 }