]> git.uio.no Git - u/mrichter/AliRoot.git/blame - PWGLF/STRANGENESS/Hypernuclei/AliAODNuclExReplicator.cxx
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
Merge branch 'feature-movesplit'
[u/mrichter/AliRoot.git] / PWGLF / STRANGENESS / Hypernuclei / AliAODNuclExReplicator.cxx
CommitLineData
46f95026 1/*************************************************************************
cd469d61 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
420938a5 22// AliAODRecoDecayLF2Prong and their daughter tracks.
cd469d61 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//
420938a5 28// The secondary vertices are recreated here, as a AliAODRecoDecayLF2Prong
cd469d61 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
36class AliESDv0;
37class AliESDVertex;
38class AliAODVertex;
39class 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"
420938a5 48//#include "AliAnalysisCuts.h"
cd469d61 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"
420938a5 60#include "AliAODRecoDecayLF.h"
61#include "AliAODRecoDecayLF2Prong.h"
cd469d61 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"
420938a5 78//#include "AliAnalysisVertexingLF.h"
cd469d61 79#include "AliAnalysisManager.h"
80#include "AliAODNuclExReplicator.h"
81#include "TH1.h"
82#include "TCanvas.h"
83#include "AliInputEventHandler.h"
84
85using std::cout;
86using std::endl;
87
88ClassImp(AliAODNuclExReplicator)
89
90//_____________________________________________________________________________
91AliAODNuclExReplicator::AliAODNuclExReplicator(const char* name, const char* title,
420938a5 92 // AliAnalysisCuts* trackCut,
93 // AliAnalysisCuts* vertexCut,
cd469d61 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),
420938a5 107//fVertexCut(vertexCut),
108 fVertices(0x0),
cd469d61 109 fNuclei(0x0),
420938a5 110// fTrackCut(trackCut),
111 fSecondaryVerices(0x0),
cd469d61 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//_____________________________________________________________________________
135AliAODNuclExReplicator::~AliAODNuclExReplicator()
136{
137 // destructor
420938a5 138 // delete fTrackCut;
139 // delete fVertexCut;
cd469d61 140 delete fList;
141}
142
143//_____________________________________________________________________________
144void 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//_____________________________________________________________________________
157Bool_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//_____________________________________________________________________________
167void 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//_____________________________________________________________________________
193Int_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//_____________________________________________________________________________
373TList* 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
420938a5 385 fSecondaryVerices = new TClonesArray("AliAODRecoDecayLF2Prong",30);
cd469d61 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//_____________________________________________________________________________
421void 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
7b08784f 435 // printf("Execute NuclEx Replicator\n");
cd469d61 436
420938a5 437 //---------------------------------
438
cd469d61 439 if (fReplicateHeader)
440 {
0a918d8d 441 AliAODHeader * header = dynamic_cast<AliAODHeader*>(source.GetHeader());
442 if(!header) AliFatal("Not a standard AOD");
443 *fHeader = *(header);
cd469d61 444 }
445
446 fVertices->Clear("C");
447
448 fNuclei->Clear("C");
449
450 fSecondaryVerices->Clear("C");
451
452 fDaughterTracks->Clear("C");
420938a5 453
454 //----------------------------------
cd469d61 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
420938a5 465 AliAODRecoDecayLF2Prong *io2Prong = 0;
cd469d61 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
420938a5 509 //AliAODRecoDecayLF *rd = 0;
46f95026 510 // AliAODRecoDecayLF2Prong*rd = 0;
cd469d61 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
7e4263f5 524 Double_t xPrimaryVertex=0.,yPrimaryVertex=0.;
cd469d61 525 xPrimaryVertex=vtx->GetX();
526 yPrimaryVertex=vtx->GetY();
cd469d61 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
f26537d8 568 // Double_t mom = aodtrack->GetDetPid()->GetTPCmomentum();
569 Double_t mom = aodtrack->P();
cd469d61 570
571 if(mom<0.150)continue;
f26537d8 572
573
574
cd469d61 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){
f26537d8 587
588 //Double_t triggerDeDx = 4*AliExternalTrackParam::BetheBlochAleph((mom*2)/(0.938*3),1.0288,31.9806,5.04114e-11,2.13096,2.38541);
cd469d61 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
46f95026 674
f26537d8 675 // AliAODTrack *trk0 = (AliAODTrack*)io2Prong->GetDaughter(0);
676 // AliAODTrack *trk1 = (AliAODTrack*)io2Prong->GetDaughter(1);
46f95026 677
1c4d3bfa 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;
46f95026 684
685 // rd = new((*fSecondaryVerices)[nsv++]) AliAODRecoDecayLF2Prong(*io2Prong);
686
687 new((*fSecondaryVerices)[nsv++]) AliAODRecoDecayLF2Prong(*io2Prong);
688
7b08784f 689 // cout<<"QUELLO CHE SALVo \npr0: "<<rd->GetProngID(0)<<" pr1: "<<rd->GetProngID(1)<<" pr2"<<rd->GetProngID(2)<<endl;
cd469d61 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 //----------------------------------------------------------
f26537d8 713
cd469d61 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 {
f26537d8 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 {
cd469d61 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
7b08784f 744 // printf("....Done NuclEx Replicator...\n");
cd469d61 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; }
7e4263f5 770 delete []indices;
cd469d61 771 //cout<<"Delete 6"<< endl;
772 delete fVertexerTracks;
773
774
775}
776
777//-----------------------------------------------------------------------------
778
779AliAODVertex *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
e690d4d0 797 Double_t vertRadius2=vertexESD->GetX()*vertexESD->GetX()+vertexESD->GetY()*vertexESD->GetY();
cd469d61 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
420938a5 843AliAODRecoDecayLF2Prong* AliAODNuclExReplicator::Make2Prong(TObjArray *twoTrackArray,const AliAODEvent &evento,
cd469d61 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
420938a5 855 // From AliAnalysisVertexingLF.cxx Method:Make2Prongs
cd469d61 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
7e4263f5 898 d0z0[0] = -999.;
899 d0z0[1] = -999.;
900 covd0z0[0] = -999.;
901 covd0z0[1] = -999.;
902 covd0z0[2] = -999.;
903
cd469d61 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
420938a5 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);
cd469d61 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//----------------------------------------------------------------------------
940void 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
420938a5 975void AliAODNuclExReplicator::AddRefs(AliAODVertex *v,AliAODRecoDecayLF *rd,
cd469d61 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
991void AliAODNuclExReplicator::Terminate(){
992
993}