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