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       AliAODHeader * header = dynamic_cast<AliAODHeader*>(source.GetHeader());
 442       if(!header) AliFatal("Not a standard AOD");
 443       *fHeader = *(header);
 444     }
 445
 446   fVertices->Clear("C");
 447
 448   fNuclei->Clear("C");
 449
 450   fSecondaryVerices->Clear("C");
 451
 452   fDaughterTracks->Clear("C");
 453
 454   //----------------------------------
 455
 456   //  cout<<"Centrality AOD source: "<<source.GetHeader()->GetCentrality()<<endl;
 457
 458   Int_t nsv(0);
 459   Int_t nnuclei(0);
 460   Int_t ntracksD(0);
 461
 462   Int_t input(0);
 463   Double_t xdummy,ydummy;
 464
 465   AliAODRecoDecayLF2Prong *io2Prong  = 0;
 466
 467   TObjArray *twoTrackArray    = new TObjArray(2);
 468   Double_t dispersion;
 469
 470   // cout<<"Qui"<<endl;
 471   //cout<<source.GetMagneticField()<<endl;
 472
 473   AliAODVertex *vtx = source.GetPrimaryVertex();
 474
 475   //  cout<<"Source "<<source<<endl;
 476   //cout<<"vtx: "<<vtx<<endl;
 477
 478   // A Set of very loose cut for a weak strange decay
 479
 480   fCosMin       = 0.99;
 481   fDCAtracksMin = 1;
 482   fRmax         = 200.;
 483   fRmin         = 0.1;
 484   fDNmin        = 0.05;
 485   fDPmin        = 0.05;
 486
 487   //----------------------------------------------------------
 488
 489   //  Int_t nindices=0;
 490   UShort_t *indices = 0;
 491   const Int_t entries = source.GetNumberOfTracks();
 492
 493   Double_t pos[3],cov[6];
 494   vtx->GetXYZ(pos);
 495   vtx->GetCovarianceMatrix(cov);
 496   fV1 = new AliESDVertex(pos,cov,100.,100,vtx->GetName());
 497   //  cout<<"fV1 pricipal loop: "<<fV1<<endl;
 498
 499   if(entries<=0) return;
 500   indices = new UShort_t[entries];
 501   memset(indices,0,sizeof(UShort_t)*entries);
 502   fAODMapSize = 100000;
 503   fAODMap = new Int_t[fAODMapSize];
 504   memset(fAODMap,0,sizeof(Int_t)*fAODMapSize);
 505   //  cent=((AliAODEvent&)source)->GetCentrality();
 506
 507   //-------------------------------------------------------------
 508
 509   //AliAODRecoDecayLF   *rd = 0;
 510   //  AliAODRecoDecayLF2Prong*rd = 0;
 511
 512
 513   if(vtx->GetNContributors()<1) {
 514
 515     // SPD vertex cut
 516     vtx =source.GetPrimaryVertexSPD();
 517
 518     if(vtx->GetNContributors()<1) {
 519       Info("AliAnalysisTaskHelium3Pi","No good vertex, skip event");
 520       return; // NO GOOD VERTEX, SKIP EVENT
 521     }
 522   }
 523
 524   Double_t xPrimaryVertex=0.,yPrimaryVertex=0.;
 525   xPrimaryVertex=vtx->GetX();
 526   yPrimaryVertex=vtx->GetY();
 527
 528   fBzkG=source.GetMagneticField();
 529   fVertexerTracks=new AliVertexerTracks(fBzkG);
 530
 531   // 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);
 532   // foPion->SetParameters(4.1,8.98482806165147636e+00,1.54000000000000005e-05,2.30445734159456084e+00,2.25624744086878559e+00);
 533
 534   Double_t TrackNumber = source.GetNumberOfTracks();
 535
 536   //Tracks arrays
 537
 538   TArrayI Track0(TrackNumber);        //Pions
 539   Int_t nTrack0=0;
 540
 541   TArrayI Track1(TrackNumber);        //Helium3
 542   Int_t nTrack1=0;
 543
 544   for(Int_t j=0; j<TrackNumber; j++){
 545
 546     //    cout<<"Inside loop tracks"<<endl;
 547
 548     AliVTrack *track = (AliVTrack*)source.GetTrack(j);
 549
 550     AliAODTrack *aodtrack =(AliAODTrack*)track;
 551
 552     //-----------------------------------------------------------
 553     //Track cuts
 554     if(aodtrack->GetTPCNcls() < 70 )continue;
 555     if(aodtrack->Chi2perNDF() > 4 )continue;
 556
 557     // cout<<"Filter Map: "<<aodtrack->GetFilterMap()<<endl;
 558     //    if (!aodtrack->TestFilterMask(BIT(4))) continue;
 559     // if(aodtrack->TestFilterMask(BIT(4))!=0)
 560     //   cout<<"!!!!!!!!!!!!!!!!!!!!!!!!!!! Filter bit mask: !!!!!!!!!!!!!!!!!!!!!!!!!!!"<<endl;
 561
 562     if (!aodtrack->IsOn(AliAODTrack::kTPCrefit)) continue;
 563     if (!aodtrack->IsOn(AliAODTrack::kTPCin)) continue;
 564     if (aodtrack->IsOn(AliAODTrack::kITSpureSA)) continue;
 565
 566     //---------------------------------------------------------------
 567
 568     //    Double_t mom = aodtrack->GetDetPid()->GetTPCmomentum();
 569     Double_t mom = aodtrack->P();
 570
 571     if(mom<0.150)continue;
 572
 573
 574
 575     Double_t nSigmaNegPion = TMath::Abs(fPIDResponse->NumberOfSigmasTPC(aodtrack,(AliPID::EParticleType)fpartType1));
 576     //    cout<<"%%% nsigma Pi: "<<nSigmaNegPion<<endl;
 577     //    Double_t nSigmaNegPion = TMath::Abs((aodtrack->GetTPCsignal() - foPion->Eval(mom))/foPion->Eval(mom))/0.07;
 578
 579     if ( nSigmaNegPion <= fnSigmaTrk1 /*&& aodtrack->Charge()==-1*/){
 580       Track0[nTrack0++]=j;
 581     }
 582
 583     Double_t nSigmaNuclei = fPIDResponse->NumberOfSigmasTPC(aodtrack,(AliPID::EParticleType)fpartType2);  //To be changed
 584     //cout<<"%%% nsigma Nuclei: "<<nSigmaNuclei<<endl;
 585
 586     if(nSigmaNuclei>-fnSigmaTrk2 && aodtrack->GetTPCsignal()<1000 && mom>0.2){
 587
 588     //Double_t triggerDeDx = 4*AliExternalTrackParam::BetheBlochAleph((mom*2)/(0.938*3),1.0288,31.9806,5.04114e-11,2.13096,2.38541);
 589       //if(aodtrack->GetTPCsignal() > triggerDeDx && aodtrack->GetTPCsignal()<5000 && mom>0.2 /*&& aodtrack->Charge()==1*/){
 590
 591       Track1[nTrack1++]=j;
 592       new((*fNuclei)[nnuclei++]) AliAODTrack(*aodtrack);
 593       //cout<<"@@@ n sigma NUCLEI"<<nSigmaNuclei<<endl;
 594     }
 595
 596   }
 597
 598   //  cout<<"Delete AodTrack..."<<endl;
 599   //  delete aodtrack;
 600   //  aodtrack->Delete();
 601
 602   //Set Track Daughters
 603   Track0.Set(nTrack0);
 604   Track1.Set(nTrack1);
 605
 606   //  cout<<"Track loops..."<<endl;
 607
 608   AliAODTrack *postrackAOD = 0;
 609   AliAODTrack *negtrackAOD = 0;
 610
 611   AliESDtrack *postrack = 0;
 612   AliESDtrack *negtrack = 0;
 613
 614   // cout <<"nTrack1 "<< nTrack1<<endl;
 615   // cout <<"nTrack0 "<< nTrack0<<endl;
 616
 617   for (Int_t i=0; i<nTrack1; i++){                            //! He Tracks Loop
 618
 619     //    cout<<"Inside track loop"<<endl;
 620
 621     Int_t Track1idx=Track1[i];
 622
 623     AliVTrack *trackpos = (AliVTrack*)source.GetTrack(Track1idx);
 624
 625     postrackAOD = (AliAODTrack*)trackpos;
 626     postrack = new AliESDtrack(trackpos);
 627
 628     for (Int_t k=0; k <nTrack0 ; k++) {                           //! Pion Tracks Loop
 629
 630       Int_t Track0idx=Track0[k];
 631
 632       AliVTrack *trackneg = (AliVTrack*)source.GetTrack(Track0idx);
 633       negtrackAOD =(AliAODTrack*)trackneg;
 634       negtrack = new AliESDtrack(trackneg);
 635
 636       twoTrackArray->AddAt(negtrack,0);
 637       twoTrackArray->AddAt(postrack,1);
 638
 639       Double_t dcap1n1 = postrack->GetDCA(negtrack,fBzkG,xdummy,ydummy);
 640
 641       Double_t dcap1toPV = TMath::Abs(postrack->GetD(xPrimaryVertex, yPrimaryVertex,fBzkG));
 642       Double_t dcan1toPV = TMath::Abs(negtrack->GetD(xPrimaryVertex, yPrimaryVertex,fBzkG));
 643
 644       if(dcap1n1>fDCAtracksMin)continue;
 645       if((xdummy+ydummy)>fRmax )continue;
 646       if((xdummy+ydummy)< fRmin)continue;
 647
 648       if ( dcan1toPV < fDNmin)
 649         if ( dcap1toPV < fDPmin) continue;
 650
 651       //      cout<<"dcap1n1: "<<dcap1n1<<endl;
 652
 653       AliAODVertex *vertexp1n1 = ReconstructSecondaryVertex(twoTrackArray,dispersion,kTRUE);
 654
 655       if(!vertexp1n1) {
 656
 657         twoTrackArray->Clear();
 658         delete negtrack;
 659         negtrack=NULL;
 660         continue;
 661       }
 662
 663       //      cout<<"vertexp1n1:  "<<vertexp1n1<<endl;
 664
 665       io2Prong = Make2Prong(twoTrackArray,source,vertexp1n1,dcap1n1);
 666
 667       //cout<<"Fuori dal loop e' quello che salvo \npr0: "<<io2Prong->GetProngID(0)<<" pr1: "<<io2Prong->GetProngID(1)<<" pr2"<<io2Prong->GetProngID(2)<<endl;
 668
 669       if(io2Prong->CosPointingAngle()<fCosMin)continue;
 670
 671       //      cout<<"if CPA \npr0: "<<io2Prong->GetProngID(0)<<" pr1: "<<io2Prong->GetProngID(1)<<" pr2"<<io2Prong->GetProngID(2)<<endl;
 672       // cout<<"pointing angle "<<io2Prong->CosPointingAngle()<<endl;
 673
 674
 675       // AliAODTrack *trk0 = (AliAODTrack*)io2Prong->GetDaughter(0);
 676       // AliAODTrack *trk1 = (AliAODTrack*)io2Prong->GetDaughter(1);
 677
 678       // cout<<"**********************************************"<<endl;
 679       // cout<<trk0/*->GetID()*/<<" "<<negtrackAOD->GetID()<<endl;
 680       // cout<<trk1/*->GetID()*/<<" "<<postrackAOD->GetID()<<endl;
 681       // cout<<"d0 io2Prong: "<<io2Prong->GetProngID(1)<<endl;
 682       // cout<<"d1 io2Prong: "<<io2Prong->GetProngID(0)<<endl;
 683       // cout<<"**********************************************"<<endl;
 684
 685       //      rd =  new((*fSecondaryVerices)[nsv++]) AliAODRecoDecayLF2Prong(*io2Prong);
 686
 687       new((*fSecondaryVerices)[nsv++]) AliAODRecoDecayLF2Prong(*io2Prong);
 688
 689       //      cout<<"QUELLO CHE SALVo \npr0: "<<rd->GetProngID(0)<<" pr1: "<<rd->GetProngID(1)<<" pr2"<<rd->GetProngID(2)<<endl;
 690
 691       // rd->SetSecondaryVtx(vertexp1n1);
 692       // vertexp1n1->SetParent(rd);
 693       // AddRefs(vertexp1n1,rd,source,twoTrackArray);
 694
 695
 696       new((*fDaughterTracks)[ntracksD++]) AliAODTrack(*negtrackAOD);
 697       new((*fDaughterTracks)[ntracksD++]) AliAODTrack(*postrackAOD);
 698
 699       delete negtrack;
 700       negtrack = NULL;
 701
 702       delete vertexp1n1;
 703       vertexp1n1= NULL;
 704       continue;
 705     }
 706
 707     delete postrack;
 708     postrack = NULL;
 709
 710   }
 711
 712   //----------------------------------------------------------
 713
 714   assert(fVertices!=0x0);
 715   fVertices->Clear("C");
 716   TIter nextV(source.GetVertices());
 717   AliAODVertex* v;
 718   Int_t nvertices(0);
 719
 720   while ( ( v = static_cast<AliAODVertex*>(nextV()) ) )
 721     {
 722        if (  v->GetType() == AliAODVertex::kPrimary     ||
 723             v->GetType() == AliAODVertex::kMainSPD     ||
 724             v->GetType() == AliAODVertex::kPileupSPD   ||
 725             v->GetType() == AliAODVertex::kPileupTracks||
 726             v->GetType() == AliAODVertex::kMainTPC  )
 727         {
 728           AliAODVertex* tmp = v->CloneWithoutRefs();
 729           AliAODVertex* copiedVertex = new((*fVertices)[nvertices++]) AliAODVertex(*tmp);
 730
 731           // to insure the main vertex retains the ncontributors information
 732           // (which is otherwise computed dynamically from
 733           // references to tracks, which we do not keep in muon aods...)
 734           // we set it here
 735
 736           copiedVertex->SetNContributors(v->GetNContributors());
 737
 738           //  fVertices->Delete();
 739           //      delete copiedVertex;
 740           delete tmp;
 741         }
 742     }
 743
 744   //  printf("....Done NuclEx Replicator...\n");
 745
 746   AliDebug(1,Form("input mu tracks=%d tracks=%d vertices=%d nnuclei=%d",
 747                   input,fSecondaryVerices->GetEntries(),fVertices->GetEntries(),fNuclei->GetEntries()));
 748
 749   // Finally, deal with MC information, if needed
 750
 751   // if ( fMCMode > 0 )
 752   //   {
 753   //     FilterMC(source);
 754   //   }
 755
 756   //cout<<"Delete..."<<endl;
 757   // delete foPion;
 758   // foPion = NULL;
 759   //cout<<"Delete 1"<<  endl;
 760
 761   if(io2Prong) {delete io2Prong; io2Prong=NULL;}
 762   //cout<<"Delete 2"<<  endl;
 763   twoTrackArray->Delete();  delete twoTrackArray;
 764   //cout<<"Delete 3"<<  endl;
 765   // vtx->Delete();  delete vtx;
 766   //cout<<"Delete 4"<<  endl;
 767   if(fV1) { delete fV1; fV1=NULL; }
 768   //cout<<"Delete 5"<<  endl;
 769   if(fAODMap) { delete [] fAODMap; fAODMap=NULL; }
 770   delete []indices;
 771   //cout<<"Delete 6"<<  endl;
 772   delete fVertexerTracks;
 773
 774
 775 }
 776
 777 //-----------------------------------------------------------------------------
 778
 779 AliAODVertex *AliAODNuclExReplicator::ReconstructSecondaryVertex(TObjArray *trkArray,
 780                                                                 Double_t &dispersion,Bool_t useTRefArray) const
 781
 782 {
 783   // Secondary vertex reconstruction with AliVertexerTracks or AliKFParticle
 784   //AliCodeTimerAuto("",0);
 785
 786   AliESDVertex *vertexESD = 0;
 787   AliAODVertex *vertexAOD = 0;
 788
 789   if(!fSecVtxWithKF) { // AliVertexerTracks
 790
 791     fVertexerTracks->SetVtxStart(fV1);
 792     vertexESD = (AliESDVertex*)fVertexerTracks->VertexForSelectedESDTracks(trkArray);
 793
 794     if(!vertexESD) return vertexAOD;
 795
 796
 797     Double_t vertRadius2=vertexESD->GetX()*vertexESD->GetX()+vertexESD->GetY()*vertexESD->GetY();
 798     if(vertRadius2>200.){
 799       // vertex outside beam pipe, reject candidate to avoid propagation through material
 800       delete vertexESD; vertexESD=NULL;
 801       return vertexAOD;
 802     }
 803
 804   } else { // Kalman Filter vertexer (AliKFParticle)
 805
 806     AliKFParticle::SetField(fBzkG);
 807
 808     AliKFVertex vertexKF;
 809
 810     Int_t nTrks = trkArray->GetEntriesFast();
 811     for(Int_t i=0; i<nTrks; i++) {
 812       AliESDtrack *esdTrack = (AliESDtrack*)trkArray->At(i);
 813       AliKFParticle daughterKF(*esdTrack,211);
 814       vertexKF.AddDaughter(daughterKF);
 815     }
 816     vertexESD = new AliESDVertex(vertexKF.Parameters(),
 817                                  vertexKF.CovarianceMatrix(),
 818                                  vertexKF.GetChi2(),
 819                                  vertexKF.GetNContributors());
 820
 821   }
 822   // convert to AliAODVertex
 823   Double_t pos[3],cov[6],chi2perNDF;
 824   vertexESD->GetXYZ(pos); // position
 825   vertexESD->GetCovMatrix(cov); //covariance matrix
 826   chi2perNDF = vertexESD->GetChi2toNDF();
 827   dispersion = vertexESD->GetDispersion();
 828   delete vertexESD;
 829   vertexESD=NULL;
 830
 831   Int_t nprongs= (useTRefArray ? 0 : trkArray->GetEntriesFast());
 832   vertexAOD = new AliAODVertex(pos,cov,chi2perNDF,0x0,-1,AliAODVertex::kUndef,nprongs);
 833
 834   //cout<<"------------------> Reconstruct vertexAOD: "<<vertexAOD<<endl;
 835
 836   return vertexAOD;
 837
 838
 839 }
 840
 841 //-----------------------------------------------------------------------------
 842
 843 AliAODRecoDecayLF2Prong* AliAODNuclExReplicator::Make2Prong(TObjArray *twoTrackArray,const AliAODEvent &evento,
 844                                                            AliAODVertex *secVert,Double_t dca)
 845
 846
 847 {
 848
 849   //cout<<"Inside make 2 prong"<<endl;
 850
 851   Int_t charge[2];
 852   charge[0]=1; //it was -1 //Ramona
 853   charge[1]=2;
 854
 855   // From AliAnalysisVertexingLF.cxx Method:Make2Prongs
 856
 857   fBzkG = evento.GetMagneticField();
 858
 859   Double_t px[2],py[2],pz[2],d0[2],d0err[2];
 860   // Also this has been changed //Ramona
 861   AliESDtrack *negtrack = (AliESDtrack*)twoTrackArray->UncheckedAt(0);
 862   AliESDtrack *postrack = (AliESDtrack*)twoTrackArray->UncheckedAt(1);
 863
 864   //cout<<"negtrack: "<<negtrack<<" postrack: "<<postrack<<endl;
 865   //cout<<"kVeryBig: "<<kVeryBig<<endl;
 866   //cout<<"secVert: "<<secVert<<endl;
 867
 868   // // propagate tracks to secondary vertex, to compute inv. mass
 869
 870   postrack->PropagateToDCA(secVert,fBzkG,kVeryBig);
 871   negtrack->PropagateToDCA(secVert,fBzkG,kVeryBig);
 872
 873   Double_t momentum[3];
 874
 875   negtrack->GetPxPyPz(momentum);
 876   px[0] = charge[0]*momentum[0];
 877   py[0] = charge[0]*momentum[1];
 878   pz[0] = charge[0]*momentum[2];
 879
 880   postrack->GetPxPyPz(momentum);
 881   px[1] = charge[1]*momentum[0];
 882   py[1] = charge[1]*momentum[1];
 883   pz[1] = charge[1]*momentum[2];
 884
 885   //cout<< px[0] <<" "<< " "<< py[0] << " "<< pz[0]<<endl;
 886   //  px[0] = momentum[0]; py[0] = momentum[1]; pz[0] = momentum[2];
 887   //cout<< px[1] <<" "<< " "<< py[1] << " "<< pz[1]<<endl;
 888   //px[1] = momentum[0]; py[1] = momentum[1]; pz[1] = momentum[2];
 889
 890
 891   // primary vertex to be used by this candidate
 892   AliAODVertex *primVertexAOD  = evento.GetPrimaryVertex();
 893   //cout<<"primVertexAOD "<<primVertexAOD<<endl;
 894   if(!primVertexAOD) return 0x0;
 895
 896   Double_t d0z0[2],covd0z0[3];
 897
 898   d0z0[0] = -999.;
 899   d0z0[1] = -999.;
 900   covd0z0[0] = -999.;
 901   covd0z0[1] = -999.;
 902   covd0z0[2] = -999.;
 903
 904   d0[0] = d0z0[0];
 905   d0err[0] = TMath::Sqrt(TMath::Abs(covd0z0[0]));
 906
 907   d0[1] = d0z0[0];
 908   d0err[1] = TMath::Sqrt(TMath::Abs(covd0z0[0]));
 909
 910   // create the object AliAODRecoDecayLF2Prong
 911   //  AliAODRecoDecayLF2Prong *the2Prong = new AliAODRecoDecayLF2Prong(secVert,px,py,pz,d0,d0err,dcap1n1);
 912   AliAODRecoDecayLF2Prong *the2Prong = new AliAODRecoDecayLF2Prong(secVert,px,py,pz,d0,d0err,dca);
 913   the2Prong->SetOwnPrimaryVtx(primVertexAOD);
 914
 915   //  the2Prong->SetProngIDs(2,{-999,999});
 916   // UShort_t id0[2]={99999,999999};
 917   // the2Prong->SetProngIDs(2,id0);
 918
 919   UShort_t id[2]={(UShort_t)postrack->GetID(),(UShort_t)negtrack->GetID()};
 920   the2Prong->SetProngIDs(2,id);
 921
 922   // cout<<"\n\n\nMake 2 Prong: id[0]"<<id[0]<<" id[1]: "<<id[1]<<endl;
 923   // cout<<"Get: 1 "<<the2Prong->GetProngID(0)<<" 2 "<<the2Prong->GetProngID(1)<<endl;
 924   // cout<<"Get: 3 "<<the2Prong->GetProngID(2)<<" 4 "<<the2Prong->GetProngID(3)<<endl<<endl<<endl;
 925   //delete primVertexAOD; primVertexAOD=NULL;
 926
 927   if(postrack->Charge()!=0 && negtrack->Charge()!=0) {
 928
 929     AddDaughterRefs(secVert,(AliAODEvent&)evento,twoTrackArray);
 930     //    AddDaughterRefs(secVert,(AliAODEvent*)evento,twoTrackArray);
 931
 932   }
 933
 934   return the2Prong;
 935
 936   delete the2Prong;
 937 }
 938
 939 //----------------------------------------------------------------------------
 940 void AliAODNuclExReplicator::AddDaughterRefs(AliAODVertex *v,
 941                                             const AliAODEvent &event,
 942                                             const TObjArray *trkArray) const
 943
 944 {
 945   // Add the AOD tracks as daughters of the vertex (TRef)
 946   // AliCodeTimerAuto("",0);
 947   // cout<<"Inside  AddDaughterRefs "<<endl;
 948
 949   Int_t nDg = v->GetNDaughters();
 950
 951   TObject *dg = 0;
 952   if(nDg) dg = v->GetDaughter(0);
 953   if(dg) return; // daughters already added
 954
 955   Int_t nTrks = trkArray->GetEntriesFast();
 956
 957   AliExternalTrackParam *track = 0;
 958   AliAODTrack *aodTrack = 0;
 959   Int_t id;
 960
 961   for(Int_t i=0; i<nTrks; i++) {
 962     track = (AliExternalTrackParam*)trkArray->UncheckedAt(i);
 963
 964     id = (Int_t)track->GetID();
 965     //    printf("---> %d\n",id);
 966     if(id<0) continue; // this track is a AliAODRecoDecay
 967
 968     aodTrack = (AliAODTrack*)event.GetTrack(fAODMap[id]);
 969     v->AddDaughter(aodTrack);
 970   }
 971   return;
 972 }
 973 //----------------------------------------------------------------------------
 974
 975 void AliAODNuclExReplicator::AddRefs(AliAODVertex *v,AliAODRecoDecayLF *rd,
 976                                     const AliAODEvent &event,
 977                                     const TObjArray *trkArray) const
 978
 979 {
 980   // Add the AOD tracks as daughters of the vertex (TRef)
 981   // Also add the references to the primary vertex and to the cuts
 982   //AliCodeTimerAuto("",0);
 983
 984   AddDaughterRefs(v,event,trkArray);
 985   rd->SetPrimaryVtxRef((AliAODVertex*)event.GetPrimaryVertex());
 986   return;
 987 }
 988
 989 //---------------------------------------------------------------------------
 990
 991 void AliAODNuclExReplicator::Terminate(){
 992
 993 }